From 4666f5cc0985a3156c3282cfdb3120d420f1fdb6 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sun, 5 Nov 2023 00:01:09 +0800 Subject: [PATCH] real-valued CG --- source/Makefile.Objects | 1 + source/module_hsolver/CMakeLists.txt | 1 + source/module_hsolver/diagh.h | 9 +- source/module_hsolver/diagh_consts.cpp | 23 ++ source/module_hsolver/diago_cg.cpp | 71 +++-- source/module_hsolver/diago_cg.h | 1 + source/module_hsolver/diago_david.cpp | 23 -- source/module_hsolver/diago_david.h | 7 - .../kernels/cuda/math_kernel_op.cu | 80 +++++- .../module_hsolver/kernels/math_kernel_op.cpp | 45 +-- .../module_hsolver/kernels/math_kernel_op.h | 44 +-- .../kernels/rocm/math_kernel_op.hip.cu | 81 +++++- .../kernels/test/math_kernel_test.cpp | 8 +- source/module_hsolver/test/CMakeLists.txt | 18 +- .../test/diago_cg_real_test.cpp | 265 ++++++++++++++++++ .../test/diago_lcao_parallel_test.sh | 1 + 16 files changed, 552 insertions(+), 126 deletions(-) create mode 100644 source/module_hsolver/diagh_consts.cpp create mode 100644 source/module_hsolver/test/diago_cg_real_test.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 46bd1549568..82a19a8c182 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -272,6 +272,7 @@ OBJS_HCONTAINER=base_matrix.o\ OBJS_HSOLVER=diago_cg.o\ diago_david.o\ + diagh_consts.o\ diago_bpcg.o\ hsolver_pw.o\ hsolver_pw_sdft.o\ diff --git a/source/module_hsolver/CMakeLists.txt b/source/module_hsolver/CMakeLists.txt index 3b61324725f..6e3c70e869c 100644 --- a/source/module_hsolver/CMakeLists.txt +++ b/source/module_hsolver/CMakeLists.txt @@ -1,4 +1,5 @@ list(APPEND objects + diagh_consts.cpp diago_cg.cpp diago_david.cpp diago_bpcg.cpp diff --git a/source/module_hsolver/diagh.h b/source/module_hsolver/diagh.h index 56ce9f2771a..2d60c3e2327 100644 --- a/source/module_hsolver/diagh.h +++ b/source/module_hsolver/diagh.h @@ -3,11 +3,18 @@ #include -#include "module_base/complexmatrix.h" #include "module_base/macros.h" #include "module_hamilt_general/hamilt.h" #include "module_psi/psi.h" +template struct consts +{ + consts(); + T zero; + T one; + T neg_one; +}; + namespace hsolver { diff --git a/source/module_hsolver/diagh_consts.cpp b/source/module_hsolver/diagh_consts.cpp new file mode 100644 index 00000000000..f82f077a50c --- /dev/null +++ b/source/module_hsolver/diagh_consts.cpp @@ -0,0 +1,23 @@ +#include "diagh.h" +template class consts; +template class consts>; +template class consts>; + +template<> consts::consts() +{ + zero = 0.0; + one = 1.0; + neg_one = -1.0; +} +template<> consts>::consts() +{ + zero = std::complex(0.0, 0.0); + one = std::complex(1.0, 0.0); + neg_one = std::complex(-1.0, 0.0); +} +template<> consts>::consts() +{ + zero = std::complex(0.0, 0.0); + one = std::complex(1.0, 0.0); + neg_one = std::complex(-1.0, 0.0); +} diff --git a/source/module_hsolver/diago_cg.cpp b/source/module_hsolver/diago_cg.cpp index 7b0e45810d9..448771f5ee3 100644 --- a/source/module_hsolver/diago_cg.cpp +++ b/source/module_hsolver/diago_cg.cpp @@ -19,9 +19,9 @@ DiagoCG::DiagoCG(const Real* precondition_in) this->precondition = precondition_in; test_cg = 0; reorder = false; - this->one = new T(1.0, 0.0); - this->zero = new T(0.0, 0.0); - this->neg_one = new T(-1.0, 0.0); + this->one = &this->cs.one; + this->zero = &this->cs.zero; + this->neg_one = &this->cs.neg_one; } template @@ -35,9 +35,6 @@ DiagoCG::~DiagoCG() { delmem_complex_op()(this->ctx, this->gradient); delmem_complex_op()(this->ctx, this->g0); delmem_complex_op()(this->ctx, this->lagrange); - delete this->one; - delete this->zero; - delete this->neg_one; } template @@ -250,7 +247,7 @@ void DiagoCG::calculate_gradient() // this->gradient[i] -= lambda * this->pphi[i]; // } // haozhihan replace this 2022-10-6 - constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->gradient, this->gradient, 1.0, this->pphi, (-lambda)); + constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->gradient, this->gradient, 1.0, this->pphi, (-lambda)); } template @@ -313,6 +310,16 @@ void DiagoCG::orthogonal_gradient(hamilt::Hamilt *phm_in, 1); } +inline void init_real(const double in, double& out) +{ + out = in; +} +template +inline void init_real(const FPTYPE in, std::complex& out) +{ + out = std::complex(in, 0.0); +} + template void DiagoCG::calculate_gamma_cg(const int iter, Real &gg_last, const Real &cg_norm, const Real &theta) { @@ -377,10 +384,11 @@ void DiagoCG::calculate_gamma_cg(const int iter, Real &gg_last, const // pcg[i] = gamma * pcg[i] + this->gradient[i]; // } // haozhihan replace this 2022-10-6 - constantvector_addORsub_constantVector_op()(this->ctx, this->dim, pcg, pcg, gamma, this->gradient, 1.0); + constantvector_addORsub_constantVector_op()(this->ctx, this->dim, pcg, pcg, gamma, this->gradient, 1.0); const Real norma = gamma * cg_norm * sin(theta); - T znorma(norma * -1, 0.0); + T znorma; + init_real(norma * -1, znorma); // haozhihan replace this 2022-10-6 // const int one = 1; @@ -389,7 +397,7 @@ void DiagoCG::calculate_gamma_cg(const int iter, Real &gg_last, const { pcg[i] -= norma * pphi_m[i]; }*/ - axpy_op()(this->ctx, this->dim, &znorma, pphi_m, 1, pcg, 1); + axpy_op()(this->ctx, this->dim, &znorma, pphi_m, 1, pcg, 1); } } @@ -436,7 +444,7 @@ bool DiagoCG::update_psi(Real &cg_norm, Real &theta, Real &eigenvalue // } // haozhihan replace this 2022-10-6 - constantvector_addORsub_constantVector_op()(this->ctx, this->dim, phi_m_pointer, phi_m_pointer, cost, pcg, sint_norm); + constantvector_addORsub_constantVector_op()(this->ctx, this->dim, phi_m_pointer, phi_m_pointer, cost, pcg, sint_norm); // std::cout << "\n overlap2 = " << this->ddot(dim, phi_m, phi_m); @@ -455,12 +463,36 @@ bool DiagoCG::update_psi(Real &cg_norm, Real &theta, Real &eigenvalue // } // haozhihan replace this 2022-10-6 - constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->sphi, this->sphi, cost, this->scg, sint_norm); - constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->hphi, this->hphi, cost, this->pphi, sint_norm); + constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->sphi, this->sphi, cost, this->scg, sint_norm); + constantvector_addORsub_constantVector_op()(this->ctx, this->dim, this->hphi, this->hphi, cost, this->pphi, sint_norm); return 0; } } +inline double get_real(const double& x) +{ + return x; +} +inline double get_real(const std::complex& x) +{ + return x.real(); +} +inline float get_real(const std::complex& x) +{ + return x.real(); +} +inline double get_conj(const double& x) +{ + return x; +} +inline std::complex get_conj(const std::complex& x) +{ + return std::conj(x); +} +inline std::complex get_conj(const std::complex& x) +{ + return std::conj(x); +} template void DiagoCG::schmit_orth( const int &m, // end @@ -499,9 +531,9 @@ void DiagoCG::schmit_orth( // be careful , here reduce m+1 Parallel_Reduce::reduce_pool(lagrange_so, m + 1); - T var(0, 0); + T var = cs.zero; syncmem_complex_d2h_op()(this->cpu_ctx, this->ctx, &var, lagrange_so + m, 1); - Real psi_norm = var.real(); + Real psi_norm = get_real(var); //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< // haozhihan replace 2022-10-6 @@ -536,7 +568,7 @@ void DiagoCG::schmit_orth( std::cout << " m = " << m << std::endl; for (int j = 0; j <= m; ++j) { - std::cout << "j = " << j << " lagrange norm = " << (conj(lagrange_so[j]) * lagrange_so[j]).real() + std::cout << "j = " << j << " lagrange norm = " << get_real(get_conj(lagrange_so[j]) * lagrange_so[j]) << std::endl; } std::cout << " in DiagoCG, psi norm = " << psi_norm << std::endl; @@ -599,4 +631,11 @@ template class DiagoCG, psi::DEVICE_CPU>; template class DiagoCG, psi::DEVICE_GPU>; template class DiagoCG, psi::DEVICE_GPU>; #endif + +#ifdef __LCAO +template class DiagoCG; +#if ((defined __CUDA) || (defined __ROCM)) +template class DiagoCG; +#endif +#endif } // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/diago_cg.h b/source/module_hsolver/diago_cg.h index 36fe4b6fc43..13d0a87e08c 100644 --- a/source/module_hsolver/diago_cg.h +++ b/source/module_hsolver/diago_cg.h @@ -111,6 +111,7 @@ class DiagoCG : public DiagH using setmem_var_h_op = psi::memory::set_memory_op; using syncmem_var_h2d_op = psi::memory::synchronize_memory_op; + consts cs; const T * one = nullptr, * zero = nullptr, * neg_one = nullptr; }; diff --git a/source/module_hsolver/diago_david.cpp b/source/module_hsolver/diago_david.cpp index 820e8a7f81e..3986c39704f 100644 --- a/source/module_hsolver/diago_david.cpp +++ b/source/module_hsolver/diago_david.cpp @@ -17,29 +17,6 @@ using namespace hsolver; -template class consts; -template class consts>; -template class consts>; - -template<> consts::consts() -{ - zero = 0.0; - one = 1.0; - neg_one = -1.0; -} -template<> consts>::consts() -{ - zero = std::complex(0.0, 0.0); - one = std::complex(1.0, 0.0); - neg_one = std::complex(-1.0, 0.0); -} -template<> consts>::consts() -{ - zero = std::complex(0.0, 0.0); - one = std::complex(1.0, 0.0); - neg_one = std::complex(-1.0, 0.0); -} - inline double get_real(const double& x) { return x; diff --git a/source/module_hsolver/diago_david.h b/source/module_hsolver/diago_david.h index b3edcde9e16..0940932769c 100644 --- a/source/module_hsolver/diago_david.h +++ b/source/module_hsolver/diago_david.h @@ -16,13 +16,6 @@ #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include "module_psi/kernels/device.h" -template struct consts -{ - consts(); - T zero; - T one; - T neg_one; -}; namespace hsolver { diff --git a/source/module_hsolver/kernels/cuda/math_kernel_op.cu b/source/module_hsolver/kernels/cuda/math_kernel_op.cu index c52af8a5994..5a0f74eb1d1 100644 --- a/source/module_hsolver/kernels/cuda/math_kernel_op.cu +++ b/source/module_hsolver/kernels/cuda/math_kernel_op.cu @@ -260,14 +260,14 @@ __global__ void vector_div_vector_kernel( } } -template +template __global__ void constantvector_addORsub_constantVector_kernel( const int size, - thrust::complex* result, - const thrust::complex* vector1, - const FPTYPE constant1, - const thrust::complex* vector2, - const FPTYPE constant2) + T* result, + const T* vector1, + const typename GetTypeReal::type constant1, + const T* vector2, + const typename GetTypeReal::type constant2) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i < size) @@ -563,8 +563,22 @@ void vector_div_vector_op, psi::DEVICE_GPU>::operator()( vector_div_vector_complex_wrapper(d, dim, result, vector1, vector2); } // vector operator: result[i] = vector1[i] * constant1 + vector2[i] * constant2 -template -void constantvector_addORsub_constantVector_op::operator()( +template <> +void constantvector_addORsub_constantVector_op::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + double* result, + const double* vector1, + const double constant1, + const double* vector2, + const double constant2) +{ + int thread = 1024; + int block = (dim + thread - 1) / thread; + constantvector_addORsub_constantVector_kernel << > > (dim, result, vector1, constant1, vector2, constant2); +} +template +inline void constantvector_addORsub_constantVector_complex_wrapper( const psi::DEVICE_GPU* d, const int& dim, std::complex* result, @@ -579,11 +593,48 @@ void constantvector_addORsub_constantVector_op::operato int thread = 1024; int block = (dim + thread - 1) / thread; - constantvector_addORsub_constantVector_kernel<<>>(dim, result_tmp, vector1_tmp,constant1, vector2_tmp, constant2); + constantvector_addORsub_constantVector_kernel> << > > (dim, result_tmp, vector1_tmp, constant1, vector2_tmp, constant2); +} +template <> +void constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + std::complex* result, + const std::complex* vector1, + const float constant1, + const std::complex* vector2, + const float constant2) +{ + constantvector_addORsub_constantVector_complex_wrapper(d, dim, result, vector1, constant1, vector2, constant2); +} +template <> +void constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + std::complex* result, + const std::complex* vector1, + const double constant1, + const std::complex* vector2, + const double constant2) +{ + constantvector_addORsub_constantVector_complex_wrapper(d, dim, result, vector1, constant1, vector2, constant2); } -template <> -void axpy_op::operator()( +template <> +void axpy_op::operator()( + const psi::DEVICE_GPU* d, + const int& N, + const double* alpha, + const double* X, + const int& incX, + double* Y, + const int& incY) +{ + cublasErrcheck(cublasDaxpy(cublas_handle, N, alpha, X, incX, Y, incY)); +} + +template <> +void axpy_op, psi::DEVICE_GPU>::operator()( const psi::DEVICE_GPU* d, const int& N, const std::complex *alpha, @@ -596,7 +647,7 @@ void axpy_op::operator()( } template <> -void axpy_op::operator()( +void axpy_op, psi::DEVICE_GPU>::operator()( const psi::DEVICE_GPU* d, const int& N, const std::complex *alpha, @@ -969,7 +1020,7 @@ template struct line_minimize_with_block_op, psi::DEVICE_GPU template struct vector_div_constant_op, psi::DEVICE_GPU>; template struct vector_mul_vector_op, psi::DEVICE_GPU>; template struct vector_div_vector_op, psi::DEVICE_GPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>; template struct matrixSetToAnother, psi::DEVICE_GPU>; template struct dot_real_op, psi::DEVICE_GPU>; @@ -978,7 +1029,7 @@ template struct line_minimize_with_block_op, psi::DEVICE_GP template struct vector_div_constant_op, psi::DEVICE_GPU>; template struct vector_mul_vector_op, psi::DEVICE_GPU>; template struct vector_div_vector_op, psi::DEVICE_GPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>; template struct matrixSetToAnother, psi::DEVICE_GPU>; #ifdef __LCAO @@ -987,5 +1038,6 @@ template struct vector_div_constant_op; template struct vector_mul_vector_op; template struct vector_div_vector_op; template struct matrixSetToAnother; +template struct constantvector_addORsub_constantVector_op; #endif } // namespace hsolver diff --git a/source/module_hsolver/kernels/math_kernel_op.cpp b/source/module_hsolver/kernels/math_kernel_op.cpp index 2d3d44c78cb..597d9cd5eae 100644 --- a/source/module_hsolver/kernels/math_kernel_op.cpp +++ b/source/module_hsolver/kernels/math_kernel_op.cpp @@ -189,18 +189,19 @@ template struct vector_div_vector_op } }; -template struct constantvector_addORsub_constantVector_op +template struct constantvector_addORsub_constantVector_op { + using Real = typename GetTypeReal::type; void operator()(const psi::DEVICE_CPU* d, - const int& dim, - std::complex* result, - const std::complex* vector1, - const FPTYPE constant1, - const std::complex* vector2, - const FPTYPE constant2) + const int& dim, + T* result, + const T* vector1, + const Real constant1, + const T* vector2, + const Real constant2) { #ifdef _OPENMP -#pragma omp parallel for schedule(static, 4096/sizeof(FPTYPE)) +#pragma omp parallel for schedule(static, 8192/sizeof(T)) #endif for (int i = 0; i < dim; i++) { @@ -242,16 +243,16 @@ struct gemv_op { } }; -template -struct axpy_op { +template +struct axpy_op { void operator()( - const psi::DEVICE_CPU * /*ctx*/, - const int &dim, - const std::complex *alpha, - const std::complex *X, - const int &incX, - std::complex *Y, - const int &incY) + const psi::DEVICE_CPU* /*ctx*/, + const int& dim, + const T* alpha, + const T* X, + const int& incX, + T* Y, + const int& incY) { BlasConnector::axpy(dim, *alpha, X, incX, Y, incY); } @@ -336,34 +337,35 @@ template struct matrixSetToAnother // Explicitly instantiate functors for the types of functor registered. template struct scal_op; -template struct axpy_op; +template struct axpy_op, psi::DEVICE_CPU>; template struct gemv_op, psi::DEVICE_CPU>; template struct gemm_op, psi::DEVICE_CPU>; template struct dot_real_op, psi::DEVICE_CPU>; template struct vector_div_constant_op, psi::DEVICE_CPU>; template struct vector_mul_vector_op, psi::DEVICE_CPU>; template struct vector_div_vector_op, psi::DEVICE_CPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_CPU>; template struct matrixTranspose_op, psi::DEVICE_CPU>; template struct matrixSetToAnother, psi::DEVICE_CPU>; template struct calc_grad_with_block_op, psi::DEVICE_CPU>; template struct line_minimize_with_block_op, psi::DEVICE_CPU>; template struct scal_op; -template struct axpy_op; +template struct axpy_op, psi::DEVICE_CPU>; template struct gemv_op, psi::DEVICE_CPU>; template struct gemm_op, psi::DEVICE_CPU>; template struct dot_real_op, psi::DEVICE_CPU>; template struct vector_div_constant_op, psi::DEVICE_CPU>; template struct vector_mul_vector_op, psi::DEVICE_CPU>; template struct vector_div_vector_op, psi::DEVICE_CPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_CPU>; template struct matrixTranspose_op, psi::DEVICE_CPU>; template struct matrixSetToAnother, psi::DEVICE_CPU>; template struct calc_grad_with_block_op, psi::DEVICE_CPU>; template struct line_minimize_with_block_op, psi::DEVICE_CPU>; #ifdef __LCAO +template struct axpy_op; template struct gemv_op; template struct gemm_op; template struct dot_real_op; @@ -372,5 +374,6 @@ template struct vector_div_constant_op; template struct vector_div_vector_op; template struct matrixTranspose_op; template struct matrixSetToAnother; +template struct constantvector_addORsub_constantVector_op; #endif } // namespace hsolver \ No newline at end of file diff --git a/source/module_hsolver/kernels/math_kernel_op.h b/source/module_hsolver/kernels/math_kernel_op.h index 1bff8dfa858..55f39fde340 100644 --- a/source/module_hsolver/kernels/math_kernel_op.h +++ b/source/module_hsolver/kernels/math_kernel_op.h @@ -217,8 +217,9 @@ template struct vector_div_vector_op }; // vector operator: result[i] = vector1[i] * constant1 + vector2[i] * constant2 -template struct constantvector_addORsub_constantVector_op +template struct constantvector_addORsub_constantVector_op { + using Real = typename GetTypeReal::type; /// @brief result[i] = vector1[i] * constant1 + vector2[i] * constant2 /// /// Input Parameters @@ -232,16 +233,16 @@ template struct constantvector_addORsub_const /// Output Parameters /// \param result : output array void operator()(const Device* d, - const int& dim, - std::complex* result, - const std::complex* vector1, - const FPTYPE constant1, - const std::complex* vector2, - const FPTYPE constant2); + const int& dim, + T* result, + const T* vector1, + const Real constant1, + const T* vector2, + const Real constant2); }; // compute Y = alpha * X + Y -template struct axpy_op +template struct axpy_op { /// @brief Y = alpha * X + Y /// @@ -257,12 +258,12 @@ template struct axpy_op /// Output Parameters /// \param Y : output array Y void operator()(const Device* d, - const int& N, - const std::complex* alpha, - const std::complex* X, - const int& incX, - std::complex* Y, - const int& incY); + const int& N, + const T* alpha, + const T* X, + const int& incX, + T* Y, + const int& incY); }; // compute y = alpha * op(A) * x + beta * y @@ -457,15 +458,16 @@ template struct vector_div_vector_op }; // vector operator: result[i] = vector1[i] * constant1 + vector2[i] * constant2 -template struct constantvector_addORsub_constantVector_op +template struct constantvector_addORsub_constantVector_op { + using Real = typename GetTypeReal::type; void operator()(const psi::DEVICE_GPU* d, - const int& dim, - std::complex* result, - const std::complex* vector1, - const FPTYPE constant1, - const std::complex* vector2, - const FPTYPE constant2); + const int& dim, + T* result, + const T* vector1, + const Real constant1, + const T* vector2, + const Real constant2); }; template struct matrixSetToAnother diff --git a/source/module_hsolver/kernels/rocm/math_kernel_op.hip.cu b/source/module_hsolver/kernels/rocm/math_kernel_op.hip.cu index 8ea97f3742f..a03ac0bc7ac 100644 --- a/source/module_hsolver/kernels/rocm/math_kernel_op.hip.cu +++ b/source/module_hsolver/kernels/rocm/math_kernel_op.hip.cu @@ -248,15 +248,15 @@ __global__ void vector_div_vector_kernel( } } -template +template __launch_bounds__(1024) __global__ void constantvector_addORsub_constantVector_kernel( const int size, - thrust::complex* result, - const thrust::complex* vector1, - const FPTYPE constant1, - const thrust::complex* vector2, - const FPTYPE constant2) + T* result, + const T* vector1, + const typename GetTypeReal::type constant1, + const T* vector2, + const typename GetTypeReal::type constant2) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i < size) @@ -553,8 +553,23 @@ void vector_div_vector_op, psi::DEVICE_GPU>::operator()( } // vector operator: result[i] = vector1[i] * constant1 + vector2[i] * constant2 -template -void constantvector_addORsub_constantVector_op::operator()( +template <> +void constantvector_addORsub_constantVector_op::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + double* result, + const double* vector1, + const double constant1, + const double* vector2, + const double constant2) +{ + int thread = 1024; + int block = (dim + thread - 1) / thread; + hipLaunchKernelGGL(HIP_KERNEL_NAME(constantvector_addORsub_constantVector_kernel), dim3(block), dim3(thread), 0, 0, dim, result, vector1, constant1, vector2, constant2); +} + +template +inline void constantvector_addORsub_constantVector_complex_wrapper( const psi::DEVICE_GPU* d, const int& dim, std::complex* result, @@ -569,11 +584,48 @@ void constantvector_addORsub_constantVector_op::operato int thread = 1024; int block = (dim + thread - 1) / thread; - hipLaunchKernelGGL(HIP_KERNEL_NAME(constantvector_addORsub_constantVector_kernel), dim3(block), dim3(thread), 0, 0, dim, result_tmp, vector1_tmp,constant1, vector2_tmp, constant2); + hipLaunchKernelGGL(HIP_KERNEL_NAME(constantvector_addORsub_constantVector_kernel>), dim3(block), dim3(thread), 0, 0, dim, result_tmp, vector1_tmp, constant1, vector2_tmp, constant2); +} +template <> +void constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + std::complex* result, + const std::complex* vector1, + const float constant1, + const std::complex* vector2, + const float constant2) +{ + constantvector_addORsub_constantVector_complex_wrapper(d, dim, result, vector1, constant1, vector2, constant2); +} +template <> +void constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>::operator()( + const psi::DEVICE_GPU* d, + const int& dim, + std::complex* result, + const std::complex* vector1, + const double constant1, + const std::complex* vector2, + const double constant2) +{ + constantvector_addORsub_constantVector_complex_wrapper(d, dim, result, vector1, constant1, vector2, constant2); } -template <> -void axpy_op::operator()( +template <> +void axpy_op::operator()( + const psi::DEVICE_GPU* d, + const int& N, + const double* alpha, + const double* X, + const int& incX, + double* Y, + const int& incY) +{ + hipblasDaxpy(cublas_handle, N, alpha, X, incX, Y, incY); +} + +template <> +void axpy_op, psi::DEVICE_GPU>::operator()( const psi::DEVICE_GPU* d, const int& N, const std::complex *alpha, @@ -586,7 +638,7 @@ void axpy_op::operator()( } template <> -void axpy_op::operator()( +void axpy_op, psi::DEVICE_GPU>::operator()( const psi::DEVICE_GPU* d, const int& N, const std::complex *alpha, @@ -956,7 +1008,7 @@ template struct line_minimize_with_block_op, psi::DEVICE_GPU template struct vector_div_constant_op, psi::DEVICE_GPU>; template struct vector_mul_vector_op, psi::DEVICE_GPU>; template struct vector_div_vector_op, psi::DEVICE_GPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>; template struct matrixSetToAnother, psi::DEVICE_GPU>; template struct dot_real_op, psi::DEVICE_GPU>; @@ -965,7 +1017,7 @@ template struct line_minimize_with_block_op, psi::DEVICE_GP template struct vector_div_constant_op, psi::DEVICE_GPU>; template struct vector_mul_vector_op, psi::DEVICE_GPU>; template struct vector_div_vector_op, psi::DEVICE_GPU>; -template struct constantvector_addORsub_constantVector_op; +template struct constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>; template struct matrixSetToAnother, psi::DEVICE_GPU>; #ifdef __LCAO @@ -974,5 +1026,6 @@ template struct vector_div_constant_op; template struct vector_mul_vector_op; template struct vector_div_vector_op; template struct matrixSetToAnother; +template struct constantvector_addORsub_constantVector_op; #endif } // namespace hsolver diff --git a/source/module_hsolver/kernels/test/math_kernel_test.cpp b/source/module_hsolver/kernels/test/math_kernel_test.cpp index 1a1c6905e82..6f578852ecb 100644 --- a/source/module_hsolver/kernels/test/math_kernel_test.cpp +++ b/source/module_hsolver/kernels/test/math_kernel_test.cpp @@ -75,8 +75,8 @@ class TestModuleHsolverMathKernel : public ::testing::Test using vector_mul_vector_op_cpu = hsolver::vector_mul_vector_op, psi::DEVICE_CPU>; using vector_div_vector_op_cpu = hsolver::vector_div_vector_op, psi::DEVICE_CPU>; using constantvector_addORsub_constantVector_op_cpu - = hsolver::constantvector_addORsub_constantVector_op; - using axpy_op_cpu = hsolver::axpy_op; + = hsolver::constantvector_addORsub_constantVector_op, psi::DEVICE_CPU>; + using axpy_op_cpu = hsolver::axpy_op, psi::DEVICE_CPU>; using scal_op_cpu = hsolver::scal_op; using gemv_op_cpu = hsolver::gemv_op, psi::DEVICE_CPU>; // gpu operator @@ -84,8 +84,8 @@ class TestModuleHsolverMathKernel : public ::testing::Test using vector_mul_vector_op_gpu = hsolver::vector_mul_vector_op, psi::DEVICE_GPU>; using vector_div_vector_op_gpu = hsolver::vector_div_vector_op, psi::DEVICE_GPU>; using constantvector_addORsub_constantVector_op_gpu - = hsolver::constantvector_addORsub_constantVector_op; - using axpy_op_gpu = hsolver::axpy_op; + = hsolver::constantvector_addORsub_constantVector_op, psi::DEVICE_GPU>; + using axpy_op_gpu = hsolver::axpy_op, psi::DEVICE_GPU>; using scal_op_gpu = hsolver::scal_op; using gemv_op_gpu = hsolver::gemv_op, psi::DEVICE_GPU>; diff --git a/source/module_hsolver/test/CMakeLists.txt b/source/module_hsolver/test/CMakeLists.txt index eef80715c17..00312d1eadd 100644 --- a/source/module_hsolver/test/CMakeLists.txt +++ b/source/module_hsolver/test/CMakeLists.txt @@ -13,7 +13,7 @@ AddTest( AddTest( TARGET HSolver_cg LIBS ${math_libs} base psi device - SOURCES diago_cg_test.cpp ../diago_cg.cpp ../diago_iter_assist.cpp + SOURCES diago_cg_test.cpp ../diago_cg.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp ../../module_basis/module_pw/test/test_tool.cpp ../../module_hamilt_general/operator.cpp ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp @@ -21,7 +21,7 @@ AddTest( AddTest( TARGET HSolver_cg_float LIBS ${math_libs} base psi device - SOURCES diago_cg_float_test.cpp ../diago_cg.cpp ../diago_iter_assist.cpp + SOURCES diago_cg_float_test.cpp ../diago_cg.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp ../../module_basis/module_pw/test/test_tool.cpp ../../module_hamilt_general/operator.cpp ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp @@ -29,7 +29,7 @@ AddTest( AddTest( TARGET HSolver_dav LIBS ${math_libs} base psi device - SOURCES diago_david_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp + SOURCES diago_david_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp ../../module_basis/module_pw/test/test_tool.cpp ../../module_hamilt_general/operator.cpp ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp @@ -37,16 +37,24 @@ AddTest( AddTest( TARGET HSolver_dav_float LIBS ${math_libs} base psi device - SOURCES diago_david_float_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp + SOURCES diago_david_float_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp ../../module_basis/module_pw/test/test_tool.cpp ../../module_hamilt_general/operator.cpp ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp ) if(ENABLE_LCAO) +AddTest( + TARGET HSolver_cg_real + LIBS ${math_libs} base psi device + SOURCES diago_cg_float_test.cpp ../diago_cg.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp + ../../module_basis/module_pw/test/test_tool.cpp + ../../module_hamilt_general/operator.cpp + ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp +) AddTest( TARGET HSolver_dav_real LIBS ${math_libs} base psi device - SOURCES diago_david_real_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp + SOURCES diago_david_real_test.cpp ../diago_david.cpp ../diago_iter_assist.cpp ../diagh_consts.cpp ../../module_basis/module_pw/test/test_tool.cpp ../../module_hamilt_general/operator.cpp ../../module_hamilt_pw/hamilt_pwdft/operator_pw/operator_pw.cpp diff --git a/source/module_hsolver/test/diago_cg_real_test.cpp b/source/module_hsolver/test/diago_cg_real_test.cpp new file mode 100644 index 00000000000..013fcaf47f8 --- /dev/null +++ b/source/module_hsolver/test/diago_cg_real_test.cpp @@ -0,0 +1,265 @@ +#include "module_base/inverse_matrix.h" +#include "module_base/lapack_connector.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" +#include "module_psi/psi.h" +#include "module_hamilt_general/hamilt.h" +#include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h" +#include "../diago_cg.h" +#include "../diago_iter_assist.h" +#include "diago_mock.h" +#include "mpi.h" +#include "module_basis/module_pw/test/test_tool.h" +#include + +#include "gtest/gtest.h" +#include + +/************************************************ + * unit test of functions in Diago_CG + ***********************************************/ + + /** + * Class Diago_CG is an approach for eigenvalue problems + * This unittest test the function Diago_CG::diag() for FPTYPE=double, Device=cpu + * with different examples. + * - the Hermite matrices (npw=500,1000) produced using random numbers and with sparsity of 0%, 60%, 80% + * - the Hamiltonian matrix read from "data-H", produced by using out_hs in INPUT of a LCAO calculation + * - a 2x2 Hermite matrix for learning and checking + * + * Note: + * The test is passed when the eignvalues are closed to these calculated by LAPACK. + * It is used together with a header file diago_mock.h. + * The default Hermite matrix generated here is real symmetric, one can add an imaginary part + * by changing two commented out lines in diago_mock.h. + * + */ + + // call lapack in order to compare to cg +void lapackEigen(int& npw, std::vector& hm, double* e, bool outtime = false) +{ + int info = 0; + auto tmp = hm; + clock_t start, end; + start = clock(); + char tmp_c1 = 'V', tmp_c2 = 'U'; + + double work_tmp; + constexpr int minus_one = -1; + dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, &work_tmp, &minus_one, &info); // get best lwork + + const int lwork = work_tmp; + double* work2 = new double[lwork]; + dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, &info); + end = clock(); + if (info) std::cout << "ERROR: Lapack solver, info=" << info << std::endl; + if (outtime) std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl; + delete[] work2; +} + +class DiagoCGPrepare +{ +public: + DiagoCGPrepare(int nband, int npw, int sparsity, bool reorder, double eps, int maxiter, double threshold) + : nband(nband), npw(npw), sparsity(sparsity), reorder(reorder), eps(eps), maxiter(maxiter), + threshold(threshold) + { +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &mypnum); +#endif + } + + int nband, npw, sparsity, maxiter, notconv; + // eps is the convergence threshold within cg_diago + double eps, avg_iter; + bool reorder; + double threshold; + int nprocs = 1, mypnum = 0; + // threshold is the comparison standard between cg and lapack + + void CompareEigen(double* precondition) + { + // calculate eigenvalues by LAPACK; + double* e_lapack = new double[npw]; + auto ev = DIAGOTEST::hmatrix_d; + if (mypnum == 0) lapackEigen(npw, ev, e_lapack, false); + // initial guess of psi by perturbing lapack psi + ModuleBase::matrix psiguess(nband, npw); + std::default_random_engine p(1); + std::uniform_int_distribution u(1, 10); + for (int i = 0; i < nband; i++) + { + for (int j = 0; j < npw; j++) + { + double rand = static_cast(u(p)) / 10.; + // psiguess(i,j) = ev(j,i)*(1+rand); + psiguess(i, j) = ev[j * DIAGOTEST::h_nc + i] * rand; + } + } + // run cg + //====================================================================== + double* en = new double[npw]; + int ik = 1; + hamilt::Hamilt* ha; + ha = new hamilt::HamiltPW(nullptr, nullptr, nullptr); + int* ngk = new int[1]; + psi::Psi psi; + psi.resize(ik, nband, npw); + //psi.fix_k(0); + for (int i = 0; i < nband; i++) + { + for (int j = 0; j < npw; j++) + { + psi(i, j) = psiguess(i, j); + } + } + + psi::Psi psi_local; + double* precondition_local; + DIAGOTEST::npw_local = new int[nprocs]; +#ifdef __MPI + DIAGOTEST::cal_division(DIAGOTEST::npw); + DIAGOTEST::divide_hpsi(psi, psi_local, DIAGOTEST::hmatrix_d, DIAGOTEST::hmatrix_local_d); //will distribute psi and Hmatrix to each process + precondition_local = new double[DIAGOTEST::npw_local[mypnum]]; + DIAGOTEST::divide_psi(precondition, precondition_local); +#else + DIAGOTEST::hmatrix_local_d = DIAGOTEST::hmatrix_d; + DIAGOTEST::npw_local[0] = DIAGOTEST::npw; + psi_local = psi; + precondition_local = new double[DIAGOTEST::npw]; + for (int i = 0;i < DIAGOTEST::npw;i++) precondition_local[i] = precondition[i]; +#endif + hsolver::DiagoCG cg(precondition_local); + psi_local.fix_k(0); + double start, end; + start = MPI_Wtime(); + cg.diag(ha, psi_local, en); + end = MPI_Wtime(); + //if(mypnum == 0) printf("diago time:%7.3f\n",end-start); + delete[] DIAGOTEST::npw_local; + delete[] precondition_local; + //====================================================================== + for (int i = 0; i < nband; i++) + { + EXPECT_NEAR(en[i], e_lapack[i], threshold); + } + + delete[] en; + delete[] e_lapack; + delete ha; + } +}; + +class DiagoCGTest : public ::testing::TestWithParam +{ +}; + +TEST_P(DiagoCGTest, RandomHamilt) +{ + DiagoCGPrepare dcp = GetParam(); + hsolver::DiagoIterAssist::PW_DIAG_NMAX = dcp.maxiter; + hsolver::DiagoIterAssist::PW_DIAG_THR = dcp.eps; + HPsi hpsi(dcp.nband, dcp.npw, dcp.sparsity); + DIAGOTEST::hmatrix_d = hpsi.hamilt(); + + DIAGOTEST::npw = dcp.npw; + // ModuleBase::matrix psi = hpsi.psi(); + dcp.CompareEigen(hpsi.precond()); +} + +INSTANTIATE_TEST_SUITE_P(VerifyCG, + DiagoCGTest, + ::testing::Values( + // nband, npw, sparsity, reorder, eps, maxiter, threshold + DiagoCGPrepare(10, 500, 0, true, 1e-5, 300, 1e-3), + DiagoCGPrepare(20, 500, 6, true, 1e-5, 300, 1e-3), + DiagoCGPrepare(20, 1000, 8, true, 1e-5, 300, 1e-3), + DiagoCGPrepare(40, 1000, 8, true, 1e-6, 300, 1e-3))); +//DiagoCGPrepare(40, 2000, 8, true, 1e-5, 500, 1e-2))); +// the last one is passed but time-consumming. + +// check that the mock class HPsi work well +// in generating a Hermite matrix +TEST(DiagoCGTest, Hamilt) +{ + int dim = 2; + int nbnd = 2; + HPsi hpsi(nbnd, dim); + std::vector hm = hpsi.hamilt(); + EXPECT_EQ(DIAGOTEST::h_nr, 2); + EXPECT_EQ(DIAGOTEST::h_nc, 2); + EXPECT_EQ(hm[DIAGOTEST::h_nc], hm[1]); +} + +// cg for a 2x2 matrix +#ifdef __MPI +#else +TEST(DiagoCGTest, TwoByTwo) +{ + int dim = 2; + int nband = 2; + ModuleBase::matrix hm(2, 2); + hm(0, 0) = 4.0; + hm(0, 1) = 1.0; + hm(1, 0) = 1.0; + hm(1, 1) = 3.0; + // nband, npw, sub, sparsity, reorder, eps, maxiter, threshold + DiagoCGPrepare dcp(nband, dim, 0, true, 1e-4, 50, 1e-10); + hsolver::DiagoIterAssist::PW_DIAG_NMAX = dcp.maxiter; + hsolver::DiagoIterAssist::PW_DIAG_THR = dcp.eps; + HPsi hpsi; + hpsi.create(nband, dim); + DIAGOTEST::hmatrix_d = hm; + DIAGOTEST::npw = dim; + dcp.CompareEigen(hpsi.precond()); +} +#endif + +TEST(DiagoCGTest, readH) +{ + // read Hamilt matrix from file data-H + std::vector hm; + std::ifstream ifs; + ifs.open("H-GammaOnly-Si64.dat"); + DIAGOTEST::readh(ifs, hm); + ifs.close(); + int dim = DIAGOTEST::npw; + int nband = 10; // not nband < dim, here dim = 26 in data-H + // nband, npw, sub, sparsity, reorder, eps, maxiter, threshold + DiagoCGPrepare dcp(nband, dim, 0, true, 1e-5, 500, 1e-3); + hsolver::DiagoIterAssist::PW_DIAG_NMAX = dcp.maxiter; + hsolver::DiagoIterAssist::PW_DIAG_THR = dcp.eps; + HPsi hpsi; + hpsi.create(nband, dim); + DIAGOTEST::hmatrix_d = hpsi.hamilt(); + DIAGOTEST::npw = dim; + dcp.CompareEigen(hpsi.precond()); +} + +int main(int argc, char** argv) +{ + int nproc = 1, myrank = 0; + +#ifdef __MPI + int nproc_in_pool, kpar = 1, mypool, rank_in_pool; + setupmpi(argc, argv, nproc, myrank); + divide_pools(nproc, myrank, nproc_in_pool, kpar, mypool, rank_in_pool); + GlobalV::NPROC_IN_POOL = nproc; +#else + MPI_Init(&argc, &argv); +#endif + + testing::InitGoogleTest(&argc, argv); + ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners(); + if (myrank != 0) delete listeners.Release(listeners.default_result_printer()); + + int result = RUN_ALL_TESTS(); + if (myrank == 0 && result != 0) + { + std::cout << "ERROR:some tests are not passed" << std::endl; + return result; + } + + MPI_Finalize(); + return 0; +} diff --git a/source/module_hsolver/test/diago_lcao_parallel_test.sh b/source/module_hsolver/test/diago_lcao_parallel_test.sh index 276fc75ce01..30cd96f43a0 100644 --- a/source/module_hsolver/test/diago_lcao_parallel_test.sh +++ b/source/module_hsolver/test/diago_lcao_parallel_test.sh @@ -10,6 +10,7 @@ for i in 6 3 2;do echo "TEST DIAGO davidson in parallel, nprocs=$i" mpirun -np $i ./HSolver_LCAO mpirun -np $i ./HSolver_dav_real + mpirun -np $i ./HSolver_cg_real if [[ $? != 0 ]];then echo -e "\e[1;33m [ FAILED ] \e[0m"\ "execute UT with $i cores error."