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

* 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 <noreply@anthropic.com>

* Refactor: pass basis_type/calculation into diag_subspace_init, drop dead test_exit_cond

Group E, the last PARAM cluster in source_hsolver. Both remaining reads were in
DiagoIterAssist, which is an all-static class.

diag_subspace_init used PARAM.inp.basis_type and PARAM.inp.calculation to pick
how the rotation matrix is applied to psi (skip the gemm entirely, apply it over
dmax rows, or over dmin rows). Both are now parameters. They are threaded rather
than added as static members: the class already carries mutable static config
(PW_DIAG_THR, SCF_ITER, ...), and coding rule 2 in AGENTS.md discourages adding
more workflow switches that can be set from several places. There are only two
call sites, so parameters are cheap here.

- HSolverLIP gains basis_type and calculation, mirroring HSolverPW which already
  receives both.
- PSIPrepare already holds basis_type as a member; calculation is taken from
  PARAM at that call site. source_psi is outside this refactor and already reads
  PARAM in the same function, so this introduces no new dependency there -- it
  just stops source_hsolver from reaching for the global itself.

test_exit_cond is deleted rather than converted. It has no caller anywhere in
the tree: DiagoCG and Diago_DavSubspace each define their own independent
test_exit_cond, and Diago_DavSubspace's already takes scf as an argument, which
is the shape this one would have needed. Its only apparent user is
test/test_diago_assist.cpp, which no CMakeLists references and which could not
compile if it were built (it includes a misspelled "diago_iter_assis.h", places
a bare statement in a class body, and calls EXPECT_EQ with one argument). That
file is left untouched -- it was already dead before this change, and removing
it is a separate call for the maintainers.

PARAM occurrences in source_hsolver production code: 5 -> 0. The only remaining
match in non-test code is a comment in simple_pexsi.cpp recording that a
hardcoded 2 should eventually become nspin.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* chore: drop accidentally committed CI notes from the branch

CI提速方案.md is a local working note that was swept into the previous commit
by a `git add -A`. It is unrelated to this PR's refactor and should not be part
of the tree; removing it here rather than rewriting history, since a merge of
develop has since landed on top of that commit.

The file is untracked, not deleted, so it stays in the local working copy.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>

Author: Critsium-xy

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);
 }