Refactor: finish removing PARAM from source_hsolver (dense diagonalizers + DiagoIterAssist)

source/source_esolver/esolver_ks_lcao.cpp

@@ -450,7 +450,9 @@ void ESolver_KS_LCAO<TK, TR>::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<hamilt::Hamilt<TK>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, 
           this->chr, PARAM.inp.nspin, skip_charge);
     }

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -351,7 +351,9 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::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<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                    this->psi[0],
                                    this->pelec,

source/source_esolver/esolver_ks_lcaopw.cpp

@@ -135,7 +135,10 @@ namespace ModuleESolver
         hsolver::DiagoIterAssist<T>::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax;
         bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
-        hsolver::HSolverLIP<T> hsolver_lip_obj(this->pw_wfc, PARAM.globalv.use_uspp);
+        hsolver::HSolverLIP<T> hsolver_lip_obj(this->pw_wfc,
+                                               PARAM.globalv.use_uspp,
+                                               PARAM.inp.basis_type,
+                                               PARAM.inp.calculation);
         hsolver_lip_obj.solve(static_cast<hamilt::Hamilt<T>*>(this->p_hamilt), *this->stp.template get_psi_t<T, base_device::DEVICE_CPU>(), this->pelec, 
           *this->psi_local, skip_charge,ucell.tpiba,ucell.nat);
 

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 <memory>
 #include <type_traits>
@@ -22,7 +21,7 @@ template <typename T>
 int DiagoCusolver<T>::DecomposedState = 0;
 
 template <typename T>
-DiagoCusolver<T>::DiagoCusolver()
+DiagoCusolver<T>::DiagoCusolver(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in)
 {
 }
 
@@ -42,13 +41,13 @@ void DiagoCusolver<T>::diag(
     ModuleBase::TITLE("DiagoCusolver", "diag");
     ModuleBase::timer::start("DiagoCusolver", "cusolver");
     // Allocate memory for eigenvalues
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     std::vector<T> 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");
 }
 

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<T>& 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

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<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real*
     hamilt::MatrixBlock<T> h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<Real> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(h_mat.row * h_mat.col);
 
     MPI_Comm COMM_DIAG = MPI_COMM_WORLD; // use all processes
@@ -30,7 +29,7 @@ void DiagoCusolverMP<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& 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);
 }

source/source_hsolver/diago_cusolvermp.h

@@ -16,11 +16,17 @@ class DiagoCusolverMP
     using Real = typename GetTypeReal<T>::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<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
+
+  private:
+    const int nlocal;
+    const int nbands;
 };
 } // namespace hsolver
 #endif // __CUSOLVERMP

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<std::complex<double>>::diag(
     matcd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> 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<std::complex<double>>::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<double>::diag(hamilt::Hamilt<double>* phm_in,
     matd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> 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<double>::diag(hamilt::Hamilt<double>* 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<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex
     Real* eigenvalue_in,
     MPI_Comm& comm)
 {
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> 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<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex
     ModuleBase::timer::end("DiagoElpa", "elpa_solve");
     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);
 }
 
 template <>
@@ -180,14 +177,12 @@ void DiagoElpa<double>::diag_pool(hamilt::MatrixBlock<double>& h_mat,
     Real* eigenvalue_in,
     MPI_Comm& comm)
 {
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> 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<double>::diag_pool(hamilt::MatrixBlock<double>& 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");
 }

source/source_hsolver/diago_elpa.h

@@ -15,6 +15,10 @@ class DiagoElpa
     using Real = typename GetTypeReal<T>::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<T>* phm_in, psi::Psi<T>& 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 <typename T>

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<T>::diag_pool(hamilt::MatrixBlock<T>& 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<T>::diag_pool(hamilt::MatrixBlock<T>& h_mat,
     int nprows, npcols, myprow, mypcol;
 
     Cblacs_gridinfo(cblacs_ctxt, &nprows, &npcols, &myprow, &mypcol);
-    std::vector<Real> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(narows * nacols);
 
     if (elpa_init(20210430) != ELPA_OK)
@@ -107,7 +106,7 @@ void DiagoElpaNative<T>::diag_pool(hamilt::MatrixBlock<T>& 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<T>::diag_pool(hamilt::MatrixBlock<T>& 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);
 }

source/source_hsolver/diago_elpa_native.h

@@ -15,6 +15,12 @@ class DiagoElpaNative
     using Real = typename GetTypeReal<T>::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<T>* phm_in, psi::Psi<T>& 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 <typename T>

source/source_hsolver/diago_iter_assist.cpp

@@ -1,5 +1,4 @@
 #include "diago_iter_assist.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_base/complexmatrix.h"
 #include "source_base/constants.h"
 #include "source_base/global_function.h"
@@ -177,6 +176,8 @@ void DiagoIterAssist<T, Device>::diag_subspace_init(hamilt::Hamilt<T, Device>* p
     int psi_nc,
     psi::Psi<T, Device>& evc,
     Real* en,
+    const std::string& basis_type,
+    const std::string& calculation,
     const std::function<void(T*, const int)>& add_to_hcc,
     const std::function<void(const T* const, const int, const int)>& export_vcc)
 {
@@ -330,13 +331,13 @@ void DiagoIterAssist<T, Device>::diag_subspace_init(hamilt::Hamilt<T, Device>* p
     //=======================
     // diagonize the H-matrix
     //=======================
-    if ((PARAM.inp.basis_type == "lcao" || PARAM.inp.basis_type == "lcao_in_pw") && PARAM.inp.calculation == "nscf")
+    if ((basis_type == "lcao" || basis_type == "lcao_in_pw") && calculation == "nscf")
     {
         GlobalV::ofs_running << " Not do zgemm to get evc." << std::endl;
     }
-    else if ((PARAM.inp.basis_type == "lcao" || PARAM.inp.basis_type == "lcao_in_pw" || PARAM.inp.basis_type == "pw")
-             && (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "md"
-                 || PARAM.inp.calculation == "relax")) // pengfei 2014-10-13
+    else if ((basis_type == "lcao" || basis_type == "lcao_in_pw" || basis_type == "pw")
+             && (calculation == "scf" || calculation == "md"
+                 || calculation == "relax")) // pengfei 2014-10-13
     {
         // because psi and evc are different here,
         // I think if psi and evc are the same,
@@ -637,31 +638,6 @@ void DiagoIterAssist<T, Device>::diag_subspace_psi(const T* hcc,
     ModuleBase::timer::end("DiagoIterAssist", "diag_subspace_psi");
 }
 
-template <typename T, typename Device>
-bool DiagoIterAssist<T, Device>::test_exit_cond(const int& ntry, const int& notconv)
-{
-    //================================================================
-    // If this logical function is true, need to do diag_subspace
-    // and cg again.
-    //================================================================
-
-    bool scf = true;
-    if (PARAM.inp.calculation == "nscf") {
-        scf = false;
-}
-
-    // If ntry <=5, try to do it better, if ntry > 5, exit.
-    const bool f1 = (ntry <= 5);
-
-    // In non-self consistent calculation, do until totally converged.
-    const bool f2 = ((!scf && (notconv > 0)));
-
-    // if self consistent calculation, if not converged > 5,
-    // using diag_subspace and cg method again. ntry++
-    const bool f3 = ((scf && (notconv > 5)));
-    return (f1 && (f2 || f3));
-}
-
 template class DiagoIterAssist<std::complex<float>, base_device::DEVICE_CPU>;
 template class DiagoIterAssist<std::complex<double>, base_device::DEVICE_CPU>;
 #if ((defined __CUDA) || (defined __ROCM))