Refactor hsolver: drop dead global includes and GlobalV::DSIZE (#7865)

Step 1 — remove includes of global_variable.h / parameter.h from files
that no longer reference any symbol from them: diago_cusolver.cpp,
diago_pexsi.cpp, diago_elpa_native.cpp and module_pexsi/simple_pexsi.cpp.

Step 2 — remove GlobalV::DSIZE from the ScaLAPACK/LAPACK solvers:

* diago_scalapack.cpp: p?sygvx needs iclustr(2*NPROW*NPCOL) and
  gap(NPROW*NPCOL) of the grid the descriptor lives on, which is not
  necessarily the size of the diag world.  Query the grid via
  Cblacs_gridinfo() on desc[1] instead.  This also fixes a latent
  out-of-bounds read in post_processing() whenever the BLACS grid is
  smaller than DSIZE; the loop now derives its bound from iclustr.

* diago_lapack.cpp: the serial LAPACK ?sygvx/?hegvx take no iclustr/gap
  at all, so gap is dropped and iclustr becomes a single-process
  placeholder.

No functional change for the existing single-grid runs.
Verified with MODULE_HSOLVER_LCAO (serial and 4 MPI ranks) and
MODULE_HSOLVER_diago_hs_parallel.

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

Author: Critsium-xy

source/source_hsolver/diago_cusolver.cpp

@@ -2,7 +2,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_base/module_external/scalapack_connector.h"
 #include "source_base/tool_title.h"
 #include "source_base/timer.h"

source/source_hsolver/diago_elpa_native.cpp

@@ -3,7 +3,6 @@
 #include "source_base/global_function.h"
 #include "source_base/module_external/blas_connector.h"
 #include "source_base/module_external/blacs_connector.h"
-#include "source_base/global_variable.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_hsolver/module_genelpa/elpa_new.h"

source/source_hsolver/diago_lapack.cpp

@@ -120,8 +120,10 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(this->nlocal, 0);
-    std::vector<int> iclustr(2 * GlobalV::DSIZE);
-    std::vector<double> gap(GlobalV::DSIZE);
+    // iclustr/gap are ScaLAPACK-only outputs; the serial LAPACK call below never
+    // writes them, so gap is dropped and iclustr is kept merely as a zero-filled
+    // single-process placeholder for the info & 2 branch of post_processing().
+    std::vector<int> iclustr(2, 0);
 
     // LAPACK dsygvx signature:
     // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU,
@@ -231,8 +233,10 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     std::vector<double> rwork(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(this->nlocal, 0);
-    std::vector<int> iclustr(2 * GlobalV::DSIZE);
-    std::vector<double> gap(GlobalV::DSIZE);
+    // iclustr/gap are ScaLAPACK-only outputs; the serial LAPACK call below never
+    // writes them, so gap is dropped and iclustr is kept merely as a zero-filled
+    // single-process placeholder for the info & 2 branch of post_processing().
+    std::vector<int> iclustr(2, 0);
 
     // LAPACK zhegvx signature:
     // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU,
@@ -393,7 +397,8 @@ void DiagoLapack<T>::post_processing(const int info, const std::vector<int>& vec
     else if (info / 2 % 2)
     {
         int degeneracy_need = 0;
-        for (int irank = 0; irank < GlobalV::DSIZE; ++irank) {
+        // `vec` is iclustr, holding one [begin, end] pair per process
+        for (std::size_t irank = 0; 2 * irank + 1 < vec.size(); ++irank) {
             degeneracy_need = std::max(degeneracy_need, vec[2 * irank + 1] - vec[2 * irank]);
         }
         const std::string str_need = "degeneracy_need = " + ModuleBase::GlobalFunc::TO_STRING(degeneracy_need) + ".\n";

source/source_hsolver/diago_pexsi.cpp

@@ -4,7 +4,6 @@
 #ifdef __PEXSI
 #include "diago_pexsi.h"
 #include "source_base/tool_title.h"
-#include "source_base/global_variable.h"
 #include "source_base/tool_quit.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
 #include "module_pexsi/pexsi_solver.h"

source/source_hsolver/diago_scalapack.cpp

@@ -12,6 +12,7 @@
 
 #include "source_base/global_function.h"
 #include "source_base/global_variable.h"
+#include "source_base/module_external/blacs_connector.h"
 #include "source_base/module_external/scalapack_connector.h"
 #include "source_hamilt/matrixblock.h"
 
@@ -20,6 +21,18 @@ typedef hamilt::MatrixBlock<std::complex<double>> matcd;
 
 namespace hsolver
 {
+namespace
+{
+/// Number of processes in the BLACS process grid that `desc` lives on.
+/// p?sygvx requires iclustr(2*NPROW*NPCOL) and gap(NPROW*NPCOL).
+int blacs_grid_size(const int* const desc)
+{
+    int nprow = 0, npcol = 0, myprow = 0, mypcol = 0;
+    Cblacs_gridinfo(desc[1], &nprow, &npcol, &myprow, &mypcol);
+    return nprow * npcol;
+}
+} // namespace
+
     template<>
     void DiagoScalapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>& psi, Real* eigenvalue_in)
 {
@@ -98,8 +111,9 @@ namespace hsolver
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(this->nlocal, 0);
-    std::vector<int> iclustr(2 * GlobalV::DSIZE);
-    std::vector<double> gap(GlobalV::DSIZE);
+    const int ngrid = blacs_grid_size(desc);
+    std::vector<int> iclustr(2 * ngrid);
+    std::vector<double> gap(ngrid);
 
     pdsygvx_(&itype,
              &jobz,
@@ -227,8 +241,9 @@ namespace hsolver
     std::vector<double> rwork(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(this->nlocal, 0);
-    std::vector<int> iclustr(2 * GlobalV::DSIZE);
-    std::vector<double> gap(GlobalV::DSIZE);
+    const int ngrid = blacs_grid_size(desc);
+    std::vector<int> iclustr(2 * ngrid);
+    std::vector<double> gap(ngrid);
 
     pzhegvx_(&itype,
              &jobz,
@@ -408,7 +423,8 @@ namespace hsolver
     else if (info / 2 % 2)
     {
         int degeneracy_need = 0;
-        for (int irank = 0; irank < GlobalV::DSIZE; ++irank) {
+        // `vec` is iclustr, sized 2*NPROW*NPCOL by the caller
+        for (std::size_t irank = 0; 2 * irank + 1 < vec.size(); ++irank) {
             degeneracy_need = std::max(degeneracy_need, vec[2 * irank + 1] - vec[2 * irank]);
 }
         const std::string str_need = "degeneracy_need = " + ModuleBase::GlobalFunc::TO_STRING(degeneracy_need) + ".\n";

source/source_hsolver/module_pexsi/simple_pexsi.cpp

@@ -1,6 +1,5 @@
 // use PEXSI to solve a Kohn-Sham equation
 // the H and S matrices are given by 2D block cyclic distribution
-#include "source_io/module_parameter/parameter.h"
 // the Density Matrix and Energy Density Matrix calculated by PEXSI are transformed to 2D block cyclic distribution
 // #include "mpi.h"
 #ifdef __PEXSI
@@ -20,7 +19,6 @@
 #include "source_base/module_external/lapack_connector.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
-#include "source_base/global_variable.h"
 #include "source_hsolver/diago_pexsi.h"
 
 namespace pexsi