From 0152800f6eeb007b6a1611566dffa3137d41ff2f Mon Sep 17 00:00:00 2001 From: Critsium Date: Fri, 31 Jul 2026 15:55:48 +0800 Subject: [PATCH] Refactor: inject nlocal/nbands into the LCAO dense diagonalizers Group D of the source_hsolver PARAM removal, and the largest single cluster: the six dense LCAO diagonalizers all sized their eigenvalue buffers from PARAM.globalv.nlocal and copied out PARAM.inp.nbands eigenvalues. Each of DiagoLapack, DiagoScalapack, DiagoElpa, DiagoElpaNative, DiagoCusolver and DiagoCusolverMP now takes (nlocal, nbands) through its constructor and stores them as members. DiagoElpaNative additionally takes use_gpu, replacing its PARAM.inp.device == "gpu" test. HSolverLCAO gains nbands and use_gpu (it already received nlocal in #7711) and supplies all eleven construction sites across hamiltSolvePsiK, parakSolve and parakSolve_cusolver. The three call contexts agree on both values: parakSolve distributes the same nlocal x nlocal matrix over a smaller pool grid, and parakSolve_cusolver gathers it per k-point, so a single pair covers all of them. nbands is injected rather than read from ParaV->get_nbands(). The two are equal for the LCAO path (LCAO_init_basis.cpp seeds ParaV from PARAM.inp.nbands), but module_lr sets paraMat_.nbands to nocc + nvirt, so sourcing it from ParaV would plant a trap for any future caller. The existing ParaV->get_nbands() reads inside parakSolve are left alone; they are not PARAM reads and are out of scope. In diago_scalapack.cpp, four call sites pass the address of the value to the Fortran routines. PARAM.globalv returns a const reference, so &PARAM.globalv.nlocal was already a const int*, matching pdsygvx_/pzhegvx_'s const int* n; &this->nlocal has the same type. Also resolves the two error-message strings in diago_lapack.cpp and diago_scalapack.cpp that embedded "PARAM.inp.nbands = ". These were deferred from #7706 because removing them needed exactly the nbands member added here. Two dead commented-out ELPA_Solver calls in diago_elpa.cpp are removed rather than updated; each duplicated the live statement two lines below it. The three LCAO diagonalizer tests supply their own nlocal/nbands, which they already read from the H/S input files, so their PARAM.sys.nlocal / PARAM.input.nbands writes and the now-dead `#define private public` include blocks are dropped. PARAM occurrences in source_hsolver production code: 60 -> 5. What remains is DiagoIterAssist's basis_type/calculation switches (group E) and one comment in simple_pexsi.cpp. Co-Authored-By: Claude Opus 5 --- source/source_esolver/esolver_ks_lcao.cpp | 4 +- .../source_esolver/esolver_ks_lcao_tddft.cpp | 4 +- source/source_hsolver/diago_cusolver.cpp | 7 ++-- source/source_hsolver/diago_cusolver.h | 9 ++++- source/source_hsolver/diago_cusolvermp.cpp | 5 +-- source/source_hsolver/diago_cusolvermp.h | 8 +++- source/source_hsolver/diago_elpa.cpp | 29 ++++++--------- source/source_hsolver/diago_elpa.h | 7 ++++ source/source_hsolver/diago_elpa_native.cpp | 9 ++--- source/source_hsolver/diago_elpa_native.h | 10 +++++ source/source_hsolver/diago_lapack.cpp | 33 ++++++++--------- source/source_hsolver/diago_lapack.h | 8 ++++ source/source_hsolver/diago_scalapack.cpp | 37 +++++++++---------- source/source_hsolver/diago_scalapack.h | 7 ++++ source/source_hsolver/hsolver_lcao.cpp | 22 +++++------ source/source_hsolver/hsolver_lcao.h | 11 ++++-- .../source_hsolver/test/diago_lapack_test.cpp | 12 +----- .../test/diago_lcao_cusolver_test.cpp | 9 +---- .../source_hsolver/test/diago_lcao_test.cpp | 11 ++---- source/source_lcao/LCAO_set.cpp | 4 +- .../module_deltaspin/cal_mw_from_lambda.cpp | 4 +- 21 files changed, 137 insertions(+), 113 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 7b44a2758c..7ebeda0f4c 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -438,7 +438,9 @@ void ESolver_KS_LCAO::hamilt2rho_single(UnitCell& ucell, int istep, int PARAM.inp.ks_solver, PARAM.globalv.kpar_lcao, PARAM.globalv.nlocal, - PARAM.inp.nelec); + PARAM.inp.nbands, + PARAM.inp.nelec, + PARAM.inp.device == "gpu"); hsolver_lcao_obj.solve(static_cast*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, this->chr, PARAM.inp.nspin, skip_charge); } diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index de35a40cf8..4ba1c5fb87 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -333,7 +333,9 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, PARAM.inp.ks_solver, PARAM.globalv.kpar_lcao, PARAM.globalv.nlocal, - PARAM.inp.nelec); + PARAM.inp.nbands, + PARAM.inp.nelec, + PARAM.inp.device == "gpu"); hsolver_lcao_obj.solve(static_cast>*>(this->p_hamilt), this->psi[0], this->pelec, diff --git a/source/source_hsolver/diago_cusolver.cpp b/source/source_hsolver/diago_cusolver.cpp index b02c413257..34ced44cb8 100644 --- a/source/source_hsolver/diago_cusolver.cpp +++ b/source/source_hsolver/diago_cusolver.cpp @@ -6,7 +6,6 @@ #include "source_base/module_external/scalapack_connector.h" #include "source_base/tool_title.h" #include "source_base/timer.h" -#include "source_io/module_parameter/parameter.h" #include #include @@ -22,7 +21,7 @@ template int DiagoCusolver::DecomposedState = 0; template -DiagoCusolver::DiagoCusolver() +DiagoCusolver::DiagoCusolver(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) { } @@ -42,13 +41,13 @@ void DiagoCusolver::diag( ModuleBase::TITLE("DiagoCusolver", "diag"); ModuleBase::timer::start("DiagoCusolver", "cusolver"); // Allocate memory for eigenvalues - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); std::vector eigenvectors(h_mat.row * h_mat.col); this->dc.Dngvd(h_mat.row, h_mat.col, h_mat.p, s_mat.p, eigen.data(), eigenvectors.data()); const int size = psi.get_nbands() * psi.get_nbasis(); BlasConnector::copy(size, eigenvectors.data(), 1, psi.get_pointer(), 1); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); ModuleBase::timer::end("DiagoCusolver", "cusolver"); } diff --git a/source/source_hsolver/diago_cusolver.h b/source/source_hsolver/diago_cusolver.h index 28ebdc852c..50d9a393c6 100644 --- a/source/source_hsolver/diago_cusolver.h +++ b/source/source_hsolver/diago_cusolver.h @@ -19,9 +19,11 @@ class DiagoCusolver public: - DiagoCusolver(); + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + DiagoCusolver(const int nlocal_in, const int nbands_in); ~DiagoCusolver(); - + // Override the diag function for CUSOLVER diagonalization void diag( hamilt::MatrixBlock& h_mat, @@ -40,6 +42,9 @@ class DiagoCusolver // Function to check if ELPA handle needs to be created or reused in MPI settings bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const; #endif + + const int nlocal; + const int nbands; }; } // namespace hsolver diff --git a/source/source_hsolver/diago_cusolvermp.cpp b/source/source_hsolver/diago_cusolvermp.cpp index b5bb80eb5f..395add83d2 100644 --- a/source/source_hsolver/diago_cusolvermp.cpp +++ b/source/source_hsolver/diago_cusolvermp.cpp @@ -1,6 +1,5 @@ #ifdef __CUSOLVERMP -#include "source_io/module_parameter/parameter.h" #include "diago_cusolvermp.h" #include "source_base/module_external/blas_connector.h" @@ -18,7 +17,7 @@ void DiagoCusolverMP::diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* hamilt::MatrixBlock h_mat, s_mat; phm_in->matrix(h_mat, s_mat); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); std::vector eigenvectors(h_mat.row * h_mat.col); MPI_Comm COMM_DIAG = MPI_COMM_WORLD; // use all processes @@ -30,7 +29,7 @@ void DiagoCusolverMP::diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* ModuleBase::timer::end("DiagoCusolverMP", "Diag_CusolverMP_gvd"); } const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); const int size = psi.get_nbands() * psi.get_nbasis(); BlasConnector::copy(size, eigenvectors.data(), inc, psi.get_pointer(), inc); } diff --git a/source/source_hsolver/diago_cusolvermp.h b/source/source_hsolver/diago_cusolvermp.h index d173315e31..155a4bf94d 100644 --- a/source/source_hsolver/diago_cusolvermp.h +++ b/source/source_hsolver/diago_cusolvermp.h @@ -16,11 +16,17 @@ class DiagoCusolverMP using Real = typename GetTypeReal::type; public: - DiagoCusolverMP() + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + DiagoCusolverMP(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) { } // the diag function for CUSOLVERMP diagonalization void diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* eigenvalue_in); + + private: + const int nlocal; + const int nbands; }; } // namespace hsolver #endif // __CUSOLVERMP diff --git a/source/source_hsolver/diago_elpa.cpp b/source/source_hsolver/diago_elpa.cpp index a264a41b9c..55ca982539 100644 --- a/source/source_hsolver/diago_elpa.cpp +++ b/source/source_hsolver/diago_elpa.cpp @@ -2,7 +2,6 @@ #include "source_base/global_function.h" #include "source_base/module_external/blas_connector.h" -#include "source_io/module_parameter/parameter.h" #include "module_genelpa/elpa_solver.h" #include "source_base/module_external/blacs_connector.h" #include "source_base/global_variable.h" @@ -75,13 +74,13 @@ void DiagoElpa>::diag( matcd h_mat, s_mat; phm_in->matrix(h_mat, s_mat); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); bool isReal = false; MPI_Comm COMM_DIAG = setmpicomm(); // set mpi_comm needed ELPA_Solver es((const bool)isReal, COMM_DIAG, - (const int)PARAM.inp.nbands, + (const int)this->nbands, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc); @@ -97,7 +96,7 @@ void DiagoElpa>::diag( es.exit(); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); #else ModuleBase::WARNING_QUIT("DiagoElpa", "DiagoElpa only can be used with macro __MPI"); @@ -113,15 +112,13 @@ void DiagoElpa::diag(hamilt::Hamilt* phm_in, matd h_mat, s_mat; phm_in->matrix(h_mat, s_mat); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); bool isReal = true; MPI_Comm COMM_DIAG = setmpicomm(); // set mpi_comm needed - // ELPA_Solver es(isReal, COMM_DIAG, PARAM.inp.nbands, h_mat.row, h_mat.col, - // h_mat.desc); ELPA_Solver es((const bool)isReal, COMM_DIAG, - (const int)PARAM.inp.nbands, + (const int)this->nbands, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc); @@ -135,7 +132,7 @@ void DiagoElpa::diag(hamilt::Hamilt* phm_in, es.exit(); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); #else ModuleBase::WARNING_QUIT("DiagoElpa", "DiagoElpa only can be used with macro __MPI"); @@ -151,11 +148,11 @@ void DiagoElpa>::diag_pool(hamilt::MatrixBlock eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); bool isReal = false; ELPA_Solver es((const bool)isReal, comm, - (const int)PARAM.inp.nbands, + (const int)this->nbands, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc); @@ -170,7 +167,7 @@ void DiagoElpa>::diag_pool(hamilt::MatrixBlocknbands, eigen.data(), inc, eigenvalue_in, inc); } template <> @@ -180,14 +177,12 @@ void DiagoElpa::diag_pool(hamilt::MatrixBlock& h_mat, Real* eigenvalue_in, MPI_Comm& comm) { - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); bool isReal = true; - // ELPA_Solver es(isReal, COMM_DIAG, PARAM.inp.nbands, h_mat.row, h_mat.col, - // h_mat.desc); ELPA_Solver es((const bool)isReal, comm, - (const int)PARAM.inp.nbands, + (const int)this->nbands, (const int)h_mat.row, (const int)h_mat.col, (const int*)h_mat.desc); @@ -203,7 +198,7 @@ void DiagoElpa::diag_pool(hamilt::MatrixBlock& h_mat, const int inc = 1; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "K-S equation was solved by genelpa2"); - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "eigenvalues were copied to ekb"); } diff --git a/source/source_hsolver/diago_elpa.h b/source/source_hsolver/diago_elpa.h index c6050af224..2bee7a14fa 100644 --- a/source/source_hsolver/diago_elpa.h +++ b/source/source_hsolver/diago_elpa.h @@ -15,6 +15,10 @@ class DiagoElpa using Real = typename GetTypeReal::type; public: + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + DiagoElpa(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {}; + void diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* eigenvalue_in); #ifdef __MPI // diagnolization used in parallel-k case @@ -30,6 +34,9 @@ class DiagoElpa bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const; static int lastmpinum; // last using mpi; #endif + + const int nlocal; + const int nbands; }; template diff --git a/source/source_hsolver/diago_elpa_native.cpp b/source/source_hsolver/diago_elpa_native.cpp index ee003a08e0..fb4cfbc448 100644 --- a/source/source_hsolver/diago_elpa_native.cpp +++ b/source/source_hsolver/diago_elpa_native.cpp @@ -4,7 +4,6 @@ #include "source_base/module_external/blas_connector.h" #include "source_base/module_external/blacs_connector.h" #include "source_base/global_variable.h" -#include "source_io/module_parameter/parameter.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_hsolver/module_genelpa/elpa_new.h" @@ -59,7 +58,7 @@ void DiagoElpaNative::diag_pool(hamilt::MatrixBlock& h_mat, ModuleBase::timer::start("DiagoElpaNative", "elpa_solve"); - int nev = PARAM.inp.nbands; + int nev = this->nbands; int narows = h_mat.row; int nacols = h_mat.col; @@ -70,7 +69,7 @@ void DiagoElpaNative::diag_pool(hamilt::MatrixBlock& h_mat, int nprows, npcols, myprow, mypcol; Cblacs_gridinfo(cblacs_ctxt, &nprows, &npcols, &myprow, &mypcol); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); std::vector eigenvectors(narows * nacols); if (elpa_init(20210430) != ELPA_OK) @@ -107,7 +106,7 @@ void DiagoElpaNative::diag_pool(hamilt::MatrixBlock& h_mat, #define ELPA_WITH_SYCL_GPU_VERSION 0 */ #if ELPA_WITH_NVIDIA_GPU_VERSION - if (PARAM.inp.device == "gpu") + if (this->use_gpu) { elpa_set(handle, "nvidia-gpu", 1, &success); elpa_set(handle, "real_kernel", ELPA_2STAGE_REAL_NVIDIA_GPU, &success); @@ -138,7 +137,7 @@ void DiagoElpaNative::diag_pool(hamilt::MatrixBlock& h_mat, } const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); const int size = psi.get_nbands() * psi.get_nbasis(); BlasConnector::copy(size, eigenvectors.data(), inc, psi.get_pointer(), inc); } diff --git a/source/source_hsolver/diago_elpa_native.h b/source/source_hsolver/diago_elpa_native.h index 2c556254e8..2cbd2c27d6 100644 --- a/source/source_hsolver/diago_elpa_native.h +++ b/source/source_hsolver/diago_elpa_native.h @@ -15,6 +15,12 @@ class DiagoElpaNative using Real = typename GetTypeReal::type; public: + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + /// @param use_gpu_in offload to the NVIDIA-GPU ELPA kernels when ELPA was built with GPU support + DiagoElpaNative(const int nlocal_in, const int nbands_in, const bool use_gpu_in) + : nlocal(nlocal_in), nbands(nbands_in), use_gpu(use_gpu_in) {}; + void diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* eigenvalue_in); #ifdef __MPI // diagnolization used in parallel-k case @@ -27,6 +33,10 @@ class DiagoElpaNative static int DecomposedState; + private: + const int nlocal; + const int nbands; + const bool use_gpu; }; template diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp index 6bdba41ae9..0020efc1e8 100644 --- a/source/source_hsolver/diago_lapack.cpp +++ b/source/source_hsolver/diago_lapack.cpp @@ -1,7 +1,5 @@ // Refactored according to diago_scalapack // This code will be futher refactored to remove the dependency of psi and hamilt -#include "source_io/module_parameter/parameter.h" - #include "diago_lapack.h" #include "source_base/global_variable.h" @@ -42,14 +40,14 @@ void DiagoLapack::diag(hamilt::Hamilt* phm_in, psi::Psi& phm_in->matrix(h_mat, s_mat); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); check_lapack_layout(h_mat, s_mat, eigen.size()); // Diag this->dsygvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); // Copy result const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } template <> @@ -61,11 +59,11 @@ void DiagoLapack>::diag(hamilt::Hamilt matcd h_mat, s_mat; phm_in->matrix(h_mat, s_mat); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); check_lapack_layout(h_mat, s_mat, eigen.size()); this->zhegvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } #ifdef __MPI @@ -78,11 +76,11 @@ void DiagoLapack>::diag(hamilt::Hamilt { ModuleBase::TITLE("DiagoLapack", "diag_pool"); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); check_lapack_layout(h_mat, s_mat, eigen.size()); this->dsygvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } template<> void DiagoLapack>::diag_pool(hamilt::MatrixBlock>& h_mat, @@ -93,11 +91,11 @@ void DiagoLapack>::diag(hamilt::Hamilt { ModuleBase::TITLE("DiagoLapack", "diag_pool"); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); check_lapack_layout(h_mat, s_mat, eigen.size()); this->zhegvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } #endif @@ -115,20 +113,20 @@ std::pair> DiagoLapack::dsygvx_once(const int ncol, memcpy(s_tmp.c, s_mat, sizeof(double) * ncol * nrow); const char jobz = 'V', range = 'I', uplo = 'U'; - const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1; + const int itype = 1, il = 1, iu = this->nbands, one = 1; int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0; double vl = 0, vu = 0; const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC; std::vector work(3, 0); std::vector iwork(1, 0); - std::vector ifail(PARAM.globalv.nlocal, 0); + std::vector ifail(this->nlocal, 0); std::vector iclustr(2 * GlobalV::DSIZE); std::vector gap(GlobalV::DSIZE); // LAPACK dsygvx signature: // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU, // ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO) - int n = PARAM.globalv.nlocal; + int n = this->nlocal; int lda = n, ldb = n, ldz = n; dsygvx_(&itype, &jobz, @@ -224,7 +222,7 @@ std::pair> DiagoLapack::zhegvx_once(const int ncol, memcpy(s_tmp.c, s_mat, sizeof(std::complex) * ncol * nrow); const char jobz = 'V', range = 'I', uplo = 'U'; - const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1; + const int itype = 1, il = 1, iu = this->nbands, one = 1; int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0; const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC; @@ -232,14 +230,14 @@ std::pair> DiagoLapack::zhegvx_once(const int ncol, std::vector> work(1, 0); std::vector rwork(3, 0); std::vector iwork(1, 0); - std::vector ifail(PARAM.globalv.nlocal, 0); + std::vector ifail(this->nlocal, 0); std::vector iclustr(2 * GlobalV::DSIZE); std::vector gap(GlobalV::DSIZE); // LAPACK zhegvx signature: // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU, // ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO) - int n = PARAM.globalv.nlocal; + int n = this->nlocal; int lda = n, ldb = n, ldz = n; zhegvx_(&itype, &jobz, @@ -416,8 +414,7 @@ void DiagoLapack::post_processing(const int info, const std::vector& vec { const std::string str_M = "M = " + ModuleBase::GlobalFunc::TO_STRING(vec[0]) + ".\n"; const std::string str_NZ = "NZ = " + ModuleBase::GlobalFunc::TO_STRING(vec[1]) + ".\n"; - const std::string str_NBANDS - = "PARAM.inp.nbands = " + ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.nbands) + ".\n"; + const std::string str_NBANDS = "nbands = " + ModuleBase::GlobalFunc::TO_STRING(this->nbands) + ".\n"; throw std::runtime_error(str_info_FILE + str_M + str_NZ + str_NBANDS); } else if (info / 16 % 2) diff --git a/source/source_hsolver/diago_lapack.h b/source/source_hsolver/diago_lapack.h index bfdf78ac34..9f674fc20f 100644 --- a/source/source_hsolver/diago_lapack.h +++ b/source/source_hsolver/diago_lapack.h @@ -26,6 +26,10 @@ class DiagoLapack using Real = typename GetTypeReal::type; public: + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + DiagoLapack(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {}; + void diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* eigenvalue_in); #ifdef __MPI // diagnolization used in parallel-k case @@ -61,6 +65,10 @@ class DiagoLapack int degeneracy_max = 12; // For reorthogonalized memory. 12 followes siesta. void post_processing(const int info, const std::vector& vec); + + private: + const int nlocal; + const int nbands; }; } // namespace hsolver diff --git a/source/source_hsolver/diago_scalapack.cpp b/source/source_hsolver/diago_scalapack.cpp index 366478fdc2..f2ddb8289b 100644 --- a/source/source_hsolver/diago_scalapack.cpp +++ b/source/source_hsolver/diago_scalapack.cpp @@ -1,6 +1,5 @@ //===================== // AUTHOR : Peize Lin -#include "source_io/module_parameter/parameter.h" // DATE : 2021-11-02 // REFACTORING AUTHOR : Daye Zheng // DATE : 2022-04-14 @@ -28,10 +27,10 @@ namespace hsolver matd h_mat, s_mat; phm_in->matrix(h_mat, s_mat); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); this->pdsygvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } template<> void DiagoScalapack>::diag(hamilt::Hamilt>* phm_in, psi::Psi>& psi, Real* eigenvalue_in) @@ -40,10 +39,10 @@ namespace hsolver matcd h_mat, s_mat; phm_in->matrix(h_mat, s_mat); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); this->pzhegvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } #ifdef __MPI @@ -56,10 +55,10 @@ namespace hsolver { ModuleBase::TITLE("DiagoScalapack", "diag_pool"); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); this->pdsygvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } template<> void DiagoScalapack>::diag_pool(hamilt::MatrixBlock>& h_mat, @@ -70,10 +69,10 @@ namespace hsolver { ModuleBase::TITLE("DiagoScalapack", "diag_pool"); assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc); - std::vector eigen(PARAM.globalv.nlocal, 0.0); + std::vector eigen(this->nlocal, 0.0); this->pzhegvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi); const int inc = 1; - BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc); + BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc); } #endif @@ -92,13 +91,13 @@ namespace hsolver memcpy(s_tmp.c, s_mat, sizeof(double) * ncol * nrow); const char jobz = 'V', range = 'I', uplo = 'U'; - const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1; + const int itype = 1, il = 1, iu = this->nbands, one = 1; int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0; double vl = 0, vu = 0; const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC; std::vector work(3, 0); std::vector iwork(1, 0); - std::vector ifail(PARAM.globalv.nlocal, 0); + std::vector ifail(this->nlocal, 0); std::vector iclustr(2 * GlobalV::DSIZE); std::vector gap(GlobalV::DSIZE); @@ -106,7 +105,7 @@ namespace hsolver &jobz, &range, &uplo, - &PARAM.globalv.nlocal, + &this->nlocal, h_tmp.c, &one, &one, @@ -152,7 +151,7 @@ namespace hsolver &jobz, &range, &uplo, - &PARAM.globalv.nlocal, + &this->nlocal, h_tmp.c, &one, &one, @@ -217,7 +216,7 @@ namespace hsolver memcpy(s_tmp.c, s_mat, sizeof(std::complex) * ncol * nrow); const char jobz = 'V', range = 'I', uplo = 'U'; - const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1; + const int itype = 1, il = 1, iu = this->nbands, one = 1; int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0; const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC; //Note: pzhegvx_ has a bug @@ -227,7 +226,7 @@ namespace hsolver std::vector> work(1, 0); std::vector rwork(3, 0); std::vector iwork(1, 0); - std::vector ifail(PARAM.globalv.nlocal, 0); + std::vector ifail(this->nlocal, 0); std::vector iclustr(2 * GlobalV::DSIZE); std::vector gap(GlobalV::DSIZE); @@ -235,7 +234,7 @@ namespace hsolver &jobz, &range, &uplo, - &PARAM.globalv.nlocal, + &this->nlocal, h_tmp.c, &one, &one, @@ -276,7 +275,7 @@ namespace hsolver // GlobalV::ofs_running<<"lwork="<degeneracy_max * PARAM.globalv.nlocal; + lrwork = rwork[0] + this->degeneracy_max * this->nlocal; int maxlrwork = std::max(lrwork,3); rwork.resize(maxlrwork, 0); liwork = iwork[0]; @@ -286,7 +285,7 @@ namespace hsolver &jobz, &range, &uplo, - &PARAM.globalv.nlocal, + &this->nlocal, h_tmp.c, &one, &one, @@ -431,7 +430,7 @@ namespace hsolver const std::string str_M = "M = " + ModuleBase::GlobalFunc::TO_STRING(vec[0]) + ".\n"; const std::string str_NZ = "NZ = " + ModuleBase::GlobalFunc::TO_STRING(vec[1]) + ".\n"; const std::string str_NBANDS - = "PARAM.inp.nbands = " + ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.nbands) + ".\n"; + = "nbands = " + ModuleBase::GlobalFunc::TO_STRING(this->nbands) + ".\n"; throw std::runtime_error(str_info_FILE + str_M + str_NZ + str_NBANDS); } else if (info / 16 % 2) diff --git a/source/source_hsolver/diago_scalapack.h b/source/source_hsolver/diago_scalapack.h index 05b1144836..7dce030616 100644 --- a/source/source_hsolver/diago_scalapack.h +++ b/source/source_hsolver/diago_scalapack.h @@ -27,6 +27,10 @@ namespace hsolver private: using Real = typename GetTypeReal::type; public: + /// @param nlocal_in global dimension of the NAO Hamiltonian + /// @param nbands_in number of lowest eigenpairs to compute + DiagoScalapack(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {}; + void diag(hamilt::Hamilt* phm_in, psi::Psi& psi, Real* eigenvalue_in); #ifdef __MPI // diagnolization used in parallel-k case @@ -67,6 +71,9 @@ namespace hsolver int degeneracy_max = 12; // For reorthogonalized memory. 12 followes siesta. void post_processing(const int info, const std::vector &vec); + + const int nlocal; + const int nbands; }; } // namespace hsolver diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp index d31658d8ad..3c8bcc3ff0 100644 --- a/source/source_hsolver/hsolver_lcao.cpp +++ b/source/source_hsolver/hsolver_lcao.cpp @@ -141,19 +141,19 @@ void HSolverLCAO::hamiltSolvePsiK(hamilt::Hamilt* hm, psi::Psi& if (this->method == "scalapack_gvx") { #ifdef __MPI - DiagoScalapack sa; + DiagoScalapack sa(this->nlocal, this->nbands); sa.diag(hm, psi, eigenvalue); #endif } #ifdef __ELPA else if (this->method == "genelpa") { - DiagoElpa el; + DiagoElpa el(this->nlocal, this->nbands); el.diag(hm, psi, eigenvalue); } else if (this->method == "elpa") { - DiagoElpaNative el; + DiagoElpaNative el(this->nlocal, this->nbands, this->use_gpu); el.diag(hm, psi, eigenvalue); } #endif @@ -161,7 +161,7 @@ void HSolverLCAO::hamiltSolvePsiK(hamilt::Hamilt* hm, psi::Psi& else if (this->method == "cusolver") { // Note: This branch will only be executed in the single-process case - DiagoCusolver cu; + DiagoCusolver cu(this->nlocal, this->nbands); hamilt::MatrixBlock hk, sk; hm->matrix(hk, sk); cu.diag(hk, sk, psi, eigenvalue); @@ -169,14 +169,14 @@ void HSolverLCAO::hamiltSolvePsiK(hamilt::Hamilt* hm, psi::Psi& #ifdef __CUSOLVERMP else if (this->method == "cusolvermp") { - DiagoCusolverMP cm; + DiagoCusolverMP cm(this->nlocal, this->nbands); cm.diag(hm, psi, eigenvalue); } #endif #endif else if (this->method == "lapack") // only for single core { - DiagoLapack la; + DiagoLapack la(this->nlocal, this->nbands); la.diag(hm, psi, eigenvalue); } else @@ -253,23 +253,23 @@ void HSolverLCAO::parakSolve(hamilt::Hamilt* pHamilt, /// solve eigenvector and eigenvalue for H(k) if (this->method == "scalapack_gvx") { - DiagoScalapack sa; + DiagoScalapack sa(this->nlocal, this->nbands); sa.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D); } else if (this->method == "lapack") { - DiagoLapack la; + DiagoLapack la(this->nlocal, this->nbands); la.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D); } #ifdef __ELPA else if (this->method == "genelpa") { - DiagoElpa el; + DiagoElpa el(this->nlocal, this->nbands); el.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D); } else if (this->method == "elpa") { - DiagoElpaNative el; + DiagoElpaNative el(this->nlocal, this->nbands, this->use_gpu); el.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D); } #endif @@ -427,7 +427,7 @@ void HSolverLCAO::parakSolve_cusolver(hamilt::Hamilt* pHamilt, if(kpt_assigned != -1) { psi_local.resize(1, ncol, nrow); - DiagoCusolver cu{}; + DiagoCusolver cu(this->nlocal, this->nbands); hamilt::MatrixBlock hk_local = hamilt::MatrixBlock{ hk_mat.data(), (size_t)nrow, (size_t)ncol, mat_para_local.desc}; diff --git a/source/source_hsolver/hsolver_lcao.h b/source/source_hsolver/hsolver_lcao.h index 60a777e3a4..3b54d92d1f 100644 --- a/source/source_hsolver/hsolver_lcao.h +++ b/source/source_hsolver/hsolver_lcao.h @@ -19,8 +19,11 @@ class HSolverLCAO const std::string method_in, const int kpar_lcao_in, const int nlocal_in, - const double nelec_in) - : ParaV(ParaV_in), method(method_in), kpar_lcao(kpar_lcao_in), nlocal(nlocal_in), nelec(nelec_in) {}; + const int nbands_in, + const double nelec_in, + const bool use_gpu_in) + : ParaV(ParaV_in), method(method_in), kpar_lcao(kpar_lcao_in), nlocal(nlocal_in), nbands(nbands_in), + nelec(nelec_in), use_gpu(use_gpu_in) {}; void solve(hamilt::Hamilt* pHamilt, psi::Psi& psi, @@ -49,8 +52,10 @@ class HSolverLCAO const std::string method; const int kpar_lcao; // number of pools for LCAO diagonalization - const int nlocal; // global dimension of the NAO Hamiltonian, only used by the pexsi branch + const int nlocal; // global dimension of the NAO Hamiltonian + const int nbands; // number of bands to be solved for const double nelec; // total number of electrons, only used by the pexsi branch + const bool use_gpu; // true if running on GPU, only used by the native-ELPA branch }; } // namespace hsolver diff --git a/source/source_hsolver/test/diago_lapack_test.cpp b/source/source_hsolver/test/diago_lapack_test.cpp index 8620eaafa4..e3d358052e 100644 --- a/source/source_hsolver/test/diago_lapack_test.cpp +++ b/source/source_hsolver/test/diago_lapack_test.cpp @@ -1,8 +1,5 @@ // Author: Zhang Xiaoyang // A modified version of diago_lcao_test.cpp -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private // Remove some useless functions and dependencies. Serialized the full code // and refactored some function. @@ -184,21 +181,14 @@ class DiagoLapackPrepare hmtest.ncol = nlocal; } - void set_env() - { - PARAM.sys.nlocal = nlocal; - PARAM.input.nbands = nbands; - } - void diago() { this->pb2d(); this->print_hs(); - this->set_env(); for (int i = 0; i < REPEATRUN; i++) { - hsolver::DiagoLapack dh; + hsolver::DiagoLapack dh(nlocal, nbands); dh.diag(&hmtest, psi, e_solver.data()); // dh->diag(&hmtest, psi, e_solver.data()); } diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp index 9f340419b7..3a5154af10 100644 --- a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp +++ b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp @@ -1,8 +1,5 @@ #include "source_hsolver/diago_scalapack.h" #include "source_hsolver/test/diago_elpa_utils.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private #include "mpi.h" #include "string.h" @@ -204,8 +201,6 @@ class DiagoPrepare void set_env() { - PARAM.sys.nlocal = nlocal; - PARAM.input.nbands = nbands; GlobalV::DSIZE = dsize; } @@ -225,13 +220,13 @@ class DiagoPrepare hmtest.s_local = this->s_local; if (ks_solver == "scalapack_gvx") { - hsolver::DiagoScalapack dh; + hsolver::DiagoScalapack dh(nlocal, nbands); dh.diag(&hmtest, psi, e_solver.data()); } #ifdef __CUDA else if (ks_solver == "cusolver") { - hsolver::DiagoCusolver dh; + hsolver::DiagoCusolver dh(nlocal, nbands); hamilt::MatrixBlock h_mat, s_mat; hmtest.matrix(h_mat, s_mat); dh.diag(h_mat, s_mat, psi, e_solver.data()); diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp index d0ed0fdae6..6d66269415 100644 --- a/source/source_hsolver/test/diago_lcao_test.cpp +++ b/source/source_hsolver/test/diago_lcao_test.cpp @@ -1,9 +1,6 @@ #include "source_hsolver/diago_scalapack.h" #include "source_hsolver/diago_lapack.h" #include "source_hsolver/test/diago_elpa_utils.h" -#define private public -#include "source_io/module_parameter/parameter.h" -#undef private #include "mpi.h" #include "string.h" @@ -205,8 +202,6 @@ class DiagoPrepare void set_env() { - PARAM.sys.nlocal = nlocal; - PARAM.input.nbands = nbands; GlobalV::DSIZE = dsize; } @@ -226,18 +221,18 @@ class DiagoPrepare hmtest.s_local = this->s_local; if (ks_solver == "scalapack_gvx") { - hsolver::DiagoScalapack dh; + hsolver::DiagoScalapack dh(nlocal, nbands); dh.diag(&hmtest, psi, e_solver.data()); } else if (ks_solver == "lapack") { - hsolver::DiagoLapack la; + hsolver::DiagoLapack la(nlocal, nbands); la.diag(&hmtest, psi, e_solver.data()); } #ifdef __ELPA else if (ks_solver == "genelpa") { - hsolver::DiagoElpa dh; + hsolver::DiagoElpa dh(nlocal, nbands); dh.diag(&hmtest, psi, e_solver.data()); } #endif diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp index 584bfc6cb2..79ec359171 100644 --- a/source/source_lcao/LCAO_set.cpp +++ b/source/source_lcao/LCAO_set.cpp @@ -234,7 +234,9 @@ void LCAO_domain::init_chg_hr( ks_solver, PARAM.globalv.kpar_lcao, PARAM.globalv.nlocal, - PARAM.inp.nelec); + PARAM.inp.nbands, + PARAM.inp.nelec, + PARAM.inp.device == "gpu"); hsolver_lcao_obj.solve(p_hamilt, psi, pelec, dm, chr, nspin, 0); } diff --git a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp index 65cbd88608..fc2fdc06b3 100644 --- a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp +++ b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp @@ -525,7 +525,9 @@ void spinconstrain::SpinConstrain>::cal_mw_from_lambda( PARAM.inp.ks_solver, PARAM.globalv.kpar_lcao, PARAM.globalv.nlocal, - PARAM.inp.nelec); + PARAM.inp.nbands, + PARAM.inp.nelec, + PARAM.inp.device == "gpu"); if (this->nspin_ == 2) { dynamic_cast, double>>*>(this->p_operator)