diff --git a/source/module_base/lapack_connector.h b/source/module_base/lapack_connector.h index 96a7a782456..aa965bb4530 100644 --- a/source/module_base/lapack_connector.h +++ b/source/module_base/lapack_connector.h @@ -159,6 +159,20 @@ class LapackConnector return aux; } + static inline + std::complex* transpose(const std::complex* a, const int n, const int lda, const int nbase_x) + { + std::complex* aux = new std::complex[lda*n]; + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < lda; ++j) + { + aux[j * n + i] = a[i * nbase_x + j]; + } + } + return aux; + } + // Transpose the fortran-form real-std::complex array to the std::complex matrix. static inline void transpose(const std::complex* aux, ModuleBase::ComplexMatrix& a, const int n, const int lda) @@ -224,6 +238,31 @@ class LapackConnector delete[] aux; delete[] bux; } + // wrap function of fortran lapack routine zhegv. + static inline + void zhegv( const int itype, const char jobz, const char uplo, const int n, std::complex* a, + const int lda, const std::complex* b, const int ldb, double* w, std::complex* work, + int lwork, double* rwork, int info, int ld_real) + { // Transpose the std::complex matrix to the fortran-form real-std::complex array. + std::complex* aux = LapackConnector::transpose(a, n, lda, ld_real); + std::complex* bux = LapackConnector::transpose(b, n, ldb, ld_real); + + // call the fortran routine + zhegv_(&itype, &jobz, &uplo, &n, aux, &lda, bux, &ldb, w, work, &lwork, rwork, &info); + // Transpose the fortran-form real-std::complex array to the std::complex matrix. + // LapackConnector::transpose(aux, a, n, lda); + // LapackConnector::transpose(bux, b, n, ldb); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < lda; ++j) + { + a[j * ld_real + i] = aux[i*lda+j]; + } + } + // free the memory. + delete[] aux; + delete[] bux; + } // wrap function of fortran lapack routine zheev. static inline @@ -252,6 +291,38 @@ class LapackConnector delete[] bux; delete[] zux; + } + // wrap function of fortran lapack routine zheev. + static inline + void zhegvx( const int itype, const char jobz, const char range, const char uplo, + const int n, const std::complex* a, const int lda, const std::complex* b, + const int ldb, const double vl, const double vu, const int il, const int iu, + const double abstol, const int m, double* w, std::complex* z, const int ldz, + std::complex* work, const int lwork, double* rwork, int* iwork, + int* ifail, int info, int nbase_x) + { + // Transpose the std::complex matrix to the fortran-form real-std::complex array. + std::complex* aux = LapackConnector::transpose(a, n, lda, nbase_x); + std::complex* bux = LapackConnector::transpose(b, n, ldb, nbase_x); + std::complex* zux = new std::complex[n*iu];// mohan modify 2009-08-02 + + // call the fortran routine + zhegvx_(&itype, &jobz, &range, &uplo, &n, aux, &lda, bux, &ldb, &vl, + &vu, &il,&iu, &abstol, &m, w, zux, &ldz, work, &lwork, rwork, iwork, ifail, &info); + + // Transpose the fortran-form real-std::complex array to the std::complex matrix + for (int i = 0; i < iu; ++i) + { + for (int j = 0; j < n; ++j) + { + z[j * nbase_x + i] = zux[i*n+j]; + } + } + + // free the memory. + delete[] aux; + delete[] bux; + delete[] zux; } // calculate the eigenvalues and eigenfunctions of a real symmetric matrix. diff --git a/source/module_hsolver/CMakeLists.txt b/source/module_hsolver/CMakeLists.txt index bbc61e8a35d..ebbacc48458 100644 --- a/source/module_hsolver/CMakeLists.txt +++ b/source/module_hsolver/CMakeLists.txt @@ -5,6 +5,7 @@ list(APPEND objects hsolver_pw_sdft.cpp diago_iter_assist.cpp src/math_kernel.cpp + src/dngvd_op.cpp ) if(ENABLE_LCAO) @@ -16,7 +17,10 @@ if(ENABLE_LCAO) endif() if (USE_CUDA) - list(APPEND objects src/cuda/math_kernel.cu) + list(APPEND objects + src/cuda/math_kernel.cu + src/cuda/dngvd_op.cu + ) elseif(USE_ROCM) list(APPEND objects src/rocm/math_kernel.cu) endif() diff --git a/source/module_hsolver/diago_cg.cpp b/source/module_hsolver/diago_cg.cpp index 1f2a39698ca..cc8a96f0c02 100644 --- a/source/module_hsolver/diago_cg.cpp +++ b/source/module_hsolver/diago_cg.cpp @@ -683,23 +683,7 @@ void DiagoCG::diag(hamilt::Hamilt *phm_in, psi:: { if(DiagoIterAssist::need_subspace || ntry > 0) { - #if defined(__CUDA) || defined(__ROCM) - psi::Psi, psi::DEVICE_CPU> cpu_psi = psi; - // hamilt::Hamilt* h_phm_in = new hamilt::HamiltPW(reinterpret_cast*>(phm_in)); - hamilt::Hamilt* h_phm_in = - new hamilt::HamiltPW( - reinterpret_cast*>(phm_in)); - DiagoIterAssist::diagH_subspace(h_phm_in, cpu_psi, cpu_psi, eigenvalue_in); - psi::memory::synchronize_memory_op, Device, psi::DEVICE_CPU>()( - psi.get_device(), - cpu_psi.get_device(), - psi.get_pointer() - psi.get_psi_bias(), - cpu_psi.get_pointer() - cpu_psi.get_psi_bias(), - psi.size()); - delete reinterpret_cast*>(h_phm_in); - #else - DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_in); - #endif + DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_in); } DiagoIterAssist::avg_iter += 1.0; diff --git a/source/module_hsolver/diago_iter_assist.cpp b/source/module_hsolver/diago_iter_assist.cpp index be44920b1a1..d122386bf22 100644 --- a/source/module_hsolver/diago_iter_assist.cpp +++ b/source/module_hsolver/diago_iter_assist.cpp @@ -7,6 +7,8 @@ #include "module_base/lapack_connector.h" #include "module_base/timer.h" #include "src_parallel/parallel_reduce.h" +#include "module_hsolver/include/math_kernel.h" +#include "module_hsolver/include/dngvd_op.h" using namespace hsolver; @@ -46,67 +48,94 @@ void DiagoIterAssist::diagH_subspace( n_band = nstart; assert(n_band <= nstart); - ModuleBase::ComplexMatrix hc(nstart, nstart); - ModuleBase::ComplexMatrix sc(nstart, nstart); - ModuleBase::ComplexMatrix hvec(nstart, n_band); + Device* ctx = {}; + + std::complex* hcc = nullptr; + std::complex* scc = nullptr; + std::complex* vcc = nullptr; + psi::memory::resize_memory_op, Device>()(ctx, hcc, nstart * nstart); + psi::memory::resize_memory_op, Device>()(ctx, scc, nstart * nstart); + psi::memory::resize_memory_op, Device>()(ctx, vcc, nstart * nstart); + psi::memory::set_memory_op, Device>()(ctx, hcc, 0, nstart * nstart); + psi::memory::set_memory_op, Device>()(ctx, scc, 0, nstart * nstart); + psi::memory::set_memory_op, Device>()(ctx, vcc, 0, nstart * nstart); const int dmin = psi.get_current_nbas(); const int dmax = psi.get_nbasis(); // qianrui improve this part 2021-3-14 - //std::complex *aux = new std::complex[dmax * nstart]; - //const std::complex *paux = aux; const std::complex* ppsi = psi.get_pointer(); - //allocated hpsi - std::vector> hpsi(psi.get_nbands() * psi.get_nbasis()); - //do hPsi for all bands + // allocated hpsi + // std::vector> hpsi(psi.get_nbands() * psi.get_nbasis()); + std::complex* hphi = nullptr; + psi::memory::resize_memory_op, Device>()(ctx, hphi, psi.get_nbands() * psi.get_nbasis()); + psi::memory::set_memory_op, Device>()(ctx, hphi, 0, psi.get_nbands() * psi.get_nbasis()); + // do hPsi for all bands psi::Range all_bands_range(1, psi.get_current_k(), 0, psi.get_nbands()-1); - hpsi_info hpsi_in(&psi, all_bands_range, hpsi.data()); + hpsi_info hpsi_in(&psi, all_bands_range, hphi); pHamilt->ops->hPsi(hpsi_in); - //use aux as a data pointer for hpsi - const std::complex *aux = hpsi.data(); - char trans1 = 'C'; - char trans2 = 'N'; - zgemm_(&trans1, - &trans2, - &nstart, - &nstart, - &dmin, - &ModuleBase::ONE, - ppsi, - &dmax, - aux, - &dmax, - &ModuleBase::ZERO, - hc.c, - &nstart); - hc = transpose(hc, false); - - zgemm_(&trans1, - &trans2, - &nstart, - &nstart, - &dmin, - &ModuleBase::ONE, - ppsi, - &dmax, - ppsi, - &dmax, - &ModuleBase::ZERO, - sc.c, - &nstart); - sc = transpose(sc, false); + gemm_op()( + ctx, + 'C', + 'N', + nstart, + nstart, + dmin, + &ModuleBase::ONE, + ppsi, + dmax, + hphi, + dmax, + &ModuleBase::ZERO, + hcc, + nstart + ); + matrixTranspose_op()( + ctx, + nstart, + nstart, + hcc, + hcc + ); + + gemm_op()( + ctx, + 'C', + 'N', + nstart, + nstart, + dmin, + &ModuleBase::ONE, + ppsi, + dmax, + ppsi, + dmax, + &ModuleBase::ZERO, + scc, + nstart + ); + matrixTranspose_op()( + ctx, + nstart, + nstart, + scc, + scc + ); if (GlobalV::NPROC_IN_POOL > 1) { - Parallel_Reduce::reduce_complex_double_pool(hc.c, nstart * nstart); - Parallel_Reduce::reduce_complex_double_pool(sc.c, nstart * nstart); + Parallel_Reduce::reduce_complex_double_pool(hcc, nstart * nstart); + Parallel_Reduce::reduce_complex_double_pool(scc, nstart * nstart); } // after generation of H and S matrix, diag them - DiagoIterAssist::diagH_LAPACK(nstart, n_band, hc, sc, nstart, en, hvec); + DiagoIterAssist::diagH_LAPACK(nstart, n_band, hcc, scc, nstart, en, vcc); + +#if defined(__CUDA) || defined(__ROCM) + hsolver::createBLAShandle(); +#endif //======================= // diagonize the H-matrix @@ -122,52 +151,65 @@ void DiagoIterAssist::diagH_subspace( // because psi and evc are different here, // I think if psi and evc are the same, // there may be problems, mohan 2011-01-01 - char transa = 'N'; - char transb = 'T'; - zgemm_(&transa, - &transb, - &dmax, // m: row of A,C - &n_band, // n: col of B,C - &nstart, // k: col of A, row of B - &ModuleBase::ONE, // alpha - ppsi, // A - &dmax, // LDA: if(N) max(1,m) if(T) max(1,k) - hvec.c, // B - &n_band, // LDB: if(N) max(1,k) if(T) max(1,n) - &ModuleBase::ZERO, // belta - evc.get_pointer(), // C - &dmax); // LDC: if(N) max(1, m) + gemm_op()( + ctx, + 'N', + 'T', + dmax, + n_band, + nstart, + &ModuleBase::ONE, + ppsi, + dmax, + vcc, + n_band, + &ModuleBase::ZERO, + evc.get_pointer(), + dmax + ); } else { // As the evc and psi may refer to the same matrix, we first // create a temporary matrix to store the result. (by wangjp) // qianrui improve this part 2021-3-13 - char transa = 'N'; - char transb = 'T'; - ModuleBase::ComplexMatrix evctmp(n_band, dmin, false); - zgemm_(&transa, - &transb, - &dmin, - &n_band, - &nstart, - &ModuleBase::ONE, - ppsi, - &dmax, - hvec.c, - &n_band, - &ModuleBase::ZERO, - evctmp.c, - &dmin); - for (int ib = 0; ib < n_band; ib++) - { - for (int ig = 0; ig < dmin; ig++) - { - evc(ib, ig) = evctmp(ib, ig); - } - } + std::complex* evctemp = nullptr; + psi::memory::resize_memory_op, Device>()(ctx, evctemp, n_band * dmin); + psi::memory::set_memory_op, Device>()(ctx, evctemp, 0, n_band * dmin); + + gemm_op()( + ctx, + 'N', + 'T', + dmin, + n_band, + nstart, + &ModuleBase::ONE, + ppsi, + dmax, + vcc, + n_band, + &ModuleBase::ZERO, + evctemp, + dmin + ); + + matrixSetToAnother()(ctx, n_band, evctemp, dmin, evc.get_pointer(), dmax); + // for (int ib = 0; ib < n_band; ib++) + // { + // for (int ig = 0; ig < dmin; ig++) + // { + // evc(ib, ig) = evctmp(ib, ig); + // } + // } + psi::memory::delete_memory_op, Device>()(ctx, evctemp); } + psi::memory::delete_memory_op, Device>()(ctx, hcc); + psi::memory::delete_memory_op, Device>()(ctx, scc); + psi::memory::delete_memory_op, Device>()(ctx, vcc); + psi::memory::delete_memory_op, Device>()(ctx, hphi); + ModuleBase::timer::tick("DiagoIterAssist", "diagH_subspace"); return; } @@ -334,10 +376,10 @@ void DiagoIterAssist::diagH_subspace_init( // by nstart states psi (atomic or random wavefunctions). // Produces on output n_band eigenvectors (n_band <= nstart) in evc. //---------------------------------------------------------------------- -template<> -void DiagoIterAssist::diagH_subspace(hamilt::Hamilt* pHamilt, const psi::Psi, psi::DEVICE_GPU> &psi, psi::Psi, psi::DEVICE_GPU> &evc, double *en, int n_band) { - ModuleBase::WARNING_QUIT("DiagoIterAssist::diagH_subspace","GPU's implementation is not supported currently!"); -} +// template<> +// void DiagoIterAssist::diagH_subspace(hamilt::Hamilt* pHamilt, const psi::Psi, psi::DEVICE_GPU> &psi, psi::Psi, psi::DEVICE_GPU> &evc, double *en, int n_band) { +// ModuleBase::WARNING_QUIT("DiagoIterAssist::diagH_subspace","GPU's implementation is not supported currently!"); +// } template<> void DiagoIterAssist::diagH_subspace_init(hamilt::Hamilt* pHamilt, const ModuleBase::ComplexMatrix &psi, psi::Psi, psi::DEVICE_GPU> &evc, double *en) { @@ -469,6 +511,86 @@ void DiagoIterAssist::diagH_LAPACK( return; } + +template +void DiagoIterAssist::diagH_LAPACK( + const int nstart, + const int nbands, + const std::complex* hcc, + const std::complex* scc, + const int ldh, // nstart + FPTYPE *e, // always in CPU + std::complex* vcc) +{ + ModuleBase::TITLE("DiagoIterAssist", "LAPACK_subspace"); + ModuleBase::timer::tick("DiagoIterAssist", "LAPACK_subspace"); + + Device* ctx = {}; + + const bool all_eigenvalues = (nstart == nbands); + +#if defined(__CUDA) || defined(__ROCM) + const psi::DEVICE_CPU * cpu_ctx = {}; + const psi::DEVICE_GPU * gpu_ctx = {}; + + FPTYPE* e_gpu = nullptr; + psi::memory::resize_memory_op()(gpu_ctx, e_gpu, nbands); + // set e in CPU value to e_gpu + psi::memory::synchronize_memory_op()( + gpu_ctx, + cpu_ctx, + e_gpu, + e, + nbands + ); + + if (all_eigenvalues) + { + //=========================== + // calculate all eigenvalues + //=========================== + dngv_op()(ctx, nstart, ldh, hcc, scc, e_gpu, vcc); + } + else + { + //===================================== + // calculate only m lowest eigenvalues + //===================================== + dngvx_op()(ctx, nstart, ldh, hcc, scc, nbands, e_gpu, vcc); + } + + // set e_gpu value to e in CPU + psi::memory::synchronize_memory_op()( + cpu_ctx, + gpu_ctx, + e, + e_gpu, + nbands + ); + psi::memory::delete_memory_op()(gpu_ctx, e_gpu); +#else + + if (all_eigenvalues) + { + //=========================== + // calculate all eigenvalues + //=========================== + dngv_op()(ctx, nstart, ldh, hcc, scc, e, vcc); + } + else + { + //===================================== + // calculate only m lowest eigenvalues + //===================================== + dngvx_op()(ctx, nstart, ldh, hcc, scc, nbands, e, vcc); + } + +#endif + + ModuleBase::timer::tick("DiagoIterAssist", "LAPACK_subspace"); + return; +} + template bool DiagoIterAssist::test_exit_cond(const int &ntry, const int ¬conv) { diff --git a/source/module_hsolver/diago_iter_assist.h b/source/module_hsolver/diago_iter_assist.h index 79bc4baa30f..09692974f03 100644 --- a/source/module_hsolver/diago_iter_assist.h +++ b/source/module_hsolver/diago_iter_assist.h @@ -42,9 +42,18 @@ class DiagoIterAssist FPTYPE *e, ModuleBase::ComplexMatrix &hvec); + static void diagH_LAPACK( + const int nstart, + const int nbands, + const std::complex* hcc, + const std::complex* sc, + const int ldh, // nstart + FPTYPE *e, + std::complex* vcc); + static bool test_exit_cond(const int &ntry, const int ¬conv); - using hpsi_info = typename hamilt::Operator, psi::DEVICE_CPU>::hpsi_info; + using hpsi_info = typename hamilt::Operator, Device>::hpsi_info; }; } // namespace hsolver diff --git a/source/module_hsolver/include/dngvd_op.h b/source/module_hsolver/include/dngvd_op.h new file mode 100644 index 00000000000..2eb6643388b --- /dev/null +++ b/source/module_hsolver/include/dngvd_op.h @@ -0,0 +1,63 @@ +// TODO: This is a temperary location for these functions. +// And will be moved to a global module(module base) later. +#ifndef MODULE_HSOLVER_DNGVD_H +#define MODULE_HSOLVER_DNGVD_H + +#include "module_base/lapack_connector.h" +#include "module_hsolver/include/math_kernel.h" +#include "module_psi/include/memory.h" + +namespace hsolver +{ + +template struct dngvx_op +{ + /// @brief DNGVX computes the first m eigenvalues ​​and their corresponding eigenvectors of + /// a complex generalized Hermitian-definite eigenproblem + /// + /// Input Parameters + /// @param d : the type of device + /// @param row : the number of rows of the matrix + /// @param col : the number of cols of the matrix + /// @param A : the hermitian matrix A in A x=lambda B x (row major) + /// @param B : the overlap matrix B in A x=lambda B x (row major) + /// @param m : the number of the first m eigenvalues to calculate + /// Output Parameter + /// @param W : calculated eigenvalues + /// @param V : calculated eigenvectors (row major) + void operator()(const Device* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + const int m, + double* W, + std::complex* V); +}; + +template struct dngv_op +{ + /// @brief DNGV computes all the eigenvalues and eigenvectors of a complex generalized + /// Hermitian-definite eigenproblem + /// + /// Input Parameters + /// @param d : the type of device + /// @param row : the number of rows of the matrix + /// @param col : the number of cols of the matrix + /// @param A : the hermitian matrix A in A x=lambda B x (row major) + /// @param B : the overlap matrix B in A x=lambda B x (row major) + /// Output Parameter + /// @param W : calculated eigenvalues + /// @param V : calculated eigenvectors (row major) + void operator()(const Device* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + double* W, + std::complex* V); +}; + +} // namespace hsolver + +#endif // !MODULE_HSOLVER_DNGVD_H \ No newline at end of file diff --git a/source/module_hsolver/include/math_kernel.h b/source/module_hsolver/include/math_kernel.h index d343a0f3dcf..1190c043018 100644 --- a/source/module_hsolver/include/math_kernel.h +++ b/source/module_hsolver/include/math_kernel.h @@ -6,15 +6,16 @@ #include "module_base/blas_connector.h" #include "module_psi/psi.h" #include "src_parallel/parallel_reduce.h" +#include "module_psi/include/memory.h" #if defined(__CUDA) || defined(__UT_USE_CUDA) #include #include "cublas_v2.h" -#define cublasErrcheck(res) { cudaAssert((res), __FILE__, __LINE__); } +#define cublasErrcheck(res) { cublasAssert((res), __FILE__, __LINE__); } -static const char *_cudaGetErrorEnum(cublasStatus_t error) { +static const char *_cublasGetErrorEnum(cublasStatus_t error) { switch (error) { case CUBLAS_STATUS_SUCCESS: return "CUBLAS_STATUS_SUCCESS"; @@ -43,9 +44,9 @@ static const char *_cudaGetErrorEnum(cublasStatus_t error) { return ""; } -inline void cudaAssert(cublasStatus_t code, const char *file, int line, bool abort=true) { +inline void cublasAssert(cublasStatus_t code, const char *file, int line, bool abort=true) { if (code != CUBLAS_STATUS_SUCCESS) { - fprintf(stderr,"cuBLAS Assert: %s %s %d\n", _cudaGetErrorEnum(code), file, line); + fprintf(stderr,"cuBLAS Assert: %s %s %d\n", _cublasGetErrorEnum(code), file, line); if (abort) exit(code); } } @@ -165,6 +166,25 @@ template struct gemm_op const int& ldc); }; +template struct matrixTranspose_op +{ + void operator()(const Device* d, + const int& row, + const int& col, + const std::complex* input_matrix, + std::complex* output_matrix); +}; + +template struct matrixSetToAnother +{ + void operator()(const Device* d, + const int& n, + const std::complex* A, + const int& LDA, + std::complex* B, + const int& LDB); +}; + #if __CUDA || __UT_USE_CUDA || __ROCM || __UT_USE_ROCM // Partially specialize functor for psi::GpuDevice. @@ -220,6 +240,17 @@ template struct constantvector_addORsub_constantVector_op struct matrixSetToAnother +{ + void operator()(const psi::DEVICE_GPU* d, + const int& n, + const std::complex* A, + const int& LDA, + std::complex* B, + const int& LDB); +}; + + void createBLAShandle(); void destoryBLAShandle(); diff --git a/source/module_hsolver/src/cuda/dngvd_op.cu b/source/module_hsolver/src/cuda/dngvd_op.cu new file mode 100644 index 00000000000..e67cb4cef44 --- /dev/null +++ b/source/module_hsolver/src/cuda/dngvd_op.cu @@ -0,0 +1,235 @@ +#include "module_hsolver/include/dngvd_op.h" +#include "src_pdiag/helper_cuda.h" + +#include + +#define cusolverErrcheck(res) \ + { \ + cusolverAssert((res), __FILE__, __LINE__); \ + } + +// cuSOLVER API errors +static const char* _cusolverGetErrorEnum(cusolverStatus_t error) +{ + switch (error) + { + case CUSOLVER_STATUS_SUCCESS: + return "CUSOLVER_STATUS_SUCCESS"; + case CUSOLVER_STATUS_NOT_INITIALIZED: + return "CUSOLVER_STATUS_NOT_INITIALIZED"; + case CUSOLVER_STATUS_ALLOC_FAILED: + return "CUSOLVER_STATUS_ALLOC_FAILED"; + case CUSOLVER_STATUS_INVALID_VALUE: + return "CUSOLVER_STATUS_INVALID_VALUE"; + case CUSOLVER_STATUS_ARCH_MISMATCH: + return "CUSOLVER_STATUS_ARCH_MISMATCH"; + case CUSOLVER_STATUS_MAPPING_ERROR: + return "CUSOLVER_STATUS_MAPPING_ERROR"; + case CUSOLVER_STATUS_EXECUTION_FAILED: + return "CUSOLVER_STATUS_EXECUTION_FAILED"; + case CUSOLVER_STATUS_INTERNAL_ERROR: + return "CUSOLVER_STATUS_INTERNAL_ERROR"; + case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED: + return "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED"; + case CUSOLVER_STATUS_NOT_SUPPORTED: + return "CUSOLVER_STATUS_NOT_SUPPORTED "; + case CUSOLVER_STATUS_ZERO_PIVOT: + return "CUSOLVER_STATUS_ZERO_PIVOT"; + case CUSOLVER_STATUS_INVALID_LICENSE: + return "CUSOLVER_STATUS_INVALID_LICENSE"; + } + return ""; +} + +inline void cusolverAssert(cusolverStatus_t code, const char* file, int line, bool abort = true) +{ + if (code != CUBLAS_STATUS_SUCCESS) + { + fprintf(stderr, "cuSOLVER Assert: %s %s %d\n", _cusolverGetErrorEnum(code), file, line); + if (abort) + exit(code); + } +} + +namespace hsolver +{ + +template <> +void dngvx_op::operator()(const psi::DEVICE_GPU* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + const int m, + double* W, + std::complex* V) +{ + createBLAShandle(); + + // init A_eigenvectors, transpose_B and all_W + double2 *A_eigenvectors, *transpose_B; + checkCudaErrors(cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col)); + checkCudaErrors(cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col)); + + matrixTranspose_op()(d, col, row, A, (std::complex*)A_eigenvectors); + matrixTranspose_op()(d, col, row, B, (std::complex*)transpose_B); + + double* all_W; + checkCudaErrors(cudaMalloc((void**)&all_W, sizeof(double) * col)); + + // prepare some values for cusolverDnZhegvd_bufferSize + cusolverDnHandle_t cusolverH; + cusolverErrcheck(cusolverDnCreate(&cusolverH)); + int* devInfo; + checkCudaErrors(cudaMalloc((void**)&devInfo, sizeof(int))); + + // calculate the sizes needed for pre-allocated buffer. + int lwork = 0; + cusolverErrcheck(cusolverDnZhegvd_bufferSize( + cusolverH, + CUSOLVER_EIG_TYPE_1, // itype = CUSOLVER_EIG_TYPE_1: A*x = (lambda)*B*x. + CUSOLVER_EIG_MODE_VECTOR, // jobz = CUSOLVER_EIG_MODE_VECTOR : Compute eigenvalues and eigenvectors. + CUBLAS_FILL_MODE_LOWER, + col, + A_eigenvectors, + // (double2)A, + row, + transpose_B, + // (double2)B, + row, + all_W, + &lwork)); + + // allocate memery + cuDoubleComplex* d_work; + checkCudaErrors(cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork)); + + // compute eigenvalues and eigenvectors. + cusolverErrcheck(cusolverDnZhegvd( + cusolverH, + CUSOLVER_EIG_TYPE_1, // itype = CUSOLVER_EIG_TYPE_1: A*x = (lambda)*B*x. + CUSOLVER_EIG_MODE_VECTOR, // jobz = CUSOLVER_EIG_MODE_VECTOR : Compute eigenvalues and eigenvectors. + CUBLAS_FILL_MODE_LOWER, + col, + A_eigenvectors, + // (double2)A, + row, + transpose_B, + // (double2)B, + row, + all_W, + d_work, + lwork, + devInfo)); + + checkCudaErrors(cudaDeviceSynchronize()); + + // get eigenvalues and eigenvectors. only m ! + checkCudaErrors(cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice)); + + checkCudaErrors(cudaMemcpy(V, A_eigenvectors, sizeof(std::complex) * col * m, cudaMemcpyDeviceToDevice)); + + matrixTranspose_op()(d, col, row, V, V); + + int info_gpu; + checkCudaErrors(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost)); + assert(0 == info_gpu); + + // free the buffer + checkCudaErrors(cudaFree(d_work)); + // free resources and destroy + checkCudaErrors(cudaFree(A_eigenvectors)); + checkCudaErrors(cudaFree(all_W)); + checkCudaErrors(cudaFree(devInfo)); + cusolverErrcheck(cusolverDnDestroy(cusolverH)); + + destoryBLAShandle(); +} + +template <> +void dngv_op::operator()(const psi::DEVICE_GPU* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + double* W, + std::complex* V) +{ + createBLAShandle(); + + // init A_eigenvectors & transpose_B + double2 *A_eigenvectors, *transpose_B; + checkCudaErrors(cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col)); + checkCudaErrors(cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col)); + + // transpose A, B to A_eigenvectors, transpose_B + matrixTranspose_op()(d, row, col, A, (std::complex*)A_eigenvectors); + matrixTranspose_op()(d, row, col, B, (std::complex*)transpose_B); + + // init all_W + double* all_W; + checkCudaErrors(cudaMalloc((void**)&all_W, sizeof(double) * row)); + + // prepare some values for cusolverDnZhegvd_bufferSize + cusolverDnHandle_t cusolverH; + cusolverErrcheck(cusolverDnCreate(&cusolverH)); + int* devInfo; + checkCudaErrors(cudaMalloc((void**)&devInfo, sizeof(int))); + + // calculate the sizes needed for pre-allocated buffer. + int lwork = 0; + cusolverErrcheck(cusolverDnZhegvd_bufferSize( + cusolverH, + CUSOLVER_EIG_TYPE_1, // itype = CUSOLVER_EIG_TYPE_1: A*x = (lambda)*B*x. + CUSOLVER_EIG_MODE_VECTOR, // jobz = CUSOLVER_EIG_MODE_VECTOR : Compute eigenvalues and eigenvectors. + CUBLAS_FILL_MODE_UPPER, + row, + A_eigenvectors, + col, + transpose_B, + col, + all_W, + &lwork)); + + // allocate memery + cuDoubleComplex* d_work; + checkCudaErrors(cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork)); + + // compute eigenvalues and eigenvectors. + cusolverErrcheck(cusolverDnZhegvd( + cusolverH, + CUSOLVER_EIG_TYPE_1, // itype = CUSOLVER_EIG_TYPE_1: A*x = (lambda)*B*x. + CUSOLVER_EIG_MODE_VECTOR, // jobz = CUSOLVER_EIG_MODE_VECTOR : Compute eigenvalues and eigenvectors. + CUBLAS_FILL_MODE_UPPER, + row, + A_eigenvectors, + col, + transpose_B, + col, + all_W, + d_work, + lwork, + devInfo)); + + checkCudaErrors(cudaDeviceSynchronize()); + + // get all eigenvalues and eigenvectors. + checkCudaErrors(cudaMemcpy(W, all_W, sizeof(double) * row, cudaMemcpyDeviceToDevice)); + matrixTranspose_op()(d, row, col, (std::complex*)A_eigenvectors, V); + + int info_gpu; + checkCudaErrors(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost)); + assert(0 == info_gpu); + + // free the buffer + checkCudaErrors(cudaFree(d_work)); + // free resources and destroy + checkCudaErrors(cudaFree(A_eigenvectors)); + checkCudaErrors(cudaFree(all_W)); + checkCudaErrors(cudaFree(devInfo)); + cusolverErrcheck(cusolverDnDestroy(cusolverH)); + + destoryBLAShandle(); +} + +} // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/src/cuda/math_kernel.cu b/source/module_hsolver/src/cuda/math_kernel.cu index 3f8c04214c9..bb39ac62428 100644 --- a/source/module_hsolver/src/cuda/math_kernel.cu +++ b/source/module_hsolver/src/cuda/math_kernel.cu @@ -1,5 +1,6 @@ #include "module_hsolver/include/math_kernel.h" #include "module_psi/psi.h" +#include "module_psi/include/memory.h" #include #include @@ -368,6 +369,88 @@ void gemm_op::operator()(const psi::DEVICE_GPU* d, } +template +__global__ void matrix_transpose_kernel( + const int row, + const int col, + const thrust::complex* in, + thrust::complex* out) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < row) + { + for (int j = 0; j < col; j++) + { + out[j * row + i] = in[i * col + j]; + } + } +} + +template <> +void matrixTranspose_op::operator()(const psi::DEVICE_GPU* d, + const int& row, + const int& col, + const std::complex* input_matrix, + std::complex* output_matrix) +{ + std::complex* device_temp = nullptr; + psi::memory::resize_memory_op, psi::DEVICE_GPU>()(d, device_temp, row * col); + + if (row == col) + { + double2 ONE, ZERO; + ONE.x = 1.0; + ONE.y = 0.0; + ZERO.x = ZERO.y = 0.0; + + // use 'geam' API todo transpose. + cublasErrcheck(cublasZgeam(diag_handle, CUBLAS_OP_T, CUBLAS_OP_N, col, row, &ONE, (double2*)input_matrix, col, &ZERO, (double2*)input_matrix, col, (double2*)device_temp, col)); + } else + { + int thread = 1024; + int block = (row + col + thread - 1) / thread; + matrix_transpose_kernel<<>>(row, col, (thrust::complex*)input_matrix, (thrust::complex*)device_temp); + } + + psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_GPU>()(d, d, output_matrix, device_temp, row * col); + +} + + +template +__global__ void matrix_setTo_another_kernel( + const int n, + const int LDA, + const int LDB, + const thrust::complex* matrix_A, + thrust::complex* matrix_B) +{ + int j = blockIdx.x * blockDim.x + threadIdx.x; + if (j < LDA) + { + for (int i = 0; i < n; i++) + { + matrix_B[i * LDB + j] = matrix_A[i * LDA + j]; + } + } +} + + +template +void matrixSetToAnother::operator()( + const psi::DEVICE_GPU* d, + const int& n, + const std::complex* A, + const int& LDA, + std::complex* B, + const int& LDB) +{ + int thread = 1024; + int block = (LDA + thread - 1) / thread; + matrix_setTo_another_kernel<<>>(n, LDA, LDB, (thrust::complex*)A, (thrust::complex*)B); +} + + // Explicitly instantiate functors for the types of functor registered. template struct zdot_real_op; @@ -376,5 +459,6 @@ template struct vector_mul_vector_op; template struct vector_div_vector_op; template struct constantvector_addORsub_constantVector_op; +template struct matrixSetToAnother; } diff --git a/source/module_hsolver/src/dngvd_op.cpp b/source/module_hsolver/src/dngvd_op.cpp new file mode 100644 index 00000000000..62ef44c326b --- /dev/null +++ b/source/module_hsolver/src/dngvd_op.cpp @@ -0,0 +1,205 @@ +#include "module_hsolver/include/dngvd_op.h" + +#include + +namespace hsolver +{ + +template <> +void dngvx_op::operator()(const psi::DEVICE_CPU* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + const int m, + double* W, + std::complex* V) +{ + int lwork; + int info = 0; + + std::string name1 = "ZHETRD"; + std::string name2 = "L"; + + int nb = LapackConnector::ilaenv(1, name1.c_str(), name2.c_str(), col, -1, -1, -1); + if (nb < 1) + { + nb = std::max(1, col); + } + + if (nb == 1 || nb >= col) + { + lwork = 2 * col; // qianrui fix a bug 2021-7-25 : lwork should be at least max(1,2*n) + } + else + { + lwork = (nb + 1) * col; + } + + std::complex* work = new std::complex[2 * lwork]; + assert(work != 0); + double* rwork = new double[7 * col]; + assert(rwork != 0); + int* iwork = new int[5 * col]; + assert(iwork != 0); + int* ifail = new int[col]; + assert(ifail != 0); + ModuleBase::GlobalFunc::ZEROS(work, lwork); // qianrui change it, only first lwork numbers are used in zhegvx + ModuleBase::GlobalFunc::ZEROS(rwork, 7 * col); + ModuleBase::GlobalFunc::ZEROS(iwork, 5 * col); + ModuleBase::GlobalFunc::ZEROS(ifail, col); + + LapackConnector::zhegvx( + 1, // ITYPE = 1: A*x = (lambda)*B*x + 'V', // JOBZ = 'V': Compute eigenvalues and eigenvectors. + 'I', // RANGE = 'I': the IL-th through IU-th eigenvalues will be found. + 'L', // UPLO = 'L': Lower triangles of A and B are stored. + col, // N = base + A, // A is COMPLEX*16 array dimension (LDA, N) + col, // LDA = base + B, // B is COMPLEX*16 array, dimension (LDB, N) + col, // LDB = base + 0.0, // Not referenced if RANGE = 'A' or 'I'. + 0.0, // Not referenced if RANGE = 'A' or 'I'. + 1, // IL: If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, + m, // IU: If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, + 0.0, // ABSTOL + m, // M: The total number of eigenvalues found. 0 <= M <= N. if RANGE = 'I', M = IU-IL+1. + W, // W store eigenvalues + V, // store eigenvector + col, // LDZ: The leading dimension of the array Z. + work, + lwork, + rwork, + iwork, + ifail, + info, + row); + + assert(0 == info); + + delete[] work; + delete[] rwork; + delete[] iwork; + delete[] ifail; +}; + +template <> +void dngv_op::operator()(const psi::DEVICE_CPU* d, + const int row, + const int col, + const std::complex* A, + const std::complex* B, + double* W, + std::complex* V) +{ + int lwork = 0; + int nb = LapackConnector::ilaenv(1, "ZHETRD", "U", col, -1, -1, -1); + if (nb < 1) + { + nb = std::max(1, col); + } + if (nb == 1 || nb >= col) + { + lwork = 2 * col; // mohan modify 2009-08-02 + } + else + { + lwork = (nb + 1) * col; + } + + std::complex* work = new std::complex[lwork]; + ModuleBase::GlobalFunc::ZEROS(work, lwork); + + //===================================================================== + // input s and (see below) h are copied so that they are not destroyed + //===================================================================== + + int info = 0; + int rwork_dim; + rwork_dim = 3 * col - 2; + double* rwork = new double[rwork_dim]; + ModuleBase::GlobalFunc::ZEROS(rwork, rwork_dim); + + for (int i = 0; i < row * col; i++) + { + V[i] = A[i]; + } + LapackConnector::zhegv(1, 'V', 'U', col, V, col, B, col, W, work, lwork, rwork, info, row); + + assert(0 == info); + + delete[] rwork; + delete[] work; +} + +// template <> +// void dngvx_op::operator()( +// const psi::DEVICE_CPU* d, +// const int row, +// const int col, +// const std::complex* A, +// const std::complex* B, +// const int m, +// float* W, +// std::complex* V) +// { +// int lwork; +// int info = 0; + +// std::string name1 = "ZHETRD"; +// std::string name2 = "L"; + +// int nb = LapackConnector::ilaenv(1, name1.c_str(), name2.c_str(), col, -1, -1, -1); +// if (nb < 1) +// { +// nb = std::max(1, col); +// } + +// if (nb == 1 || nb >= col) +// { +// lwork = 2 * col; // qianrui fix a bug 2021-7-25 : lwork should be at least max(1,2*n) +// } else +// { +// lwork = (nb + 1) * col; +// } + +// std::complex *work = new std::complex[2 * lwork]; +// assert(work != 0); +// float *rwork = new float[7 * col]; +// assert(rwork != 0); +// int *iwork = new int[5 * col]; +// assert(iwork != 0); +// int *ifail = new int[col]; +// assert(ifail != 0); +// ModuleBase::GlobalFunc::ZEROS(work, lwork); // qianrui change it, only first lwork numbers are used in zhegvx +// ModuleBase::GlobalFunc::ZEROS(rwork, 7 * col); +// ModuleBase::GlobalFunc::ZEROS(iwork, 5 * col); +// ModuleBase::GlobalFunc::ZEROS(ifail, col); + +// LapackConnector::chegvx(1, // ITYPE = 1: A*x = (lambda)*B*x +// 'V', // JOBZ = 'V': Compute eigenvalues and eigenvectors. +// 'I', // RANGE = 'I': the IL-th through IU-th eigenvalues will be found. +// 'L', // UPLO = 'L': Lower triangles of A and B are stored. +// col, // N = base +// A, // A is COMPLEX*16 array dimension (LDA, N) +// col, // LDA = base +// B, // B is COMPLEX*16 array, dimension (LDB, N) +// col, // LDB = base +// 0.0, // Not referenced if RANGE = 'A' or 'I'. +// 0.0, // Not referenced if RANGE = 'A' or 'I'. +// 1, // IL: If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= +// IL <= IU <= N, m, // IU: If RANGE='I', the index of the largest eigenvalue to be +// returned. 1 <= IL <= IU <= N, 0.0, // ABSTOL m, // M: The total number of +// eigenvalues found. 0 <= M <= N. if RANGE = 'I', M = IU-IL+1. W, // W store +// eigenvalues V, // store eigenvector col, // LDZ: The leading dimension of the +// array Z. work, lwork, rwork, iwork, ifail, info, row); + +// delete[] work; +// delete[] rwork; +// delete[] iwork; +// delete[] ifail; + +// }; + +} // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/src/math_kernel.cpp b/source/module_hsolver/src/math_kernel.cpp index c9fec3e8470..d813af348b4 100644 --- a/source/module_hsolver/src/math_kernel.cpp +++ b/source/module_hsolver/src/math_kernel.cpp @@ -208,6 +208,52 @@ void gemm_op::operator()(const psi::DEVICE_CPU* d, zgemm_(&transa, &transb, &m, &n ,&k, alpha, a ,&lda, b, &ldb, beta, c, &ldc); } +template struct matrixTranspose_op +{ + void operator()(const psi::DEVICE_CPU* d, + const int& row, + const int& col, + const std::complex* input_matrix, + std::complex* output_matrix) + { + std::complex* temp = nullptr; + psi::memory::resize_memory_op, psi::DEVICE_CPU>()(d, temp, row * col); + for (int i = 0; i < row; i++) + { + for (int j = 0; j < col; j++) + { + temp[j * row + i] = input_matrix[i * col + j]; + } + } + for (int i = 0; i < row * col; i++) + { + output_matrix[i] = temp[i]; + } + psi::memory::delete_memory_op, psi::DEVICE_CPU>()(d, temp); + } +}; + +template struct matrixSetToAnother +{ + void operator()(const psi::DEVICE_CPU* d, + const int& n, + const std::complex* A, + const int& LDA, + std::complex* B, + const int& LDB) + { + for (int i = 0; i < n; i++) + { + for (int j = 0; j < LDA; j++) + { + B[i * LDB + j] = A[i * LDA + j]; + } + } + } +}; + + + // Explicitly instantiate functors for the types of functor registered. template struct zdot_real_op; template struct vector_div_constant_op; @@ -215,5 +261,8 @@ template struct vector_mul_vector_op; template struct vector_div_vector_op; template struct constantvector_addORsub_constantVector_op; +template struct matrixTranspose_op; +template struct matrixSetToAnother; + } // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/test/CMakeLists.txt b/source/module_hsolver/test/CMakeLists.txt index 247a78f5519..b7b71d4ab0a 100644 --- a/source/module_hsolver/test/CMakeLists.txt +++ b/source/module_hsolver/test/CMakeLists.txt @@ -9,6 +9,7 @@ AddTest( ../../module_hamilt/operator.cpp ../../module_hamilt/ks_pw/operator_pw.cpp ../src/math_kernel.cpp + ../src/dngvd_op.cpp ) AddTest( TARGET HSolver_dav @@ -19,6 +20,7 @@ AddTest( ../../module_hamilt/operator.cpp ../../module_hamilt/ks_pw/operator_pw.cpp ../src/math_kernel.cpp + ../src/dngvd_op.cpp ) if(ENABLE_LCAO) @@ -36,16 +38,32 @@ AddTest( LIBS ${math_libs} base cublas SOURCES math_kernel_test.cpp ../../src_parallel/parallel_global.cpp ../../src_parallel/parallel_common.cpp ../../src_parallel/parallel_reduce.cpp + ../../module_psi/src/memory_psi.cpp ../../module_hsolver/src/math_kernel.cpp ../../module_psi/src/cuda/memory.cu ../src/cuda/math_kernel.cu ) +AddTest( + TARGET Dngvd_UTs + LIBS ${math_libs} base cublas cusolver + SOURCES math_dngvd_test.cpp + ../../src_parallel/parallel_global.cpp + ../../src_parallel/parallel_common.cpp ../../src_parallel/parallel_reduce.cpp + ../../module_psi/src/memory_psi.cpp + ../../module_psi/src/cuda/memory.cu + ../../module_hsolver/src/math_kernel.cpp + ../src/cuda/math_kernel.cu + ../../module_hsolver/src/dngvd_op.cpp + ../src/cuda/dngvd_op.cu + +) else() AddTest( TARGET Hsolver_UTs LIBS ${math_libs} base SOURCES math_kernel_test.cpp ../../src_parallel/parallel_global.cpp ../../src_parallel/parallel_common.cpp ../../src_parallel/parallel_reduce.cpp + ../../module_psi/src/memory_psi.cpp ../src/math_kernel.cpp ) endif() diff --git a/source/module_hsolver/test/math_dngvd_test.cpp b/source/module_hsolver/test/math_dngvd_test.cpp new file mode 100644 index 00000000000..805987f23a9 --- /dev/null +++ b/source/module_hsolver/test/math_dngvd_test.cpp @@ -0,0 +1,648 @@ +#include "module_base/complexmatrix.h" +#include "module_base/lapack_connector.h" +#include "module_hsolver/include/dngvd_op.h" +#include "module_hsolver/include/math_kernel.h" +#include "module_psi/include/memory.h" + +#include +#include +#include +#include +#include +#include + +class TestModuleHsolverMathDngvd : public ::testing::Test +{ + protected: + using resize_memory_op_Z = psi::memory::resize_memory_op, psi::DEVICE_GPU>; + using delete_memory_op_Z = psi::memory::delete_memory_op, psi::DEVICE_GPU>; + using resize_memory_op_D = psi::memory::resize_memory_op; + using delete_memory_op_D = psi::memory::delete_memory_op; + // from CPU to GPU + using synchronize_memory_op_C2G_Z + = psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_CPU>; + using synchronize_memory_op_C2G_D = psi::memory::synchronize_memory_op; + using synchronize_memory_op_G2C_Z + = psi::memory::synchronize_memory_op, psi::DEVICE_CPU, psi::DEVICE_GPU>; + using synchronize_memory_op_G2C_D = psi::memory::synchronize_memory_op; + + const psi::DEVICE_CPU* cpu_ctx = {}; + const psi::DEVICE_GPU* gpu_ctx = {}; + + // prepare A & B in CPU + std::vector> matrix_A = { + {-0.351417, -1.73472}, + {-8.32667, 2.3744}, + {4.16334, 3.64292}, + {5.20417, -3.85976}, + {-8.32667, -2.3744}, + {0.551651, -2.60209}, + {2.08167, 1.9082}, + {-6.93889, 1.04083}, + {4.16334, -3.64292}, + {2.08167, -1.9082}, + {0.551651, -2.25514}, + {-1.31839, 5.20417}, + {5.20417, 3.85976}, + {-6.93889, -1.04083}, + {-1.31839, -5.20417}, + {0.551651, -3.64292} + + // {-4.280e-01,3.084e-17}, {-1.288e-16,9.021e-17}, {5.204e-18,-3.990e-17}, {-5.204e-17,-2.776e-17}, + // {-1.288e-16,-9.021e-17}, {4.574e-01,-8.687e-17}, {-6.939e-18,1.908e-17}, {8.327e-17,-1.041e-17}, + // {5.204e-18,3.990e-17}, {-6.939e-18,-1.908e-17}, {4.574e-01,-5.783e-17},{-6.939e-18,4.163e-17}, + // {-5.204e-17,2.776e-17}, {8.327e-17,1.041e-17}, {-6.939e-18,-4.163e-17}, {4.574e-01,-3.451e-17} + + // {-4.280e-01,3.084e-17}, {-1.288e-16,9.021e-17}, {5.204e-18,-3.990e-17}, {-1.145e-16,2.255e-17}, + // {3.946e-17,-2.949e-17}, {4.574e-01,-8.687e-17}, {2.082e-17,1.908e-17}, {4.163e-17,-4.163e-17}, + // {3.296e-17,-9.454e-17}, {-6.939e-18,-1.908e-17}, {4.574e-01,-5.783e-17},{-1.388e-17,6.939e-18}, + // {-5.204e-17,2.776e-17}, {8.327e-17,1.041e-17}, {-6.939e-18,-4.163e-17}, {4.574e-01,-3.451e-17} + }; + std::vector> matrix_B = { + {1, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {1, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {1, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {0, 0}, + {1, 0} + // {1.000e+00,0.000e+00}, {-7.069e-17,-3.123e-17}, {-5.204e-17,0.000e+00}, {5.551e-17,-2.082e-17}, + // {-7.069e-17,3.123e-17}, {1.000e+00,0.000e+00}, {1.110e-16,-7.286e-17}, {8.327e-17,-1.110e-16}, + // {-5.204e-17,0.000e+00}, {1.110e-16,7.286e-17}, {1.000e+00,0.000e+00}, {-9.714e-17,2.776e-17}, + // {5.551e-17,2.082e-17}, {8.327e-17,1.110e-16}, {-9.714e-17,-2.776e-17}, {1.000e+00,0.000e+00} + }; + const int matrix_size = 16; + + // prepare W & V in CPU in dngv_op + std::vector W_dngv_op = {0.0, 0.0, 0.0, 0.0}; + std::vector> matrix_V_dngv_op = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; + + // prepare W & V in CPU in dngvx_op + std::vector W_DNGVX = {0.0, 0.0}; + std::vector> matrix_V_DNGVX = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; +}; + +TEST_F(TestModuleHsolverMathDngvd, dngv_cpu) +{ + // (1) + ModuleBase::ComplexMatrix B(4, 4); + ModuleBase::ComplexMatrix V(4, 4); + // ================================== + int col = 4; + int row = 4; + int lwork = 0; + int nb = LapackConnector::ilaenv(1, "ZHETRD", "L", col, -1, -1, -1); + if (nb < 1) + { + nb = std::max(1, col); + } + if (nb == 1 || nb >= col) + { + lwork = 2 * col; // mohan modify 2009-08-02 + } + else + { + lwork = (nb + 1) * col; + } + std::complex* work = new std::complex[lwork]; + ModuleBase::GlobalFunc::ZEROS(work, lwork); + int info = 0; + int rwork_dim; + rwork_dim = 3 * col - 2; + double* rwork = new double[rwork_dim]; + ModuleBase::GlobalFunc::ZEROS(rwork, rwork_dim); + // V.c = matrix_A.data(); + // B.c = matrix_B.data(); + for (int i = 0; i < matrix_size; i++) + { + V.c[i] = matrix_A[i]; + } + for (int i = 0; i < matrix_size; i++) + { + B.c[i] = matrix_B[i]; + } + LapackConnector::zhegv(1, 'V', 'U', col, V, col, B, col, W_dngv_op.data(), work, lwork, rwork, info); + // ================================== + + // (2) + std::vector W_result = {0.0, 0.0, 0.0, 0.0}; + std::vector> V_result = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; + hsolver::dngv_op()(cpu_ctx, + 4, + 4, + matrix_A.data(), + matrix_B.data(), + W_result.data(), + V_result.data()); + + // output + // std::cout << W_dngv_op[0] << "\t" << W_dngv_op[1] << W_dngv_op[2] << W_dngv_op[3] << std::endl; + // std::cout << W_result[0] << "\t" << W_result[1] << W_result[2] << W_result[3] << std::endl; + // std::cout << "CPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V_result[i * 4 + j] << ", "; + // } + // std::cout << std::endl; + // } + // std::cout << "GPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V(i, j) << ", "; + // } + // std::cout << std::endl; + // } + + // test + for (int i = 0; i < W_dngv_op.size(); i++) + { + EXPECT_LT(fabs(W_dngv_op[i] - W_result[i]), 1e-8); + } + for (int i = 0; i < V_result.size(); i++) + { + EXPECT_LT(fabs(V.c[i].imag() - V_result[i].imag()), 1e-8); + EXPECT_LT(fabs(V.c[i].real() - V_result[i].real()), 1e-8); + } + delete[] rwork; + delete[] work; +} + +TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) +{ + // (1) + ModuleBase::ComplexMatrix A(4, 4); + ModuleBase::ComplexMatrix B(4, 4); + ModuleBase::ComplexMatrix V(4, 4); + // ================================== + int lwork; + int info = 0; + std::string name1 = "ZHETRD"; + std::string name2 = "L"; + int col = 4; + int nb = LapackConnector::ilaenv(1, name1.c_str(), name2.c_str(), col, -1, -1, -1); + if (nb < 1) + { + nb = std::max(1, col); + } + if (nb == 1 || nb >= col) + { + lwork = 2 * col; // qianrui fix a bug 2021-7-25 : lwork should be at least max(1,2*n) + } + else + { + lwork = (nb + 1) * col; + } + std::complex* work = new std::complex[2 * lwork]; + assert(work != 0); + double* rwork = new double[7 * col]; + assert(rwork != 0); + int* iwork = new int[5 * col]; + assert(iwork != 0); + int* ifail = new int[col]; + assert(ifail != 0); + ModuleBase::GlobalFunc::ZEROS(work, lwork); // qianrui change it, only first lwork numbers are used in zhegvx + ModuleBase::GlobalFunc::ZEROS(rwork, 7 * col); + ModuleBase::GlobalFunc::ZEROS(iwork, 5 * col); + ModuleBase::GlobalFunc::ZEROS(ifail, col); + for (int i = 0; i < matrix_size; i++) + { + A.c[i] = matrix_A[i]; + } + for (int i = 0; i < matrix_size; i++) + { + B.c[i] = matrix_B[i]; + } + for (int i = 0; i < matrix_size; i++) + { + V.c[i] = 0.0; + } + int m = 2; + LapackConnector::zhegvx( + 1, // ITYPE = 1: A*x = (lambda)*B*x + 'V', // JOBZ = 'V': Compute eigenvalues and eigenvectors. + 'I', // RANGE = 'I': the IL-th through IU-th eigenvalues will be found. + 'L', // UPLO = 'L': Lower triangles of A and B are stored. + col, // N = base + A, // A is COMPLEX*16 array dimension (LDA, N) + col, // LDA = base + B, // B is COMPLEX*16 array, dimension (LDB, N) + col, // LDB = base + 0.0, // Not referenced if RANGE = 'A' or 'I'. + 0.0, // Not referenced if RANGE = 'A' or 'I'. + 1, // IL: If RANGE='I', the index of the smallest eigenvalue to be returned. 1 <= IL <= IU <= N, + m, // IU: If RANGE='I', the index of the largest eigenvalue to be returned. 1 <= IL <= IU <= N, + 0.0, // ABSTOL + m, // M: The total number of eigenvalues found. 0 <= M <= N. if RANGE = 'I', M = IU-IL+1. + W_DNGVX.data(), // W store eigenvalues + V, // store eigenvector + col, // LDZ: The leading dimension of the array Z. + work, + lwork, + rwork, + iwork, + ifail, + info); + + // ================================== + + // (2) + std::vector W_result = {0.0, 0.0}; + std::vector> V_result = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; + hsolver::dngvx_op()(cpu_ctx, + 4, + 4, + matrix_A.data(), + matrix_B.data(), + 2, + W_result.data(), + V_result.data()); + + // output + // std::cout << W_result[0] << "\t" << W_result[1] << std::endl; + // std::cout << W_DNGVX[0] << "\t" << W_DNGVX[1]<< std::endl; + // std::cout << "CPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V_result[i * 4 + j] << ", "; + // } + // std::cout << std::endl; + // } + // std::cout << "GPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V(i, j) << ", "; + // } + // std::cout << std::endl; + // } + + // test + for (int i = 0; i < W_result.size(); i++) + { + EXPECT_LT(fabs(W_DNGVX[i] - W_result[i]), 1e-8); + } + for (int i = 0; i < V_result.size(); i++) + { + EXPECT_LT(fabs(V.c[i].imag() - V_result[i].imag()), 1e-8); + EXPECT_LT(fabs(V.c[i].real() - V_result[i].real()), 1e-8); + } + delete[] work; + delete[] rwork; + delete[] iwork; + delete[] ifail; +} + +#if __UT_USE_CUDA || __UT_USE_ROCM + +// computes all the eigenvalues and eigenvectors +TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) +{ + // prepare A & B in GPU + std::complex* device_matrix_A = nullptr; + std::complex* device_matrix_B = nullptr; + resize_memory_op_Z()(gpu_ctx, device_matrix_A, matrix_size); + resize_memory_op_Z()(gpu_ctx, device_matrix_B, matrix_size); + synchronize_memory_op_C2G_Z()(gpu_ctx, cpu_ctx, device_matrix_A, matrix_A.data(), matrix_size); + synchronize_memory_op_C2G_Z()(gpu_ctx, cpu_ctx, device_matrix_B, matrix_B.data(), matrix_size); + + // prepare W & V in GPU in dngv_op + double* device_W_dngv_op = nullptr; + resize_memory_op_D()(gpu_ctx, device_W_dngv_op, W_dngv_op.size()); + psi::memory::set_memory_op()(gpu_ctx, device_W_dngv_op, 0, W_dngv_op.size()); + std::complex* device_matrix_V_dngv_op = nullptr; + resize_memory_op_Z()(gpu_ctx, device_matrix_V_dngv_op, matrix_V_dngv_op.size()); + psi::memory::set_memory_op, psi::DEVICE_GPU>()(gpu_ctx, + device_matrix_V_dngv_op, + 0, + matrix_V_dngv_op.size()); + + // run in GPU + hsolver::dngv_op()(gpu_ctx, + 4, + 4, + device_matrix_A, + device_matrix_B, + device_W_dngv_op, + device_matrix_V_dngv_op); + // copy W data from GPU to CPU + std::vector W_result = {0.0, 0.0, 0.0, 0.0}; + synchronize_memory_op_G2C_D()(cpu_ctx, gpu_ctx, W_result.data(), device_W_dngv_op, W_result.size()); + // copy V data from GPU to CPU + std::vector> V_result = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; + synchronize_memory_op_G2C_Z()(cpu_ctx, gpu_ctx, V_result.data(), device_matrix_V_dngv_op, V_result.size()); + + // run in CPU + hsolver::dngv_op()(cpu_ctx, + 4, + 4, + matrix_A.data(), + matrix_B.data(), + W_dngv_op.data(), + matrix_V_dngv_op.data()); + + // // CPU + // std::cout << W_dngv_op[0] << "\t" << W_dngv_op[1] << "\t" << W_dngv_op[2] << "\t" << W_dngv_op[3] << + // std::endl; + // // GPU + // std::cout << W_result[0] << "\t" << W_result[1] << "\t" << W_result[2] << "\t" << W_result[3] << std::endl; + // std::cout << "CPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << matrix_V_dngv_op[i * 4 + j] << ", "; + + // } + // std::cout << std::endl; + // } + // std::cout << "GPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V_result[i * 4 + j] << ", "; + // } + // std::cout << std::endl; + // } + + // we need to compare + // 1. W with W_result + // 2. matrix_V with V_result + for (int i = 0; i < W_dngv_op.size(); i++) + { + EXPECT_LT(fabs(W_dngv_op[i] - W_result[i]), 1e-8); + } + for (int i = 0; i < matrix_V_dngv_op.size(); i++) + { + EXPECT_LT(fabs(matrix_V_dngv_op[i].imag()) - fabs(V_result[i].imag()), 1e-8); + EXPECT_LT(fabs(matrix_V_dngv_op[i].real()) - fabs(V_result[i].real()), 1e-8); + } + + // delete values in GPU + delete_memory_op_Z()(gpu_ctx, device_matrix_A); + delete_memory_op_Z()(gpu_ctx, device_matrix_B); + delete_memory_op_Z()(gpu_ctx, device_matrix_V_dngv_op); + delete_memory_op_D()(gpu_ctx, device_W_dngv_op); +} + +// computes the first m eigenvalues ​​and their corresponding eigenvectors +TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) +{ + // prepare A & B in GPU + std::complex* device_matrix_A = nullptr; + std::complex* device_matrix_B = nullptr; + resize_memory_op_Z()(gpu_ctx, device_matrix_A, matrix_size); + resize_memory_op_Z()(gpu_ctx, device_matrix_B, matrix_size); + synchronize_memory_op_C2G_Z()(gpu_ctx, cpu_ctx, device_matrix_A, matrix_A.data(), matrix_size); + synchronize_memory_op_C2G_Z()(gpu_ctx, cpu_ctx, device_matrix_B, matrix_B.data(), matrix_size); + // prepare W & V in GPU in dngvx_op + double* device_W_DNGVX = nullptr; + resize_memory_op_D()(gpu_ctx, device_W_DNGVX, W_DNGVX.size()); + synchronize_memory_op_C2G_D()(gpu_ctx, cpu_ctx, device_W_DNGVX, W_DNGVX.data(), W_DNGVX.size()); + std::complex* device_matrix_V_DNGVX = nullptr; + resize_memory_op_Z()(gpu_ctx, device_matrix_V_DNGVX, matrix_V_DNGVX.size()); + synchronize_memory_op_C2G_Z()(gpu_ctx, + cpu_ctx, + device_matrix_V_DNGVX, + matrix_V_DNGVX.data(), + matrix_V_DNGVX.size()); + + // run in GPU + hsolver::dngvx_op()(gpu_ctx, + 4, + 4, + device_matrix_A, + device_matrix_B, + 2, + device_W_DNGVX, + device_matrix_V_DNGVX); + // copy W data from GPU to CPU + std::vector W_result = {0.0, 0.0}; + synchronize_memory_op_G2C_D()(cpu_ctx, gpu_ctx, W_result.data(), device_W_DNGVX, W_result.size()); + // copy V data from GPU to CPU + std::vector> V_result = {{0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}}; + synchronize_memory_op_G2C_Z()(cpu_ctx, gpu_ctx, V_result.data(), device_matrix_V_DNGVX, V_result.size()); + + // run in CPU + hsolver::dngvx_op()(cpu_ctx, + 4, + 4, + matrix_A.data(), + matrix_B.data(), + 2, + W_DNGVX.data(), + matrix_V_DNGVX.data()); + + // std::cout << "GPU::" << std::endl; + // std::cout << W_result[0] << "\t" << W_result[1] << std::endl; + // std::cout << "CPU::" << std::endl; + // std::cout << W_DNGVX[0] << "\t" << W_DNGVX[1] << std::endl; + // std::cout << "GPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << V_result[i * 4 + j] << ", "; + // } + // std::cout << std::endl; + // } + // std::cout << "CPU vector" << std::endl; + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 4; j++) + // { + // std::cout << matrix_V_DNGVX[i * 4 + j] << ", "; + // } + // std::cout << std::endl; + // } + + // we need to compare + // 1. W with W_result + // 2. matrix_V with V_result + for (int i = 0; i < W_DNGVX.size(); i++) + { + EXPECT_LT(fabs(W_DNGVX[i] - W_result[i]), 1e-8); + } + for (int i = 0; i < matrix_V_dngv_op.size(); i++) + { + EXPECT_LT(fabs(matrix_V_dngv_op[i].imag()) - fabs(V_result[i].imag()), 1e-8); + EXPECT_LT(fabs(matrix_V_dngv_op[i].real()) - fabs(V_result[i].real()), 1e-8); + } + + // delete values in GPU + delete_memory_op_Z()(gpu_ctx, device_matrix_A); + delete_memory_op_Z()(gpu_ctx, device_matrix_B); + delete_memory_op_Z()(gpu_ctx, device_matrix_V_DNGVX); + delete_memory_op_D()(gpu_ctx, device_W_DNGVX); +} + +TEST_F(TestModuleHsolverMathDngvd, transpose_gpu) +{ + // prepare transpose in GPU + std::vector> transpose = { + {-0.351417, -1.73472}, + {-8.32667, 2.3744}, + {4.16334, 3.64292}, + {-0.351417, -1.73472}, + {-8.32667, 2.3744}, + {4.16334, 3.64292}, + // {-0.351417,-1.73472}, {-8.32667,2.3744}, {4.16334,3.64292} + }; + std::complex* device_transpose = nullptr; + resize_memory_op_Z()(gpu_ctx, device_transpose, matrix_size); + synchronize_memory_op_C2G_Z()(gpu_ctx, cpu_ctx, device_transpose, transpose.data(), transpose.size()); + + // run + hsolver::createBLAShandle(); + hsolver::matrixTranspose_op()(gpu_ctx, 2, 3, device_transpose, device_transpose); + hsolver::destoryBLAShandle(); + + // copy transpose data from GPU to CPU + std::vector> transpose_result = { + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.0, 0.0}, + // {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0} + }; + synchronize_memory_op_G2C_Z()(cpu_ctx, gpu_ctx, transpose_result.data(), device_transpose, transpose.size()); + + // std::vector > test_result = { + // {-0.351417,-1.73472}, {-0.351417,-1.73472}, {-0.351417,-1.73472}, + // {-8.32667,2.3744}, {-8.32667,2.3744}, {-8.32667,2.3744}, + // {4.16334,3.64292}, {4.16334,3.64292}, {4.16334,3.64292}, + // }; + std::vector> test_result = { + {-0.351417, -1.73472}, + {-0.351417, -1.73472}, + {-8.32667, 2.3744}, + {-8.32667, 2.3744}, + {4.16334, 3.64292}, + {4.16334, 3.64292}, + }; + + // for (int i = 0; i < 3; i++) + // { + // for (int j = 0; j < 2; j++) + // { + // std::cout << transpose_result[i * 2 + j]; + // } + // std::cout << std::endl; + // } + + for (int i = 0; i < transpose_result.size(); i++) + { + EXPECT_LT(fabs(test_result[i].imag()) - fabs(transpose_result[i].imag()), 1e-8); + EXPECT_LT(fabs(test_result[i].real()) - fabs(transpose_result[i].real()), 1e-8); + } +} + +#endif // __UT_USE_CUDA || __UT_USE_ROCM diff --git a/source/module_hsolver/test/math_kernel_test.cpp b/source/module_hsolver/test/math_kernel_test.cpp index f24e285f7a7..1b05d5d0e43 100644 --- a/source/module_hsolver/test/math_kernel_test.cpp +++ b/source/module_hsolver/test/math_kernel_test.cpp @@ -1,44 +1,54 @@ -#include -#include -#include -#include "module_psi/include/memory.h" -#include "module_hsolver/include/math_kernel.h" #include "module_base/blas_connector.h" #include "module_base/constants.h" +#include "module_hsolver/include/math_kernel.h" +#include "module_psi/include/memory.h" +#include +#include +#include class TestModuleHsolverMathKernel : public ::testing::Test { protected: // xx = tf.random.uniform([100], minval=-4, maxval=4, dtype = tf.float64) - std::vector > psi_L = { - {-0.65412617, -0.74208893}, {-2.21731157, 0.42540039}, - {3.36373004, -2.51647562}, {-2.985111 , -0.53251562}, - {0.37908265, 0.81605825}, { 1.66281318, 2.71761869}, - {2.2010268 , 0.65498149}, { 1.51153638, 0.71501482}, - {0.53546578, 1.4564317 }, {-2.36701143, 1.23009056}, - {3.41302551, -2.3175205 }, {-0.27628221, -1.35701656} - }; - - std::vector > psi_R = { - {-1.67837557e-01, -1.70017454e-01}, {-2.92128115e-02, 2.82765887e-01}, - {-8.71641062e-02, -1.15934278e-01}, { 3.36269232e+00, -1.44692661e-02}, - {-3.81342874e-03, -1.58276988e-01}, { 2.33504238e-01, -1.93195840e-03}, - { 2.45520665e-01, 6.46854620e-01}, { 1.58255340e+00, 2.70915699e+00}, - {-1.66142311e-01, 6.27839507e-02}, { 2.17077193e+00, 4.87104731e-01}, - { 1.41257916e+00, 5.45282609e-01}, {-1.29333636e-01, -5.04228492e-03} - }; + std::vector> psi_L = {{-0.65412617, -0.74208893}, + {-2.21731157, 0.42540039}, + {3.36373004, -2.51647562}, + {-2.985111, -0.53251562}, + {0.37908265, 0.81605825}, + {1.66281318, 2.71761869}, + {2.2010268, 0.65498149}, + {1.51153638, 0.71501482}, + {0.53546578, 1.4564317}, + {-2.36701143, 1.23009056}, + {3.41302551, -2.3175205}, + {-0.27628221, -1.35701656}}; + + std::vector> psi_R = {{-1.67837557e-01, -1.70017454e-01}, + {-2.92128115e-02, 2.82765887e-01}, + {-8.71641062e-02, -1.15934278e-01}, + {3.36269232e+00, -1.44692661e-02}, + {-3.81342874e-03, -1.58276988e-01}, + {2.33504238e-01, -1.93195840e-03}, + {2.45520665e-01, 6.46854620e-01}, + {1.58255340e+00, 2.70915699e+00}, + {-1.66142311e-01, 6.27839507e-02}, + {2.17077193e+00, 4.87104731e-01}, + {1.41257916e+00, 5.45282609e-01}, + {-1.29333636e-01, -5.04228492e-03}}; const int dim = psi_L.size(); const double expected_result = -5.0016151713691288; - const psi::DEVICE_CPU * cpu_ctx = {}; - const psi::DEVICE_GPU * gpu_ctx = {}; + const psi::DEVICE_CPU* cpu_ctx = {}; + const psi::DEVICE_GPU* gpu_ctx = {}; - void SetUp() override { + void SetUp() override + { } - void TearDown() override { + void TearDown() override + { } using zdot_real_cpu_op = hsolver::zdot_real_op; @@ -47,22 +57,25 @@ class TestModuleHsolverMathKernel : public ::testing::Test using resize_memory_op = psi::memory::resize_memory_op, psi::DEVICE_GPU>; using delete_memory_op = psi::memory::delete_memory_op, psi::DEVICE_GPU>; // from CPU to GPU - using synchronize_memory_op = psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_CPU>; - + using synchronize_memory_op + = psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_CPU>; + // form GPU to CPU by haozhihan - using synchronize_memory_op_gpu = psi::memory::synchronize_memory_op, psi::DEVICE_CPU, psi::DEVICE_GPU>; + using synchronize_memory_op_gpu + = psi::memory::synchronize_memory_op, psi::DEVICE_CPU, psi::DEVICE_GPU>; // about double operator by haozhihan using resize_memory_op_double = psi::memory::resize_memory_op; using delete_memory_op_double = psi::memory::delete_memory_op; using synchronize_memory_op_double = psi::memory::synchronize_memory_op; - -// haozhihan add + + // haozhihan add // cpu operator using vector_div_constant_op_cpu = hsolver::vector_div_constant_op; using vector_mul_vector_op_cpu = hsolver::vector_mul_vector_op; using vector_div_vector_op_cpu = hsolver::vector_div_vector_op; - using constantvector_addORsub_constantVector_op_cpu = hsolver::constantvector_addORsub_constantVector_op; + using constantvector_addORsub_constantVector_op_cpu + = hsolver::constantvector_addORsub_constantVector_op; using axpy_op_cpu = hsolver::axpy_op; using scal_op_cpu = hsolver::scal_op; using gemv_op_cpu = hsolver::gemv_op; @@ -70,478 +83,630 @@ class TestModuleHsolverMathKernel : public ::testing::Test using vector_div_constant_op_gpu = hsolver::vector_div_constant_op; using vector_mul_vector_op_gpu = hsolver::vector_mul_vector_op; using vector_div_vector_op_gpu = hsolver::vector_div_vector_op; - using constantvector_addORsub_constantVector_op_gpu = hsolver::constantvector_addORsub_constantVector_op; + using constantvector_addORsub_constantVector_op_gpu + = hsolver::constantvector_addORsub_constantVector_op; using axpy_op_gpu = hsolver::axpy_op; using scal_op_gpu = hsolver::scal_op; using gemv_op_gpu = hsolver::gemv_op; - - // haozhihan add - std::vector > L = { - {-0.65412617, -0.74208893}, {-2.21731157, 0.42540039}, - {3.36373004, -2.51647562}, {-2.985111 , -0.53251562}, - {0.37908265, 0.81605825}, { 1.66281318, 2.71761869}, - {2.2010268 , 0.65498149}, { 1.51153638, 0.71501482}, - {0.53546578, 1.4564317 }, {-2.36701143, 1.23009056}, - {3.41302551, -2.3175205 }, {-0.27628221, -1.35701656} - }; - - std::vector > R = { - {-1.67837557e-01, -1.70017454e-01}, {-2.92128115e-02, 2.82765887e-01}, - {-8.71641062e-02, -1.15934278e-01}, { 3.36269232e+00, -1.44692661e-02}, - {-3.81342874e-03, -1.58276988e-01}, { 2.33504238e-01, -1.93195840e-03}, - { 2.45520665e-01, 6.46854620e-01}, { 1.58255340e+00, 2.70915699e+00}, - {-1.66142311e-01, 6.27839507e-02}, { 2.17077193e+00, 4.87104731e-01}, - { 1.41257916e+00, 5.45282609e-01}, {-1.29333636e-01, -5.04228492e-03} - }; + std::vector> L = {{-0.65412617, -0.74208893}, + {-2.21731157, 0.42540039}, + {3.36373004, -2.51647562}, + {-2.985111, -0.53251562}, + {0.37908265, 0.81605825}, + {1.66281318, 2.71761869}, + {2.2010268, 0.65498149}, + {1.51153638, 0.71501482}, + {0.53546578, 1.4564317}, + {-2.36701143, 1.23009056}, + {3.41302551, -2.3175205}, + {-0.27628221, -1.35701656}}; + + std::vector> R = {{-1.67837557e-01, -1.70017454e-01}, + {-2.92128115e-02, 2.82765887e-01}, + {-8.71641062e-02, -1.15934278e-01}, + {3.36269232e+00, -1.44692661e-02}, + {-3.81342874e-03, -1.58276988e-01}, + {2.33504238e-01, -1.93195840e-03}, + {2.45520665e-01, 6.46854620e-01}, + {1.58255340e+00, 2.70915699e+00}, + {-1.66142311e-01, 6.27839507e-02}, + {2.17077193e+00, 4.87104731e-01}, + {1.41257916e+00, 5.45282609e-01}, + {-1.29333636e-01, -5.04228492e-03}}; // (1) for test vector_div_constant_op - const std::vector > input = L; + const std::vector> input = L; const double constant = 5.5; - const std::vector > output_vector_div_constant_op = { - {-0.11893203,-0.13492526}, {-0.40314756, 0.07734553}, - {0.61158728, -0.45754102}, {-0.54274745,-0.09682102}, - {0.06892412, 0.14837423}, {0.30232967, 0.49411249}, - {0.40018669, 0.11908754}, {0.27482480, 0.13000269}, - {0.09735741, 0.26480576}, {-0.43036571, 0.22365283}, - {0.62055009, -0.42136736}, {-0.05023313,-0.24673028} - }; + const std::vector> output_vector_div_constant_op = {{-0.11893203, -0.13492526}, + {-0.40314756, 0.07734553}, + {0.61158728, -0.45754102}, + {-0.54274745, -0.09682102}, + {0.06892412, 0.14837423}, + {0.30232967, 0.49411249}, + {0.40018669, 0.11908754}, + {0.27482480, 0.13000269}, + {0.09735741, 0.26480576}, + {-0.43036571, 0.22365283}, + {0.62055009, -0.42136736}, + {-0.05023313, -0.24673028}}; // (2) for test vector_mul_vector_op & vector_div_vector_op const std::vector input_double = { - -0.65412617, -0.74208893, -2.21731157, 0.42540039, - 3.36373004, -2.51647562, -2.985111 , -0.53251562, - 0.37908265, 0.81605825, 1.66281318, 2.71761869, - }; - const std::vector >output_vector_mul_vector_op = { - {0.42788105, 0.48541979}, {1.64544237, -0.31568492}, - {-7.45843754, 5.57981051}, {-1.26986738,-0.22653235}, - {1.27513170, 2.74499965}, {-4.18442883,-6.83882118}, - {-6.57030931,-1.95519245}, {-0.80491673,-0.38075656}, - {0.20298579, 0.55210799}, {-1.93161921, 1.00382555}, - {5.67522380, -3.85360363}, {-0.75082970,-3.68785357} - }; - - const std::vector >output_vector_div_vector_op = { - {1.00000000, 1.13447369}, {2.98793242, -0.57324718}, - {-1.51703084, 1.13492197}, {-7.01717974,-1.25179862}, - {0.11269711, 0.24260516}, {-0.66077063,-1.07993047}, - {-0.73733499,-0.21941613}, {-2.83848271,-1.34271145}, - {1.41253043, 3.84198987}, {-2.90054225, 1.50735632}, - {2.05256102, -1.39373474}, {-0.10166335,-0.49934031} + -0.65412617, + -0.74208893, + -2.21731157, + 0.42540039, + 3.36373004, + -2.51647562, + -2.985111, + -0.53251562, + 0.37908265, + 0.81605825, + 1.66281318, + 2.71761869, }; + const std::vector> output_vector_mul_vector_op = {{0.42788105, 0.48541979}, + {1.64544237, -0.31568492}, + {-7.45843754, 5.57981051}, + {-1.26986738, -0.22653235}, + {1.27513170, 2.74499965}, + {-4.18442883, -6.83882118}, + {-6.57030931, -1.95519245}, + {-0.80491673, -0.38075656}, + {0.20298579, 0.55210799}, + {-1.93161921, 1.00382555}, + {5.67522380, -3.85360363}, + {-0.75082970, -3.68785357}}; + + const std::vector> output_vector_div_vector_op = {{1.00000000, 1.13447369}, + {2.98793242, -0.57324718}, + {-1.51703084, 1.13492197}, + {-7.01717974, -1.25179862}, + {0.11269711, 0.24260516}, + {-0.66077063, -1.07993047}, + {-0.73733499, -0.21941613}, + {-2.83848271, -1.34271145}, + {1.41253043, 3.84198987}, + {-2.90054225, 1.50735632}, + {2.05256102, -1.39373474}, + {-0.10166335, -0.49934031}}; // (3) for test constantvector_addORsub_constantVector_op const double constant1 = 6.6; const double constant2 = 4.4; - const std::vector > input1 = L; - const std::vector > input2 = R; - const std::vector >output_constantvector_addORsub_constantVector_op = { - {-5.05571797, -5.64586374}, {-14.76279273,4.05181248}, - {21.81709620,-17.11884992}, {-4.90588639,-3.57826786}, - {2.48516640, 4.68956570}, {12.00198564,17.92778274}, - {15.60706781, 7.16903816}, {16.93937507,16.63938857}, - {2.80304798, 9.88869860}, {-6.07087895,10.26185851}, - {28.74131667,-12.89639182}, {-2.39253058,-8.97849535} - }; - - // (4) for test axpy_op (compute Y = alpha * X + Y ) + const std::vector> input1 = L; + const std::vector> input2 = R; + const std::vector> output_constantvector_addORsub_constantVector_op + = {{-5.05571797, -5.64586374}, + {-14.76279273, 4.05181248}, + {21.81709620, -17.11884992}, + {-4.90588639, -3.57826786}, + {2.48516640, 4.68956570}, + {12.00198564, 17.92778274}, + {15.60706781, 7.16903816}, + {16.93937507, 16.63938857}, + {2.80304798, 9.88869860}, + {-6.07087895, 10.26185851}, + {28.74131667, -12.89639182}, + {-2.39253058, -8.97849535}}; + + // (4) for test axpy_op (compute Y = alpha * X + Y ) const std::complex alpha_axpy{-1.5, -2.5}; - const std::vector > X_axpy = L; - std::vector > Y_axpy = R; - - const std::vector >output_axpy_op = { - {-1.04187063, 2.57843137}, {4.36025552, 5.18794423}, - {-11.42394822,-4.75054595}, {6.50906977, 8.24708166}, - {1.46770822, -2.33007099}, {4.53333119,-8.23539294}, - {-1.41856581, -5.83818462}, {1.10278588,-2.14220619}, - {2.67173827, -3.46052805}, {8.79651547, 4.55949747}, - {-9.50076036, -4.51100042}, {-3.10745172,2.72118808} - }; + const std::vector> X_axpy = L; + std::vector> Y_axpy = R; + + const std::vector> output_axpy_op = {{-1.04187063, 2.57843137}, + {4.36025552, 5.18794423}, + {-11.42394822, -4.75054595}, + {6.50906977, 8.24708166}, + {1.46770822, -2.33007099}, + {4.53333119, -8.23539294}, + {-1.41856581, -5.83818462}, + {1.10278588, -2.14220619}, + {2.67173827, -3.46052805}, + {8.79651547, 4.55949747}, + {-9.50076036, -4.51100042}, + {-3.10745172, 2.72118808}}; // (5) for test scal_op (x = alpha * x) const std::complex alpha_scal{-1.5, -2.5}; - std::vector > X_scal = L; - - const std::vector >output_scal_op = { - {-0.87403307, 2.74844882}, {4.38946833, 4.90517834}, - {-11.33678411,-4.63461167}, {3.14637745, 8.26155093}, - {1.47152165, -2.17179400}, {4.29982696, -8.23346099}, - {-1.66408648, -6.48503924}, {-0.47976752,-4.85136318}, - {2.83788058, -3.52331200}, {6.62574354, 4.07239273}, - {-10.91333952,-5.05628302}, {-2.97811808, 2.72623036} - }; - + std::vector> X_scal = L; + + const std::vector> output_scal_op = {{-0.87403307, 2.74844882}, + {4.38946833, 4.90517834}, + {-11.33678411, -4.63461167}, + {3.14637745, 8.26155093}, + {1.47152165, -2.17179400}, + {4.29982696, -8.23346099}, + {-1.66408648, -6.48503924}, + {-0.47976752, -4.85136318}, + {2.83788058, -3.52331200}, + {6.62574354, 4.07239273}, + {-10.91333952, -5.05628302}, + {-2.97811808, 2.72623036}}; // (6) for test gemv_op ( y = alpha * op(A) * x + beta * y ) - const std::vector > A_gemv = { // 2 * 3 - {-0.87403307,2.74844882}, {4.38946833,4.90517834}, - {-11.33678411,-4.63461167}, {3.14637745,8.26155093}, - {1.47152165,-2.17179400}, {4.29982696,-8.23346099} - }; - - const std::vector > X_gemv = { // 2 * 1 - {-0.87403307,2.74844882}, {4.38946833,4.90517834} - }; - - std::vector > Y_gemv = { // 3 * 1 - {1.47152165,-2.17179400}, {4.29982696,-8.23346099}, - {3.14637745,8.26155093} - }; - - std::vector > Y_test_gemv = { // 3 * 1 - {1.47152165,-2.17179400}, {4.29982696,-8.23346099}, - {3.14637745,8.26155093} - }; - + const std::vector> A_gemv = {// 2 * 3 + {-0.87403307, 2.74844882}, + {4.38946833, 4.90517834}, + {-11.33678411, -4.63461167}, + {3.14637745, 8.26155093}, + {1.47152165, -2.17179400}, + {4.29982696, -8.23346099}}; + + const std::vector> X_gemv = {// 2 * 1 + {-0.87403307, 2.74844882}, + {4.38946833, 4.90517834}}; + + std::vector> Y_gemv = {// 3 * 1 + {1.47152165, -2.17179400}, + {4.29982696, -8.23346099}, + {3.14637745, 8.26155093}}; + + std::vector> Y_test_gemv = {// 3 * 1 + {1.47152165, -2.17179400}, + {4.29982696, -8.23346099}, + {3.14637745, 8.26155093}}; }; // template -// FPTYPE zdot_real(const int &dim, const std::complex* psi_L, const std::complex* psi_R, const psi::AbacusDevice_t device = psi::CpuDevice, const bool reduce = true); +// FPTYPE zdot_real(const int &dim, const std::complex* psi_L, const std::complex* psi_R, const +// psi::AbacusDevice_t device = psi::CpuDevice, const bool reduce = true); TEST_F(TestModuleHsolverMathKernel, zdot_real_op_cpu) { - double result = zdot_real_cpu_op()(cpu_ctx, dim, psi_L.data(), psi_R.data(), false); - EXPECT_LT(fabs(result - expected_result), 1e-12); + double result = zdot_real_cpu_op()(cpu_ctx, dim, psi_L.data(), psi_R.data(), false); + EXPECT_LT(fabs(result - expected_result), 1e-12); } TEST_F(TestModuleHsolverMathKernel, vector_div_constant_op_cpu) { - std::vector > output(input.size()); - vector_div_constant_op_cpu()(cpu_ctx, dim, output.data(), input.data(), constant); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_div_constant_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_div_constant_op[i].real()), 1e-8); - } + std::vector> output(input.size()); + vector_div_constant_op_cpu()(cpu_ctx, dim, output.data(), input.data(), constant); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_div_constant_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_div_constant_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, vector_mul_vector_op_cpu) { - std::vector > output(input.size()); - vector_mul_vector_op_cpu()(cpu_ctx, dim, output.data(), input.data(), input_double.data()); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_mul_vector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_mul_vector_op[i].real()), 1e-8); - } + std::vector> output(input.size()); + vector_mul_vector_op_cpu()(cpu_ctx, dim, output.data(), input.data(), input_double.data()); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_mul_vector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_mul_vector_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, vector_div_vector_op_cpu) { - std::vector > output(input.size()); - vector_div_vector_op_cpu()(cpu_ctx, dim, output.data(), input.data(), input_double.data()); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_div_vector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_div_vector_op[i].real()), 1e-8); - } + std::vector> output(input.size()); + vector_div_vector_op_cpu()(cpu_ctx, dim, output.data(), input.data(), input_double.data()); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_div_vector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_div_vector_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, constantvector_addORsub_constantVector_op_cpu) { - std::vector > output(input.size()); - constantvector_addORsub_constantVector_op_cpu()(cpu_ctx, dim, output.data(), input1.data(), constant1, input2.data(), constant2); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_constantvector_addORsub_constantVector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_constantvector_addORsub_constantVector_op[i].real()), 1e-8); - } + std::vector> output(input.size()); + constantvector_addORsub_constantVector_op_cpu()(cpu_ctx, + dim, + output.data(), + input1.data(), + constant1, + input2.data(), + constant2); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_constantvector_addORsub_constantVector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_constantvector_addORsub_constantVector_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, axpy_op_cpu) { - axpy_op_cpu()(cpu_ctx, dim, &alpha_axpy, X_axpy.data(), 1, Y_axpy.data(), 1); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(Y_axpy[i].imag() - output_axpy_op[i].imag()), 1e-8); - EXPECT_LT(fabs(Y_axpy[i].real() - output_axpy_op[i].real()), 1e-8); - } + axpy_op_cpu()(cpu_ctx, dim, &alpha_axpy, X_axpy.data(), 1, Y_axpy.data(), 1); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(Y_axpy[i].imag() - output_axpy_op[i].imag()), 1e-8); + EXPECT_LT(fabs(Y_axpy[i].real() - output_axpy_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, scal_op_cpu) { - scal_op_cpu()(cpu_ctx, dim, &alpha_scal, X_scal.data(), 1); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(X_scal[i].imag() - output_scal_op[i].imag()), 1e-8); - EXPECT_LT(fabs(X_scal[i].real() - output_scal_op[i].real()), 1e-8); - } + scal_op_cpu()(cpu_ctx, dim, &alpha_scal, X_scal.data(), 1); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(X_scal[i].imag() - output_scal_op[i].imag()), 1e-8); + EXPECT_LT(fabs(X_scal[i].real() - output_scal_op[i].real()), 1e-8); + } } TEST_F(TestModuleHsolverMathKernel, gemv_op_cpu) { - gemv_op_cpu()(cpu_ctx, 'C', 2, 3, &ModuleBase::ONE, A_gemv.data(), 2, X_gemv.data(), 1, &ModuleBase::ONE, Y_gemv.data(), 1); - char trans = 'C'; - int inc = 1; - int row = 2; - int col = 3; - zgemv_(&trans, &row, &col, &ModuleBase::ONE, A_gemv.data(), &row, X_gemv.data(), &inc, &ModuleBase::ONE, Y_test_gemv.data(), &inc); - for (int i = 0; i < Y_gemv.size(); i++) - { - EXPECT_LT(fabs(Y_gemv[i].imag() - Y_test_gemv[i].imag()), 1e-12); - EXPECT_LT(fabs(Y_gemv[i].real() - Y_test_gemv[i].real()), 1e-12); - } + gemv_op_cpu()(cpu_ctx, + 'C', + 2, + 3, + &ModuleBase::ONE, + A_gemv.data(), + 2, + X_gemv.data(), + 1, + &ModuleBase::ONE, + Y_gemv.data(), + 1); + char trans = 'C'; + int inc = 1; + int row = 2; + int col = 3; + zgemv_(&trans, + &row, + &col, + &ModuleBase::ONE, + A_gemv.data(), + &row, + X_gemv.data(), + &inc, + &ModuleBase::ONE, + Y_test_gemv.data(), + &inc); + for (int i = 0; i < Y_gemv.size(); i++) + { + EXPECT_LT(fabs(Y_gemv[i].imag() - Y_test_gemv[i].imag()), 1e-12); + EXPECT_LT(fabs(Y_gemv[i].real() - Y_test_gemv[i].real()), 1e-12); + } } - #if __UT_USE_CUDA || __UT_USE_ROCM TEST_F(TestModuleHsolverMathKernel, zdot_real_op_gpu) { - std::complex* psi_L_dev = NULL, * psi_R_dev = NULL; - resize_memory_op()(gpu_ctx, psi_L_dev, psi_L.size()); - resize_memory_op()(gpu_ctx, psi_R_dev, psi_R.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, psi_L_dev, psi_L.data(), psi_L.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, psi_R_dev, psi_R.data(), psi_R.size()); - double result = zdot_real_gpu_op()(gpu_ctx, dim, psi_L_dev, psi_R_dev, false); - EXPECT_LT(fabs(result - expected_result), 1e-12); - delete_memory_op()(gpu_ctx, psi_L_dev); - delete_memory_op()(gpu_ctx, psi_R_dev); + std::complex*psi_L_dev = NULL, *psi_R_dev = NULL; + resize_memory_op()(gpu_ctx, psi_L_dev, psi_L.size()); + resize_memory_op()(gpu_ctx, psi_R_dev, psi_R.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, psi_L_dev, psi_L.data(), psi_L.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, psi_R_dev, psi_R.data(), psi_R.size()); + double result = zdot_real_gpu_op()(gpu_ctx, dim, psi_L_dev, psi_R_dev, false); + EXPECT_LT(fabs(result - expected_result), 1e-12); + delete_memory_op()(gpu_ctx, psi_L_dev); + delete_memory_op()(gpu_ctx, psi_R_dev); } TEST_F(TestModuleHsolverMathKernel, vector_div_constant_op_gpu) { - // in CPU - std::vector > output(input.size()); - // in GPU - std::complex* input_dev = NULL; - std::complex* output_dev = NULL; - resize_memory_op()(gpu_ctx, input_dev, input.size()); - resize_memory_op()(gpu_ctx, output_dev, input.size()); - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); - // run - vector_div_constant_op_gpu()(gpu_ctx, dim, output_dev, input_dev, constant); - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); - - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_div_constant_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_div_constant_op[i].real()), 1e-8); - } - delete_memory_op()(gpu_ctx, input_dev); - delete_memory_op()(gpu_ctx, output_dev); + // in CPU + std::vector> output(input.size()); + // in GPU + std::complex* input_dev = NULL; + std::complex* output_dev = NULL; + resize_memory_op()(gpu_ctx, input_dev, input.size()); + resize_memory_op()(gpu_ctx, output_dev, input.size()); + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); + // run + vector_div_constant_op_gpu()(gpu_ctx, dim, output_dev, input_dev, constant); + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); + + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_div_constant_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_div_constant_op[i].real()), 1e-8); + } + delete_memory_op()(gpu_ctx, input_dev); + delete_memory_op()(gpu_ctx, output_dev); } TEST_F(TestModuleHsolverMathKernel, vector_mul_vector_op_gpu) { - // in CPU - std::vector > output(input.size()); - - // in GPU - std::complex* input_dev = NULL; - double* input_double_dev = NULL; - std::complex* output_dev = NULL; - - // resize memory for values - resize_memory_op()(gpu_ctx, input_dev, input.size()); - resize_memory_op_double()(gpu_ctx, input_double_dev, input.size()); - resize_memory_op()(gpu_ctx, output_dev, input.size()); - - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); - synchronize_memory_op_double()(gpu_ctx, cpu_ctx, input_double_dev, input_double.data(), input.size()); - - // run - vector_mul_vector_op_gpu()(gpu_ctx, dim, output_dev, input_dev, input_double_dev); - - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); - - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_mul_vector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_mul_vector_op[i].real()), 1e-8); - } + // in CPU + std::vector> output(input.size()); - delete_memory_op()(gpu_ctx, input_dev); - delete_memory_op_double()(gpu_ctx, input_double_dev); - delete_memory_op()(gpu_ctx, output_dev); -} + // in GPU + std::complex* input_dev = NULL; + double* input_double_dev = NULL; + std::complex* output_dev = NULL; + // resize memory for values + resize_memory_op()(gpu_ctx, input_dev, input.size()); + resize_memory_op_double()(gpu_ctx, input_double_dev, input.size()); + resize_memory_op()(gpu_ctx, output_dev, input.size()); -TEST_F(TestModuleHsolverMathKernel, vector_div_vector_op_gpu) -{ - // in CPU - std::vector > output(input.size()); - - // in GPU - std::complex* input_dev = NULL; - double* input_double_dev = NULL; - std::complex* output_dev = NULL; - - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, input_dev, input.size()); - resize_memory_op_double()(gpu_ctx, input_double_dev, input.size()); - resize_memory_op()(gpu_ctx, output_dev, input.size()); - - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); - synchronize_memory_op_double()(gpu_ctx, cpu_ctx, input_double_dev, input_double.data(), input.size()); - - // run - vector_div_vector_op_gpu()(gpu_ctx, dim, output_dev, input_dev, input_double_dev); - - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); - - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_vector_div_vector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_vector_div_vector_op[i].real()), 1e-8); - } + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); + synchronize_memory_op_double()(gpu_ctx, cpu_ctx, input_double_dev, input_double.data(), input.size()); - delete_memory_op()(gpu_ctx, input_dev); - delete_memory_op_double()(gpu_ctx, input_double_dev); - delete_memory_op()(gpu_ctx, output_dev); -} + // run + vector_mul_vector_op_gpu()(gpu_ctx, dim, output_dev, input_dev, input_double_dev); + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); -TEST_F(TestModuleHsolverMathKernel, constantvector_addORsub_constantVector_op_gpu) + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_mul_vector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_mul_vector_op[i].real()), 1e-8); + } + + delete_memory_op()(gpu_ctx, input_dev); + delete_memory_op_double()(gpu_ctx, input_double_dev); + delete_memory_op()(gpu_ctx, output_dev); +} + +TEST_F(TestModuleHsolverMathKernel, vector_div_vector_op_gpu) { - // in CPU - std::vector > output(input.size()); + // in CPU + std::vector> output(input.size()); - // in GPU - std::complex* input1_dev = NULL; - std::complex* input2_dev = NULL; - std::complex* output_dev = NULL; + // in GPU + std::complex* input_dev = NULL; + double* input_double_dev = NULL; + std::complex* output_dev = NULL; - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, input1_dev, input.size()); - resize_memory_op()(gpu_ctx, input2_dev, input.size()); - resize_memory_op()(gpu_ctx, output_dev, input.size()); + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, input_dev, input.size()); + resize_memory_op_double()(gpu_ctx, input_double_dev, input.size()); + resize_memory_op()(gpu_ctx, output_dev, input.size()); - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, input1_dev, input1.data(), input.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, input2_dev, input2.data(), input.size()); + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, input_dev, input.data(), input.size()); + synchronize_memory_op_double()(gpu_ctx, cpu_ctx, input_double_dev, input_double.data(), input.size()); - // run - constantvector_addORsub_constantVector_op_gpu()(gpu_ctx, dim, output_dev, input1_dev, constant1, input2_dev, constant2); + // run + vector_div_vector_op_gpu()(gpu_ctx, dim, output_dev, input_dev, input_double_dev); - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(output[i].imag() - output_constantvector_addORsub_constantVector_op[i].imag()), 1e-8); - EXPECT_LT(fabs(output[i].real() - output_constantvector_addORsub_constantVector_op[i].real()), 1e-8); - } + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_vector_div_vector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_vector_div_vector_op[i].real()), 1e-8); + } - delete_memory_op()(gpu_ctx, input1_dev); - delete_memory_op()(gpu_ctx, input2_dev); - delete_memory_op()(gpu_ctx, output_dev); + delete_memory_op()(gpu_ctx, input_dev); + delete_memory_op_double()(gpu_ctx, input_double_dev); + delete_memory_op()(gpu_ctx, output_dev); } - -TEST_F(TestModuleHsolverMathKernel, axpy_op_gpu) +TEST_F(TestModuleHsolverMathKernel, constantvector_addORsub_constantVector_op_gpu) { - // in GPU - std::complex* X_axpy_dev = NULL; - std::complex* Y_axpy_dev = NULL; - - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, X_axpy_dev, X_axpy.size()); - resize_memory_op()(gpu_ctx, Y_axpy_dev, Y_axpy.size()); - - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, X_axpy_dev, X_axpy.data(), X_axpy.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, Y_axpy_dev, Y_axpy.data(), Y_axpy.size()); - - // run - hsolver::createBLAShandle(); - axpy_op_gpu()(gpu_ctx, dim, &alpha_axpy, X_axpy_dev, 1, Y_axpy_dev, 1); - hsolver::destoryBLAShandle(); - - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, Y_axpy.data(), Y_axpy_dev, Y_axpy.size()); + // in CPU + std::vector> output(input.size()); + + // in GPU + std::complex* input1_dev = NULL; + std::complex* input2_dev = NULL; + std::complex* output_dev = NULL; + + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, input1_dev, input.size()); + resize_memory_op()(gpu_ctx, input2_dev, input.size()); + resize_memory_op()(gpu_ctx, output_dev, input.size()); + + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, input1_dev, input1.data(), input.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, input2_dev, input2.data(), input.size()); + + // run + constantvector_addORsub_constantVector_op_gpu()(gpu_ctx, + dim, + output_dev, + input1_dev, + constant1, + input2_dev, + constant2); + + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, output.data(), output_dev, output.size()); + + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(output[i].imag() - output_constantvector_addORsub_constantVector_op[i].imag()), 1e-8); + EXPECT_LT(fabs(output[i].real() - output_constantvector_addORsub_constantVector_op[i].real()), 1e-8); + } + delete_memory_op()(gpu_ctx, input1_dev); + delete_memory_op()(gpu_ctx, input2_dev); + delete_memory_op()(gpu_ctx, output_dev); +} - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(Y_axpy[i].imag() - output_axpy_op[i].imag()), 1e-8); - EXPECT_LT(fabs(Y_axpy[i].real() - output_axpy_op[i].real()), 1e-8); - } +TEST_F(TestModuleHsolverMathKernel, axpy_op_gpu) +{ + // in GPU + std::complex* X_axpy_dev = NULL; + std::complex* Y_axpy_dev = NULL; + + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, X_axpy_dev, X_axpy.size()); + resize_memory_op()(gpu_ctx, Y_axpy_dev, Y_axpy.size()); + + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, X_axpy_dev, X_axpy.data(), X_axpy.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, Y_axpy_dev, Y_axpy.data(), Y_axpy.size()); + + // run + hsolver::createBLAShandle(); + axpy_op_gpu()(gpu_ctx, dim, &alpha_axpy, X_axpy_dev, 1, Y_axpy_dev, 1); + hsolver::destoryBLAShandle(); + + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, Y_axpy.data(), Y_axpy_dev, Y_axpy.size()); + + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(Y_axpy[i].imag() - output_axpy_op[i].imag()), 1e-8); + EXPECT_LT(fabs(Y_axpy[i].real() - output_axpy_op[i].real()), 1e-8); + } - delete_memory_op()(gpu_ctx, X_axpy_dev); - delete_memory_op()(gpu_ctx, Y_axpy_dev); + delete_memory_op()(gpu_ctx, X_axpy_dev); + delete_memory_op()(gpu_ctx, Y_axpy_dev); } TEST_F(TestModuleHsolverMathKernel, scal_op_gpu) { - // in GPU - std::complex* X_scal_dev = NULL; + // in GPU + std::complex* X_scal_dev = NULL; - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, X_scal_dev, X_scal.size()); - + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, X_scal_dev, X_scal.size()); - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, X_scal_dev, X_scal.data(), X_scal.size()); + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, X_scal_dev, X_scal.data(), X_scal.size()); - // run - hsolver::createBLAShandle(); - scal_op_gpu()(gpu_ctx, dim, &alpha_scal, X_scal_dev, 1); - hsolver::destoryBLAShandle(); + // run + hsolver::createBLAShandle(); + scal_op_gpu()(gpu_ctx, dim, &alpha_scal, X_scal_dev, 1); + hsolver::destoryBLAShandle(); - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, X_scal.data(), X_scal_dev, X_scal.size()); + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, X_scal.data(), X_scal_dev, X_scal.size()); - - for (int i = 0; i < input.size(); i++) - { - EXPECT_LT(fabs(X_scal[i].imag() - output_scal_op[i].imag()), 1e-8); - EXPECT_LT(fabs(X_scal[i].real() - output_scal_op[i].real()), 1e-8); - } - delete_memory_op()(gpu_ctx, X_scal_dev); + for (int i = 0; i < input.size(); i++) + { + EXPECT_LT(fabs(X_scal[i].imag() - output_scal_op[i].imag()), 1e-8); + EXPECT_LT(fabs(X_scal[i].real() - output_scal_op[i].real()), 1e-8); + } + delete_memory_op()(gpu_ctx, X_scal_dev); } TEST_F(TestModuleHsolverMathKernel, gemv_op_gpu) { - // in GPU - std::complex* A_gemv_dev = NULL; - std::complex* X_gemv_dev = NULL; - std::complex* Y_gemv_dev = NULL; - - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, A_gemv_dev, A_gemv.size()); - resize_memory_op()(gpu_ctx, X_gemv_dev, X_gemv.size()); - resize_memory_op()(gpu_ctx, Y_gemv_dev, Y_gemv.size()); - - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, A_gemv_dev, A_gemv.data(), A_gemv.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, X_gemv_dev, X_gemv.data(), X_gemv.size()); - synchronize_memory_op()(gpu_ctx, cpu_ctx, Y_gemv_dev, Y_gemv.data(), Y_gemv.size()); - - // run - hsolver::createBLAShandle(); - gemv_op_gpu()(gpu_ctx, 'C', 2, 3, &ModuleBase::ONE, A_gemv_dev, 2, X_gemv_dev, 1, &ModuleBase::ONE, Y_gemv_dev, 1); - hsolver::destoryBLAShandle(); - // syn the output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, Y_gemv.data(), Y_gemv_dev, Y_gemv.size()); - - // cal right answer: Y_test_gemv - char trans = 'C'; - int inc = 1; - int row = 2; - int col = 3; - zgemv_(&trans, &row, &col, &ModuleBase::ONE, A_gemv.data(), &row, X_gemv.data(), &inc, &ModuleBase::ONE, Y_test_gemv.data(), &inc); - - - for (int i = 0; i < Y_gemv.size(); i++) - { - EXPECT_LT(fabs(Y_gemv[i].imag() - Y_test_gemv[i].imag()), 1e-12); - EXPECT_LT(fabs(Y_gemv[i].real() - Y_test_gemv[i].real()), 1e-12); - } + // in GPU + std::complex* A_gemv_dev = NULL; + std::complex* X_gemv_dev = NULL; + std::complex* Y_gemv_dev = NULL; + + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, A_gemv_dev, A_gemv.size()); + resize_memory_op()(gpu_ctx, X_gemv_dev, X_gemv.size()); + resize_memory_op()(gpu_ctx, Y_gemv_dev, Y_gemv.size()); + + // syn the input data in CPU to GPU + synchronize_memory_op()(gpu_ctx, cpu_ctx, A_gemv_dev, A_gemv.data(), A_gemv.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, X_gemv_dev, X_gemv.data(), X_gemv.size()); + synchronize_memory_op()(gpu_ctx, cpu_ctx, Y_gemv_dev, Y_gemv.data(), Y_gemv.size()); + + // run + hsolver::createBLAShandle(); + gemv_op_gpu()(gpu_ctx, 'C', 2, 3, &ModuleBase::ONE, A_gemv_dev, 2, X_gemv_dev, 1, &ModuleBase::ONE, Y_gemv_dev, 1); + hsolver::destoryBLAShandle(); + // syn the output data in GPU to CPU + synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, Y_gemv.data(), Y_gemv_dev, Y_gemv.size()); + + // cal right answer: Y_test_gemv + char trans = 'C'; + int inc = 1; + int row = 2; + int col = 3; + zgemv_(&trans, + &row, + &col, + &ModuleBase::ONE, + A_gemv.data(), + &row, + X_gemv.data(), + &inc, + &ModuleBase::ONE, + Y_test_gemv.data(), + &inc); + + for (int i = 0; i < Y_gemv.size(); i++) + { + EXPECT_LT(fabs(Y_gemv[i].imag() - Y_test_gemv[i].imag()), 1e-12); + EXPECT_LT(fabs(Y_gemv[i].real() - Y_test_gemv[i].real()), 1e-12); + } - delete_memory_op()(gpu_ctx, A_gemv_dev); - delete_memory_op()(gpu_ctx, X_gemv_dev); - delete_memory_op()(gpu_ctx, Y_gemv_dev); + delete_memory_op()(gpu_ctx, A_gemv_dev); + delete_memory_op()(gpu_ctx, X_gemv_dev); + delete_memory_op()(gpu_ctx, Y_gemv_dev); } +TEST_F(TestModuleHsolverMathKernel, matrixSetToAnother_op_gpu) +{ + // const std::vector > expect_result = { + // {-0.11893203,-0.13492526}, {-0.40314756, 0.07734553}, {0.06892412, 0.14837423}, {0.0, 0.0}, + // {0.61158728, -0.45754102}, {-0.54274745,-0.09682102}, {0.30232967, 0.49411249}, {0.0, 0.0} + // }; + + const std::vector> A = {{-0.11893203, -0.13492526}, + {-0.40314756, 0.07734553}, + {0.06892412, 0.14837423}, + {0.61158728, -0.45754102}, + {-0.54274745, -0.09682102}, + {0.30232967, 0.49411249}}; + + const std::vector> B(8); + + int n = 2; + int LDA = 3; + int LDB = 4; + + std::complex* device_A = nullptr; + psi::memory::resize_memory_op, psi::DEVICE_GPU>()(gpu_ctx, device_A, A.size()); + psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_CPU>()(gpu_ctx, + cpu_ctx, + device_A, + A.data(), + A.size()); + + std::complex* device_B = nullptr; + psi::memory::resize_memory_op, psi::DEVICE_GPU>()(gpu_ctx, device_B, B.size()); + psi::memory::synchronize_memory_op, psi::DEVICE_GPU, psi::DEVICE_CPU>()(gpu_ctx, + cpu_ctx, + device_B, + B.data(), + B.size()); + + // run + hsolver::matrixSetToAnother()(gpu_ctx, n, device_A, LDA, device_B, LDB); + + std::vector> B_gpu2cpu(8); + psi::memory::synchronize_memory_op, psi::DEVICE_CPU, psi::DEVICE_GPU>()(cpu_ctx, + gpu_ctx, + B_gpu2cpu.data(), + device_B, + B_gpu2cpu.size()); + + std::vector> B_cpu(8); + hsolver::matrixSetToAnother()(cpu_ctx, n, A.data(), LDA, B_cpu.data(), LDB); + + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 2; j++) + // { + // std::cout << B_gpu2cpu[i * 2 + j]; + // } + // std::cout << std::endl; + // } + // std::cout << std::endl; + + // for (int i = 0; i < 4; i++) + // { + // for (int j = 0; j < 2; j++) + // { + // std::cout << B_cpu[i * 2 + j]; + // } + // std::cout << std::endl; + // } + + for (int i = 0; i < B_cpu.size(); i++) + { + EXPECT_LT(fabs(B_gpu2cpu[i].imag() - B_cpu[i].imag()), 1e-12); + EXPECT_LT(fabs(B_gpu2cpu[i].real() - B_cpu[i].real()), 1e-12); + } - + delete_memory_op()(gpu_ctx, device_A); + delete_memory_op()(gpu_ctx, device_B); +} #endif // __UT_USE_CUDA || __UT_USE_ROCM @@ -550,7 +715,7 @@ TEST_F(TestModuleHsolverMathKernel, gemv_op_gpu) TEST_F(TestGelu, gelu_gpu_rocm) { std::vector gelu(nloc, 0.0); - + double * gelu_dev = NULL, * xx_dev = NULL; deepmd::malloc_device_memory_sync(gelu_dev, gelu); deepmd::malloc_device_memory_sync(xx_dev, xx); @@ -563,7 +728,7 @@ TEST_F(TestGelu, gelu_gpu_rocm) EXPECT_EQ(gelu.size(), expected_gelu.size()); for (int jj = 0; jj < gelu.size(); ++jj){ EXPECT_LT(fabs(gelu[jj] - expected_gelu[jj]) , 1e-5); - } + } } #endif // TENSORFLOW_USE_ROCM */ \ No newline at end of file