From 47fd75388d77314d5a09c9df0029c7f8d0b1915a Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:30:19 +0800 Subject: [PATCH 01/14] Add matrixTranspose_op & matrixSetToAnother_op --- source/module_hsolver/include/math_kernel.h | 31 +++++++ source/module_hsolver/src/cuda/math_kernel.cu | 84 +++++++++++++++++++ source/module_hsolver/src/math_kernel.cpp | 49 +++++++++++ 3 files changed, 164 insertions(+) diff --git a/source/module_hsolver/include/math_kernel.h b/source/module_hsolver/include/math_kernel.h index d343a0f3dcf..8e9e49d3960 100644 --- a/source/module_hsolver/include/math_kernel.h +++ b/source/module_hsolver/include/math_kernel.h @@ -6,6 +6,7 @@ #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) @@ -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/math_kernel.cu b/source/module_hsolver/src/cuda/math_kernel.cu index 3f8c04214c9..762c98e023f 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. + 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/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 From 9cdec40ab871cc46adfb8514614f924c18bf61ae Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:32:33 +0800 Subject: [PATCH 02/14] Add dngv_op & dngvd_op --- source/module_hsolver/include/dngvd_op.h | 73 +++++++ source/module_hsolver/src/cuda/dngvd_op.cu | 187 ++++++++++++++++++ source/module_hsolver/src/dngvd_op.cpp | 217 +++++++++++++++++++++ 3 files changed, 477 insertions(+) create mode 100644 source/module_hsolver/include/dngvd_op.h create mode 100644 source/module_hsolver/src/cuda/dngvd_op.cu create mode 100644 source/module_hsolver/src/dngvd_op.cpp diff --git a/source/module_hsolver/include/dngvd_op.h b/source/module_hsolver/include/dngvd_op.h new file mode 100644 index 00000000000..9ac15c27ed2 --- /dev/null +++ b/source/module_hsolver/include/dngvd_op.h @@ -0,0 +1,73 @@ +// 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_psi/include/memory.h" +#include "module_hsolver/include/math_kernel.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/src/cuda/dngvd_op.cu b/source/module_hsolver/src/cuda/dngvd_op.cu new file mode 100644 index 00000000000..10fc8fa10f9 --- /dev/null +++ b/source/module_hsolver/src/cuda/dngvd_op.cu @@ -0,0 +1,187 @@ +#include "module_hsolver/include/dngvd_op.h" + + +#include + +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; + cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col); + 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; + cudaMalloc((void**)&all_W, sizeof(double) * col); + + // prepare some values for cusolverDnZhegvd_bufferSize + cusolverDnHandle_t cusolverH; + cusolverDnCreate(&cusolverH); + int* devInfo; + cudaMalloc((void**)&devInfo, sizeof(int)); + + // calculate the sizes needed for pre-allocated buffer. + int lwork = 0; + 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; + cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork); + + // compute eigenvalues and eigenvectors. + 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); + + cudaDeviceSynchronize(); + + // get eigenvalues and eigenvectors. only m ! + cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice); + + cudaMemcpy(V, A_eigenvectors, sizeof(std::complex)*col*m, cudaMemcpyDeviceToDevice); + + matrixTranspose_op()(d, col, row, V, V); + + int info_gpu; + cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost); + assert(0 == info_gpu); + + // free the buffer + cudaFree(d_work); + // free resources and destroy + cudaFree(A_eigenvectors); + cudaFree(all_W); + cudaFree(devInfo); + 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; + cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col); + cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col); + + // transpose A to A_eigenvectors + matrixTranspose_op()(d, row, col, A, (std::complex*)A_eigenvectors); + // transpose B to transpose_B + matrixTranspose_op()(d, row, col, B, (std::complex*)transpose_B); + + // init all_W + double* all_W; + cudaMalloc((void**)&all_W, sizeof(double) * row); + + // prepare some values for cusolverDnZhegvd_bufferSize + cusolverDnHandle_t cusolverH; + cusolverDnCreate(&cusolverH); + int* devInfo; + cudaMalloc((void**)&devInfo, sizeof(int)); + + // calculate the sizes needed for pre-allocated buffer. + int lwork = 0; + 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; + cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork); + + // compute eigenvalues and eigenvectors. + 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); + + cudaDeviceSynchronize(); + + // get all eigenvalues and eigenvectors. + cudaMemcpy(W, all_W, sizeof(double) * row, cudaMemcpyDeviceToDevice); + matrixTranspose_op()(d, row, col, (std::complex*)A_eigenvectors, V); + + int info_gpu; + cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost); + assert(0 == info_gpu); + + // free the buffer + cudaFree(d_work); + // free resources and destroy + cudaFree(A_eigenvectors); + cudaFree(all_W); + cudaFree(devInfo); + cusolverDnDestroy(cusolverH); + destoryBLAShandle(); +} + +} // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/src/dngvd_op.cpp b/source/module_hsolver/src/dngvd_op.cpp new file mode 100644 index 00000000000..f0e699cf773 --- /dev/null +++ b/source/module_hsolver/src/dngvd_op.cpp @@ -0,0 +1,217 @@ +#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 From 94e86f36823a13c297fdf4c7ca652392d298af98 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:33:11 +0800 Subject: [PATCH 03/14] Add dngv UT --- source/module_hsolver/CMakeLists.txt | 6 +- source/module_hsolver/test/CMakeLists.txt | 18 + .../module_hsolver/test/math_dngvd_test.cpp | 567 ++++++++++++++++++ 3 files changed, 590 insertions(+), 1 deletion(-) create mode 100644 source/module_hsolver/test/math_dngvd_test.cpp 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/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..17ab54afc3d --- /dev/null +++ b/source/module_hsolver/test/math_dngvd_test.cpp @@ -0,0 +1,567 @@ +#include +#include +#include +#include +#include +#include +#include "module_psi/include/memory.h" +#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" + + +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 From 309161e4a58eaf390b9a50c7154a81587520ca9b Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:33:57 +0800 Subject: [PATCH 04/14] Add matrixSetToAnother_op_gpu UT --- .../module_hsolver/test/math_kernel_test.cpp | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/source/module_hsolver/test/math_kernel_test.cpp b/source/module_hsolver/test/math_kernel_test.cpp index f24e285f7a7..c6297d522c2 100644 --- a/source/module_hsolver/test/math_kernel_test.cpp +++ b/source/module_hsolver/test/math_kernel_test.cpp @@ -542,6 +542,97 @@ TEST_F(TestModuleHsolverMathKernel, gemv_op_gpu) +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 = { + {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} + }; + + 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 = { + {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}, + }; + 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 = { + {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::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 From 62f9e724187d9c6a15880df68027dd712d639c8f Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:37:44 +0800 Subject: [PATCH 05/14] Support GPU for diagH_LAPACK & diagH_subspace --- source/module_hsolver/diago_iter_assist.cpp | 271 +++++++++++++------- source/module_hsolver/diago_iter_assist.h | 11 +- 2 files changed, 194 insertions(+), 88 deletions(-) diff --git a/source/module_hsolver/diago_iter_assist.cpp b/source/module_hsolver/diago_iter_assist.cpp index be44920b1a1..18ebd51cc28 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,61 @@ 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, + std::complex* vcc) +{ + ModuleBase::TITLE("DiagoIterAssist", "LAPACK_subspace"); + ModuleBase::timer::tick("DiagoIterAssist", "LAPACK_subspace"); + + Device* ctx = {}; + + const bool all_eigenvalues = (nstart == nbands); + + 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 + ); + } + + + 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 From 845cd487530c5a1376ea475c9a1a9efd197cea10 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:38:57 +0800 Subject: [PATCH 06/14] Move eigenvalue_in to GPU memory --- source/module_hsolver/diago_cg.cpp | 48 +++++++++++++++++++----------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/source/module_hsolver/diago_cg.cpp b/source/module_hsolver/diago_cg.cpp index 1f2a39698ca..b109b1e070f 100644 --- a/source/module_hsolver/diago_cg.cpp +++ b/source/module_hsolver/diago_cg.cpp @@ -683,23 +683,37 @@ 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 + #if defined(__CUDA) || defined(__ROCM) + const psi::DEVICE_CPU * cpu_ctx = {}; + const psi::DEVICE_GPU * gpu_ctx = {}; + + FPTYPE* eigenvalue_gpu = nullptr; + psi::memory::resize_memory_op()(gpu_ctx, eigenvalue_gpu, psi.get_nbands()); + + // set eigenvalue_in value to eigenvalue_gpu + psi::memory::synchronize_memory_op()( + gpu_ctx, + cpu_ctx, + eigenvalue_gpu, + eigenvalue_in, + psi.get_nbands() + ); + + // run in GPU + DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_gpu); + + // set eigenvalue_in value to eigenvalue_gpu + psi::memory::synchronize_memory_op()( + cpu_ctx, + gpu_ctx, + eigenvalue_in, + eigenvalue_gpu, + psi.get_nbands() + ); + #else + // run in CPU + DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_in); + #endif } DiagoIterAssist::avg_iter += 1.0; From 85396eb657f20ffc83ae06c47e8e2e68e1f94709 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 10:42:30 +0800 Subject: [PATCH 07/14] Update lapack_connector.h --- source/module_base/lapack_connector.h | 71 +++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) 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. From f103ec76b69ac881644018c214e8fd1abedb4b10 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 15:09:25 +0800 Subject: [PATCH 08/14] Format the codes in math_kernel_test.cpp --- .../module_hsolver/test/math_kernel_test.cpp | 19 +++---------------- 1 file changed, 3 insertions(+), 16 deletions(-) diff --git a/source/module_hsolver/test/math_kernel_test.cpp b/source/module_hsolver/test/math_kernel_test.cpp index c6297d522c2..31a8a335d5c 100644 --- a/source/module_hsolver/test/math_kernel_test.cpp +++ b/source/module_hsolver/test/math_kernel_test.cpp @@ -554,10 +554,7 @@ TEST_F(TestModuleHsolverMathKernel, matrixSetToAnother_op_gpu) {0.61158728, -0.45754102}, {-0.54274745,-0.09682102}, {0.30232967, 0.49411249} }; - const std::vector > B = { - {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} - }; + const std::vector > B(8); int n = 2; int LDA = 3; @@ -582,20 +579,10 @@ TEST_F(TestModuleHsolverMathKernel, matrixSetToAnother_op_gpu) device_B, LDB); - std::vector > B_gpu2cpu = { - {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}, - }; + 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 = { - {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}, - }; + std::vector > B_cpu(8); hsolver::matrixSetToAnother()( cpu_ctx, n, From 758b960332298e1089ced0a3443c14bca49ef5a8 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 15:25:34 +0800 Subject: [PATCH 09/14] Move e memory transfer from diago_cg to subspace --- source/module_hsolver/diago_cg.cpp | 30 ---------- source/module_hsolver/diago_iter_assist.cpp | 63 ++++++++++++++------- 2 files changed, 44 insertions(+), 49 deletions(-) diff --git a/source/module_hsolver/diago_cg.cpp b/source/module_hsolver/diago_cg.cpp index b109b1e070f..cc8a96f0c02 100644 --- a/source/module_hsolver/diago_cg.cpp +++ b/source/module_hsolver/diago_cg.cpp @@ -683,37 +683,7 @@ void DiagoCG::diag(hamilt::Hamilt *phm_in, psi:: { if(DiagoIterAssist::need_subspace || ntry > 0) { - #if defined(__CUDA) || defined(__ROCM) - const psi::DEVICE_CPU * cpu_ctx = {}; - const psi::DEVICE_GPU * gpu_ctx = {}; - - FPTYPE* eigenvalue_gpu = nullptr; - psi::memory::resize_memory_op()(gpu_ctx, eigenvalue_gpu, psi.get_nbands()); - - // set eigenvalue_in value to eigenvalue_gpu - psi::memory::synchronize_memory_op()( - gpu_ctx, - cpu_ctx, - eigenvalue_gpu, - eigenvalue_in, - psi.get_nbands() - ); - - // run in GPU - DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_gpu); - - // set eigenvalue_in value to eigenvalue_gpu - psi::memory::synchronize_memory_op()( - cpu_ctx, - gpu_ctx, - eigenvalue_in, - eigenvalue_gpu, - psi.get_nbands() - ); - #else - // run in CPU DiagoIterAssist::diagH_subspace(phm_in, psi, psi, eigenvalue_in); - #endif } DiagoIterAssist::avg_iter += 1.0; diff --git a/source/module_hsolver/diago_iter_assist.cpp b/source/module_hsolver/diago_iter_assist.cpp index 18ebd51cc28..8f74a128b3f 100644 --- a/source/module_hsolver/diago_iter_assist.cpp +++ b/source/module_hsolver/diago_iter_assist.cpp @@ -529,38 +529,63 @@ void DiagoIterAssist::diagH_LAPACK( 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, - vcc - ); + 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, - vcc - ); + 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 + ); +#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; From 121de838d4f09923036fe3a91c7e02a4795e7cd3 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Tue, 15 Nov 2022 15:30:00 +0800 Subject: [PATCH 10/14] Add delete_memory_op for e_gpu in diagH_LAPACK --- source/module_hsolver/diago_iter_assist.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/module_hsolver/diago_iter_assist.cpp b/source/module_hsolver/diago_iter_assist.cpp index 8f74a128b3f..d122386bf22 100644 --- a/source/module_hsolver/diago_iter_assist.cpp +++ b/source/module_hsolver/diago_iter_assist.cpp @@ -519,7 +519,7 @@ void DiagoIterAssist::diagH_LAPACK( const std::complex* hcc, const std::complex* scc, const int ldh, // nstart - FPTYPE *e, + FPTYPE *e, // always in CPU std::complex* vcc) { ModuleBase::TITLE("DiagoIterAssist", "LAPACK_subspace"); @@ -535,7 +535,6 @@ void DiagoIterAssist::diagH_LAPACK( 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, @@ -568,6 +567,7 @@ void DiagoIterAssist::diagH_LAPACK( e_gpu, nbands ); + psi::memory::delete_memory_op()(gpu_ctx, e_gpu); #else if (all_eigenvalues) From b9022368ae686c2b2039cffa388810c8a7fa2443 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Wed, 16 Nov 2022 11:36:32 +0800 Subject: [PATCH 11/14] Add error checking macro for cuBlas/cuda API calls --- source/module_hsolver/include/math_kernel.h | 8 +- source/module_hsolver/src/cuda/dngvd_op.cu | 125 ++++++++++++------ source/module_hsolver/src/cuda/math_kernel.cu | 2 +- 3 files changed, 89 insertions(+), 46 deletions(-) diff --git a/source/module_hsolver/include/math_kernel.h b/source/module_hsolver/include/math_kernel.h index 8e9e49d3960..1190c043018 100644 --- a/source/module_hsolver/include/math_kernel.h +++ b/source/module_hsolver/include/math_kernel.h @@ -13,9 +13,9 @@ #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"; @@ -44,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); } } diff --git a/source/module_hsolver/src/cuda/dngvd_op.cu b/source/module_hsolver/src/cuda/dngvd_op.cu index 10fc8fa10f9..6e5dbc479ca 100644 --- a/source/module_hsolver/src/cuda/dngvd_op.cu +++ b/source/module_hsolver/src/cuda/dngvd_op.cu @@ -1,8 +1,50 @@ #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 { @@ -17,29 +59,28 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, std::complex* V) { createBLAShandle(); + // init A_eigenvectors, transpose_B and all_W double2 *A_eigenvectors, *transpose_B; - cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col); - cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col); + 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; - cudaMalloc((void**)&all_W, sizeof(double) * col); + checkCudaErrors( cudaMalloc((void**)&all_W, sizeof(double) * col) ); // prepare some values for cusolverDnZhegvd_bufferSize cusolverDnHandle_t cusolverH; - cusolverDnCreate(&cusolverH); + cusolverErrcheck( cusolverDnCreate(&cusolverH) ); int* devInfo; - cudaMalloc((void**)&devInfo, sizeof(int)); + checkCudaErrors( cudaMalloc((void**)&devInfo, sizeof(int)) ); // calculate the sizes needed for pre-allocated buffer. int lwork = 0; - cusolverDnZhegvd_bufferSize( + 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. @@ -52,14 +93,14 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, // (double2)B, row, all_W, - &lwork); + &lwork) ); // allocate memery cuDoubleComplex* d_work; - cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork); + checkCudaErrors( cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork) ); // compute eigenvalues and eigenvectors. - cusolverDnZhegvd( + 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. @@ -74,28 +115,29 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, all_W, d_work, lwork, - devInfo); + devInfo) ); - cudaDeviceSynchronize(); + checkCudaErrors( cudaDeviceSynchronize() ); // get eigenvalues and eigenvectors. only m ! - cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice); + checkCudaErrors( cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice) ); - cudaMemcpy(V, A_eigenvectors, sizeof(std::complex)*col*m, cudaMemcpyDeviceToDevice); + checkCudaErrors( cudaMemcpy(V, A_eigenvectors, sizeof(std::complex)*col*m, cudaMemcpyDeviceToDevice) ); matrixTranspose_op()(d, col, row, V, V); int info_gpu; - cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost); + checkCudaErrors( cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost) ); assert(0 == info_gpu); // free the buffer - cudaFree(d_work); + checkCudaErrors( cudaFree(d_work) ); // free resources and destroy - cudaFree(A_eigenvectors); - cudaFree(all_W); - cudaFree(devInfo); - cusolverDnDestroy(cusolverH); + checkCudaErrors( cudaFree(A_eigenvectors) ); + checkCudaErrors( cudaFree(all_W) ); + checkCudaErrors( cudaFree(devInfo) ); + cusolverErrcheck( cusolverDnDestroy(cusolverH) ); + destoryBLAShandle(); } @@ -109,29 +151,29 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, std::complex* V) { createBLAShandle(); + // init A_eigenvectors & transpose_B double2 *A_eigenvectors, *transpose_B; - cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col); - cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col); + checkCudaErrors( cudaMalloc((void**)&A_eigenvectors, sizeof(double2) * row * col) ); + checkCudaErrors( cudaMalloc((void**)&transpose_B, sizeof(double2) * row * col) ); - // transpose A to A_eigenvectors + // transpose A, B to A_eigenvectors, transpose_B matrixTranspose_op()(d, row, col, A, (std::complex*)A_eigenvectors); - // transpose B to transpose_B matrixTranspose_op()(d, row, col, B, (std::complex*)transpose_B); // init all_W double* all_W; - cudaMalloc((void**)&all_W, sizeof(double) * row); + checkCudaErrors( cudaMalloc((void**)&all_W, sizeof(double) * row) ); // prepare some values for cusolverDnZhegvd_bufferSize cusolverDnHandle_t cusolverH; - cusolverDnCreate(&cusolverH); + cusolverErrcheck( cusolverDnCreate(&cusolverH) ); int* devInfo; - cudaMalloc((void**)&devInfo, sizeof(int)); + checkCudaErrors( cudaMalloc((void**)&devInfo, sizeof(int)) ); // calculate the sizes needed for pre-allocated buffer. int lwork = 0; - cusolverDnZhegvd_bufferSize( + 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. @@ -142,14 +184,14 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, transpose_B, col, all_W, - &lwork); + &lwork) ); // allocate memery cuDoubleComplex* d_work; - cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork); + checkCudaErrors( cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork) ); // compute eigenvalues and eigenvectors. - cusolverDnZhegvd( + 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. @@ -162,25 +204,26 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, all_W, d_work, lwork, - devInfo); + devInfo) ); - cudaDeviceSynchronize(); + checkCudaErrors( cudaDeviceSynchronize() ); // get all eigenvalues and eigenvectors. - cudaMemcpy(W, all_W, sizeof(double) * row, cudaMemcpyDeviceToDevice); + checkCudaErrors( cudaMemcpy(W, all_W, sizeof(double) * row, cudaMemcpyDeviceToDevice) ); matrixTranspose_op()(d, row, col, (std::complex*)A_eigenvectors, V); int info_gpu; - cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost); + checkCudaErrors( cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost) ); assert(0 == info_gpu); // free the buffer - cudaFree(d_work); + checkCudaErrors( cudaFree(d_work) ); // free resources and destroy - cudaFree(A_eigenvectors); - cudaFree(all_W); - cudaFree(devInfo); - cusolverDnDestroy(cusolverH); + checkCudaErrors( cudaFree(A_eigenvectors) ); + checkCudaErrors( cudaFree(all_W) ); + checkCudaErrors( cudaFree(devInfo) ); + cusolverErrcheck( cusolverDnDestroy(cusolverH) ); + destoryBLAShandle(); } diff --git a/source/module_hsolver/src/cuda/math_kernel.cu b/source/module_hsolver/src/cuda/math_kernel.cu index 762c98e023f..ee796d981b1 100644 --- a/source/module_hsolver/src/cuda/math_kernel.cu +++ b/source/module_hsolver/src/cuda/math_kernel.cu @@ -404,7 +404,7 @@ void matrixTranspose_op::operator()(const psi::DEVICE_G ZERO.x = ZERO.y = 0.0; // use 'geam' API todo transpose. - 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); + 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; From 4143715f579cd0040e5340dc269e0d120400eec2 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Wed, 16 Nov 2022 13:24:48 +0800 Subject: [PATCH 12/14] Format dngvd_op.cu --- source/module_hsolver/src/cuda/dngvd_op.cu | 141 +++++++++++---------- 1 file changed, 73 insertions(+), 68 deletions(-) diff --git a/source/module_hsolver/src/cuda/dngvd_op.cu b/source/module_hsolver/src/cuda/dngvd_op.cu index 6e5dbc479ca..e67cb4cef44 100644 --- a/source/module_hsolver/src/cuda/dngvd_op.cu +++ b/source/module_hsolver/src/cuda/dngvd_op.cu @@ -1,47 +1,53 @@ #include "module_hsolver/include/dngvd_op.h" - #include "src_pdiag/helper_cuda.h" #include - -#define cusolverErrcheck(res) { cusolverAssert((res), __FILE__, __LINE__); } +#define cusolverErrcheck(res) \ + { \ + cusolverAssert((res), __FILE__, __LINE__); \ + } // cuSOLVER API errors -static const char *_cusolverGetErrorEnum(cusolverStatus_t error) { - switch (error) { +static const char* _cusolverGetErrorEnum(cusolverStatus_t error) +{ + switch (error) + { case CUSOLVER_STATUS_SUCCESS: - return "CUSOLVER_STATUS_SUCCESS"; + return "CUSOLVER_STATUS_SUCCESS"; case CUSOLVER_STATUS_NOT_INITIALIZED: - return "CUSOLVER_STATUS_NOT_INITIALIZED"; + return "CUSOLVER_STATUS_NOT_INITIALIZED"; case CUSOLVER_STATUS_ALLOC_FAILED: - return "CUSOLVER_STATUS_ALLOC_FAILED"; + return "CUSOLVER_STATUS_ALLOC_FAILED"; case CUSOLVER_STATUS_INVALID_VALUE: - return "CUSOLVER_STATUS_INVALID_VALUE"; + return "CUSOLVER_STATUS_INVALID_VALUE"; case CUSOLVER_STATUS_ARCH_MISMATCH: - return "CUSOLVER_STATUS_ARCH_MISMATCH"; + return "CUSOLVER_STATUS_ARCH_MISMATCH"; case CUSOLVER_STATUS_MAPPING_ERROR: - return "CUSOLVER_STATUS_MAPPING_ERROR"; + return "CUSOLVER_STATUS_MAPPING_ERROR"; case CUSOLVER_STATUS_EXECUTION_FAILED: - return "CUSOLVER_STATUS_EXECUTION_FAILED"; + return "CUSOLVER_STATUS_EXECUTION_FAILED"; case CUSOLVER_STATUS_INTERNAL_ERROR: - return "CUSOLVER_STATUS_INTERNAL_ERROR"; + return "CUSOLVER_STATUS_INTERNAL_ERROR"; case CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED: - return "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED"; + return "CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED"; case CUSOLVER_STATUS_NOT_SUPPORTED: - return "CUSOLVER_STATUS_NOT_SUPPORTED "; + return "CUSOLVER_STATUS_NOT_SUPPORTED "; case CUSOLVER_STATUS_ZERO_PIVOT: - return "CUSOLVER_STATUS_ZERO_PIVOT"; + return "CUSOLVER_STATUS_ZERO_PIVOT"; case CUSOLVER_STATUS_INVALID_LICENSE: - return "CUSOLVER_STATUS_INVALID_LICENSE"; - } - return ""; + 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); +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); } } @@ -59,28 +65,27 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, 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) ); + 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) ); + checkCudaErrors(cudaMalloc((void**)&all_W, sizeof(double) * col)); // prepare some values for cusolverDnZhegvd_bufferSize cusolverDnHandle_t cusolverH; - cusolverErrcheck( cusolverDnCreate(&cusolverH) ); + cusolverErrcheck(cusolverDnCreate(&cusolverH)); int* devInfo; - checkCudaErrors( cudaMalloc((void**)&devInfo, sizeof(int)) ); + checkCudaErrors(cudaMalloc((void**)&devInfo, sizeof(int))); // calculate the sizes needed for pre-allocated buffer. int lwork = 0; - cusolverErrcheck( cusolverDnZhegvd_bufferSize( + 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. @@ -93,14 +98,14 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, // (double2)B, row, all_W, - &lwork) ); + &lwork)); // allocate memery cuDoubleComplex* d_work; - checkCudaErrors( cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork) ); + checkCudaErrors(cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork)); // compute eigenvalues and eigenvectors. - cusolverErrcheck( cusolverDnZhegvd( + 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. @@ -115,29 +120,29 @@ void dngvx_op::operator()(const psi::DEVICE_GPU* d, all_W, d_work, lwork, - devInfo) ); + devInfo)); - checkCudaErrors( cudaDeviceSynchronize() ); + checkCudaErrors(cudaDeviceSynchronize()); // get eigenvalues and eigenvectors. only m ! - checkCudaErrors( cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice) ); + checkCudaErrors(cudaMemcpy(W, all_W, sizeof(double) * m, cudaMemcpyDeviceToDevice)); - checkCudaErrors( cudaMemcpy(V, A_eigenvectors, sizeof(std::complex)*col*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) ); + checkCudaErrors(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost)); assert(0 == info_gpu); - + // free the buffer - checkCudaErrors( cudaFree(d_work) ); + checkCudaErrors(cudaFree(d_work)); // free resources and destroy - checkCudaErrors( cudaFree(A_eigenvectors) ); - checkCudaErrors( cudaFree(all_W) ); - checkCudaErrors( cudaFree(devInfo) ); - cusolverErrcheck( cusolverDnDestroy(cusolverH) ); - + checkCudaErrors(cudaFree(A_eigenvectors)); + checkCudaErrors(cudaFree(all_W)); + checkCudaErrors(cudaFree(devInfo)); + cusolverErrcheck(cusolverDnDestroy(cusolverH)); + destoryBLAShandle(); } @@ -151,29 +156,29 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, 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) ); - + 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) ); + checkCudaErrors(cudaMalloc((void**)&all_W, sizeof(double) * row)); // prepare some values for cusolverDnZhegvd_bufferSize cusolverDnHandle_t cusolverH; - cusolverErrcheck( cusolverDnCreate(&cusolverH) ); + cusolverErrcheck(cusolverDnCreate(&cusolverH)); int* devInfo; - checkCudaErrors( cudaMalloc((void**)&devInfo, sizeof(int)) ); + checkCudaErrors(cudaMalloc((void**)&devInfo, sizeof(int))); // calculate the sizes needed for pre-allocated buffer. int lwork = 0; - cusolverErrcheck( cusolverDnZhegvd_bufferSize( + 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. @@ -184,14 +189,14 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, transpose_B, col, all_W, - &lwork) ); + &lwork)); // allocate memery cuDoubleComplex* d_work; - checkCudaErrors( cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork) ); + checkCudaErrors(cudaMalloc((void**)&d_work, sizeof(cuDoubleComplex) * lwork)); // compute eigenvalues and eigenvectors. - cusolverErrcheck( cusolverDnZhegvd( + 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. @@ -204,26 +209,26 @@ void dngv_op::operator()(const psi::DEVICE_GPU* d, all_W, d_work, lwork, - devInfo) ); + devInfo)); - checkCudaErrors( cudaDeviceSynchronize() ); + checkCudaErrors(cudaDeviceSynchronize()); // get all eigenvalues and eigenvectors. - checkCudaErrors( cudaMemcpy(W, all_W, sizeof(double) * row, cudaMemcpyDeviceToDevice) ); + 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) ); + checkCudaErrors(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost)); assert(0 == info_gpu); // free the buffer - checkCudaErrors( cudaFree(d_work) ); + checkCudaErrors(cudaFree(d_work)); // free resources and destroy - checkCudaErrors( cudaFree(A_eigenvectors) ); - checkCudaErrors( cudaFree(all_W) ); - checkCudaErrors( cudaFree(devInfo) ); - cusolverErrcheck( cusolverDnDestroy(cusolverH) ); - + checkCudaErrors(cudaFree(A_eigenvectors)); + checkCudaErrors(cudaFree(all_W)); + checkCudaErrors(cudaFree(devInfo)); + cusolverErrcheck(cusolverDnDestroy(cusolverH)); + destoryBLAShandle(); } From 5c8ad5c8bac7b67f5a5a4b1bc9ce4361daa7592f Mon Sep 17 00:00:00 2001 From: haozhihan Date: Wed, 16 Nov 2022 14:09:05 +0800 Subject: [PATCH 13/14] Format math_kernel.cu --- source/module_hsolver/src/cuda/math_kernel.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_hsolver/src/cuda/math_kernel.cu b/source/module_hsolver/src/cuda/math_kernel.cu index ee796d981b1..bb39ac62428 100644 --- a/source/module_hsolver/src/cuda/math_kernel.cu +++ b/source/module_hsolver/src/cuda/math_kernel.cu @@ -404,7 +404,7 @@ void matrixTranspose_op::operator()(const psi::DEVICE_G 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) ); + 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; From 9f54b2915b75ab55eccb8760f22bc553d94dbd61 Mon Sep 17 00:00:00 2001 From: haozhihan Date: Wed, 16 Nov 2022 14:30:26 +0800 Subject: [PATCH 14/14] Format some files --- source/module_hsolver/include/dngvd_op.h | 66 +- source/module_hsolver/src/dngvd_op.cpp | 126 +- .../module_hsolver/test/math_dngvd_test.cpp | 499 ++++---- .../module_hsolver/test/math_kernel_test.cpp | 1039 +++++++++-------- 4 files changed, 938 insertions(+), 792 deletions(-) diff --git a/source/module_hsolver/include/dngvd_op.h b/source/module_hsolver/include/dngvd_op.h index 9ac15c27ed2..2eb6643388b 100644 --- a/source/module_hsolver/include/dngvd_op.h +++ b/source/module_hsolver/include/dngvd_op.h @@ -4,70 +4,60 @@ #define MODULE_HSOLVER_DNGVD_H #include "module_base/lapack_connector.h" -#include "module_psi/include/memory.h" #include "module_hsolver/include/math_kernel.h" - +#include "module_psi/include/memory.h" namespace hsolver { - - -template -struct dngvx_op +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 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); + 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 +template struct dngv_op { - /// @brief DNGV computes all the eigenvalues and eigenvectors of a complex generalized + /// @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) + /// @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); + 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/src/dngvd_op.cpp b/source/module_hsolver/src/dngvd_op.cpp index f0e699cf773..62ef44c326b 100644 --- a/source/module_hsolver/src/dngvd_op.cpp +++ b/source/module_hsolver/src/dngvd_op.cpp @@ -6,15 +6,14 @@ 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) +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; @@ -31,49 +30,51 @@ void dngvx_op::operator()( 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 + } + else { lwork = (nb + 1) * col; } - std::complex *work = new std::complex[2 * lwork]; + std::complex* work = new std::complex[2 * lwork]; assert(work != 0); - double *rwork = new double[7 * col]; + double* rwork = new double[7 * col]; assert(rwork != 0); - int *iwork = new int[5 * col]; + int* iwork = new int[5 * col]; assert(iwork != 0); - int *ifail = new int[col]; + 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); + 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); @@ -81,19 +82,16 @@ void dngvx_op::operator()( 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) +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); @@ -110,7 +108,7 @@ void dngv_op::operator()( lwork = (nb + 1) * col; } - std::complex *work = new std::complex[lwork]; + std::complex* work = new std::complex[lwork]; ModuleBase::GlobalFunc::ZEROS(work, lwork); //===================================================================== @@ -120,7 +118,7 @@ void dngv_op::operator()( int info = 0; int rwork_dim; rwork_dim = 3 * col - 2; - double *rwork = new double[rwork_dim]; + double* rwork = new double[rwork_dim]; ModuleBase::GlobalFunc::ZEROS(rwork, rwork_dim); for (int i = 0; i < row * col; i++) @@ -143,7 +141,7 @@ void dngv_op::operator()( // const std::complex* A, // const std::complex* B, // const int m, -// float* W, +// float* W, // std::complex* V) // { // int lwork; @@ -183,27 +181,19 @@ void dngv_op::operator()( // '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 +// 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); +// 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; @@ -212,6 +202,4 @@ void dngv_op::operator()( // }; - - } // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/test/math_dngvd_test.cpp b/source/module_hsolver/test/math_dngvd_test.cpp index 17ab54afc3d..805987f23a9 100644 --- a/source/module_hsolver/test/math_dngvd_test.cpp +++ b/source/module_hsolver/test/math_dngvd_test.cpp @@ -1,38 +1,52 @@ -#include -#include -#include -#include -#include -#include -#include "module_psi/include/memory.h" #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: + 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_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_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 = {}; + 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} + 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}, @@ -44,35 +58,67 @@ class TestModuleHsolverMathDngvd : public ::testing::Test // {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}, + 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 + // 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> 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} - }; + 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) @@ -97,12 +143,12 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_cpu) { lwork = (nb + 1) * col; } - std::complex *work = new std::complex[lwork]; + 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]; + double* rwork = new double[rwork_dim]; ModuleBase::GlobalFunc::ZEROS(rwork, rwork_dim); // V.c = matrix_A.data(); // B.c = matrix_B.data(); @@ -119,22 +165,30 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_cpu) // (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() - ); - + 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; @@ -152,7 +206,7 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_cpu) // { // for (int j = 0; j < 4; j++) // { - // std::cout << V(i, j) << ", "; + // std::cout << V(i, j) << ", "; // } // std::cout << std::endl; // } @@ -171,9 +225,6 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_cpu) delete[] work; } - - - TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) { // (1) @@ -194,17 +245,18 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) 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 + } + else { lwork = (nb + 1) * col; } - std::complex *work = new std::complex[2 * lwork]; + std::complex* work = new std::complex[2 * lwork]; assert(work != 0); - double *rwork = new double[7 * col]; + double* rwork = new double[7 * col]; assert(rwork != 0); - int *iwork = new int[5 * col]; + int* iwork = new int[5 * col]; assert(iwork != 0); - int *ifail = new int[col]; + 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); @@ -223,53 +275,61 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) 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); - - + 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() - ); - + 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; @@ -287,7 +347,7 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) // { // for (int j = 0; j < 4; j++) // { - // std::cout << V(i, j) << ", "; + // std::cout << V(i, j) << ", "; // } // std::cout << std::endl; // } @@ -308,8 +368,6 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_cpu) delete[] ifail; } - - #if __UT_USE_CUDA || __UT_USE_ROCM // computes all the eigenvalues and eigenvectors @@ -322,51 +380,60 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) 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 + + // 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()); + 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 - ); + 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} - }; + // 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() - ); - + 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; + // 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; @@ -375,7 +442,7 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) // for (int j = 0; j < 4; j++) // { // std::cout << matrix_V_dngv_op[i * 4 + j] << ", "; - + // } // std::cout << std::endl; // } @@ -388,10 +455,9 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) // } // std::cout << std::endl; // } - - // we need to compare - // 1. W with W_result + // 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++) { @@ -399,10 +465,10 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) } 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); + 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); @@ -410,8 +476,6 @@ TEST_F(TestModuleHsolverMathDngvd, dngv_gpu) delete_memory_op_D()(gpu_ctx, device_W_dngv_op); } - - // computes the first m eigenvalues ​​and their corresponding eigenvectors TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) { @@ -422,48 +486,58 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) 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 + // 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()); - + 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 - ); + 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} - }; + // 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() - ); + 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; @@ -483,14 +557,13 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) // { // for (int j = 0; j < 4; j++) // { - // std::cout << matrix_V_DNGVX[i * 4 + j] << ", "; + // std::cout << matrix_V_DNGVX[i * 4 + j] << ", "; // } // std::cout << std::endl; // } - - // we need to compare - // 1. W with W_result + // we need to compare + // 1. W with W_result // 2. matrix_V with V_result for (int i = 0; i < W_DNGVX.size(); i++) { @@ -498,10 +571,10 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) } 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); + 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); @@ -509,28 +582,35 @@ TEST_F(TestModuleHsolverMathDngvd, dngvx_gpu) 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}, + 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}, + // 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()); @@ -540,10 +620,13 @@ TEST_F(TestModuleHsolverMathDngvd, transpose_gpu) // {-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}, + 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++) @@ -557,11 +640,9 @@ TEST_F(TestModuleHsolverMathDngvd, transpose_gpu) 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); - } + 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 31a8a335d5c..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,557 +83,631 @@ 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; - - // resize memory for values in GPU - resize_memory_op()(gpu_ctx, X_scal_dev, X_scal.size()); - + // in GPU + std::complex* X_scal_dev = NULL; - // syn the input data in CPU to GPU - synchronize_memory_op()(gpu_ctx, cpu_ctx, X_scal_dev, X_scal.data(), X_scal.size()); + // resize memory for values in GPU + resize_memory_op()(gpu_ctx, X_scal_dev, X_scal.size()); - // run - hsolver::createBLAShandle(); - scal_op_gpu()(gpu_ctx, dim, &alpha_scal, X_scal_dev, 1); - hsolver::destoryBLAShandle(); + // 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 output data in GPU to CPU - synchronize_memory_op_gpu()(cpu_ctx, gpu_ctx, X_scal.data(), X_scal_dev, X_scal.size()); + // 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()); - 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); - } + // 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); + delete_memory_op()(gpu_ctx, device_A); + delete_memory_op()(gpu_ctx, device_B); } - #endif // __UT_USE_CUDA || __UT_USE_ROCM /* @@ -628,7 +715,7 @@ TEST_F(TestModuleHsolverMathKernel, matrixSetToAnother_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); @@ -641,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