From d0820bf7aa5533ce2f8bf69ee72d50e3bd40b5f9 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 11:28:32 +0800
Subject: [PATCH 01/10] Refactor diago_lapack.cpp

---
 source/source_hsolver/diago_lapack.cpp | 69 +++++++++++++++++++++++---
 1 file changed, 63 insertions(+), 6 deletions(-)

diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp
index 90018c288d..a735edcc1d 100644
--- a/source/source_hsolver/diago_lapack.cpp
+++ b/source/source_hsolver/diago_lapack.cpp
@@ -142,11 +142,37 @@ int DiagoLapack<T>::dsygvx_once(const int ncol,
             ifail.data(),
             &info);*/
 
-    double *ev = new double[ncol * ncol];
+        // Use dsygvx_ to compute only the requested eigenvalues (range 'I').
+        int lwork2 = std::max(1, ncol * ncol);
+        std::vector<double> work_large(lwork2, 0.0);
+        int liwork = std::max(1, 5 * ncol);
+        std::vector<int> iwork_large(liwork, 0);
 
-    dsygv_(&itype, &jobz, &uplo, &PARAM.globalv.nlocal, h_tmp.c, &ncol, s_tmp.c, &ncol, ekb, ev, &lwork, &info);
+        dsygvx_(&itype,
+            &jobz,
+            &range,
+            &uplo,
+            &PARAM.globalv.nlocal,
+            h_tmp.c,
+            &ncol,
+            s_tmp.c,
+            &ncol,
+            &vl,
+            &vu,
+            &il,
+            &iu,
+            &abstol,
+            &M,
+            ekb,
+            wfc_2d.get_pointer(),
+            &ncol,
+            work_large.data(),
+            &lwork2,
+            iwork_large.data(),
+            ifail.data(),
+            &info);
 
-    return info;
+        return info;
 }
 
 template <typename T>
@@ -243,11 +269,42 @@ int DiagoLapack<T>::zhegvx_once(const int ncol,
 
     */
 
-    std::complex<double> *ev = new std::complex<double>[ncol * ncol];
+        // Use zhegvx_ to compute only the requested eigenvalues (range 'I').
+        // Allocate a reasonably large complex workspace and integer workspace.
+        int lwork2 = std::max(1, ncol * ncol);
+        std::complex<double> *workc = new std::complex<double>[lwork2];
+        int liwork = std::max(1, 5 * ncol);
+        std::vector<int> iwork_large(liwork, 0);
+
+        zhegvx_(&itype,
+            &jobz,
+            &range,
+            &uplo,
+            &PARAM.globalv.nlocal,
+            h_tmp.c,
+            &ncol,
+            s_tmp.c,
+            &ncol,
+            &vl,
+            &vu,
+            &il,
+            &iu,
+            &abstol,
+            &M,
+            ekb,
+            wfc_2d.get_pointer(),
+            &ncol,
+            workc,
+            &lwork2,
+            rwork,
+            iwork_large.data(),
+            ifail.data(),
+            &info);
 
-    zhegv_(&itype, &jobz, &uplo, &PARAM.globalv.nlocal, h_tmp.c, &ncol, s_tmp.c, &ncol, ekb, ev, &lwork, rwork, &info);
+        delete[] workc;
+        delete[] rwork;
 
-    return info;
+        return info;
 }
 
 template <typename T>

From 2512ca923c28d3d4d3436a1ffea02a00bfafd001 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 13:19:14 +0800
Subject: [PATCH 02/10] Totally refactor lapack sover

---
 source/source_hsolver/CMakeLists.txt   |   1 +
 source/source_hsolver/diago_lapack.cpp | 345 ++++++++++++++-----------
 source/source_hsolver/diago_lapack.h   |  28 +-
 source/source_hsolver/hsolver_lcao.cpp |   6 +-
 tests/performance/P100_si16_lcao/INPUT |   1 +
 5 files changed, 221 insertions(+), 160 deletions(-)

diff --git a/source/source_hsolver/CMakeLists.txt b/source/source_hsolver/CMakeLists.txt
index c2aa8c60e4..568510b603 100644
--- a/source/source_hsolver/CMakeLists.txt
+++ b/source/source_hsolver/CMakeLists.txt
@@ -21,6 +21,7 @@ if(ENABLE_LCAO)
         hsolver_lcao.cpp
         diago_scalapack.cpp
         parallel_k2d.cpp
+        diago_lapack.cpp
     )
   else ()
     list(APPEND objects
diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp
index a735edcc1d..f78298a488 100644
--- a/source/source_hsolver/diago_lapack.cpp
+++ b/source/source_hsolver/diago_lapack.cpp
@@ -30,12 +30,10 @@ void DiagoLapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>&
     // Diag
     this->dsygvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     // Copy result
-    int size = eigen.size();
-    for (int i = 0; i < size; i++)
-    {
-        eigenvalue_in[i] = eigen[i];
-    }
+    const int inc = 1;
+    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
+
 template <>
 void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>>* phm_in,
                                              psi::Psi<std::complex<double>>& psi,
@@ -48,52 +46,79 @@ void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>
 
     std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
     this->zhegvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
-    int size = eigen.size();
-    for (int i = 0; i < size; i++)
-    {
-        eigenvalue_in[i] = eigen[i];
-    }
+    const int inc = 1;
+    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 
+#ifdef __MPI
+ template<>
+    void DiagoLapack<double>::diag_pool(hamilt::MatrixBlock<double>& h_mat,
+    hamilt::MatrixBlock<double>& s_mat,
+    psi::Psi<double>& psi,
+    Real* eigenvalue_in,
+    MPI_Comm& comm)
+{
+    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<double> eigen(PARAM.globalv.nlocal, 0.0);
+    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);
+}
+    template<>
+    void DiagoLapack<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex<double>>& h_mat,
+    hamilt::MatrixBlock<std::complex<double>>& s_mat,
+    psi::Psi<std::complex<double>>& psi,
+    Real* eigenvalue_in,
+    MPI_Comm& comm)
+{
+    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<double> eigen(PARAM.globalv.nlocal, 0.0);
+    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);
+}
+#endif
+
 template <typename T>
-int DiagoLapack<T>::dsygvx_once(const int ncol,
+std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
                                 const int nrow,
                                 const double* const h_mat,
                                 const double* const s_mat,
                                 double* const ekb,
                                 psi::Psi<double>& wfc_2d) const
 {
-    // Copy matrix to temp variables
     ModuleBase::matrix h_tmp(ncol, nrow, false);
     memcpy(h_tmp.c, h_mat, sizeof(double) * ncol * nrow);
-
-
     ModuleBase::matrix s_tmp(ncol, nrow, false);
     memcpy(s_tmp.c, s_mat, sizeof(double) * ncol * nrow);
 
-    // Prepare caculate parameters
     const char jobz = 'V', range = 'I', uplo = 'U';
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
-    int M = 0, info = 0;
+    int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0;
     double vl = 0, vu = 0;
-    const double abstol = 0;
-
-    int lwork = (ncol + 2) * ncol;
-
+    const double abstol = 0, orfac = -1;
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(PARAM.globalv.nlocal, 0);
-
-    // Original Lapack caculate, obelsete
-    /*dsygvx_(&itype,
+    std::vector<int> iclustr(2 * GlobalV::DSIZE);
+    std::vector<double> 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 lda = n, ldb = n, ldz = n;
+        dsygvx_(&itype,
             &jobz,
             &range,
             &uplo,
-            &PARAM.globalv.nlocal,
+            &n,
             h_tmp.c,
-            &ncol,
+            &lda,
             s_tmp.c,
-            &ncol,
+            &ldb,
             &vl,
             &vu,
             &il,
@@ -102,31 +127,34 @@ int DiagoLapack<T>::dsygvx_once(const int ncol,
             &M,
             ekb,
             wfc_2d.get_pointer(),
-            &ncol,
+            &ldz,
             work.data(),
             &lwork,
             iwork.data(),
             ifail.data(),
             &info);
-
-    // Throw error if it returns info
-    if (info)
+    if (info) {
         throw std::runtime_error("info = " + ModuleBase::GlobalFunc::TO_STRING(info) + ".\n"
                                  + std::string(__FILE__) + " line "
                                  + std::to_string(__LINE__));
-    //lwork = work[0];
-    //work.resize(std::max(lwork, 3), 0);
-    //iwork.resize(iwork[0], 0);
+}
+
+    // Query returned optimal lwork in work[0]
+    lwork = static_cast<int>(work[0]);
+    work.resize(std::max(lwork, 3), 0);
+    // LAPACK integer workspace: use conservative size (5*N)
+    liwork = std::max(1, 5 * n);
+    iwork.resize(liwork, 0);
 
-    dsygvx_(&itype,
+        dsygvx_(&itype,
             &jobz,
             &range,
             &uplo,
-            &PARAM.globalv.nlocal,
+            &n,
             h_tmp.c,
-            &PARAM.globalv.nlocal,
+            &lda,
             s_tmp.c,
-            &PARAM.globalv.nlocal,
+            &ldb,
             &vl,
             &vu,
             &il,
@@ -135,48 +163,35 @@ int DiagoLapack<T>::dsygvx_once(const int ncol,
             &M,
             ekb,
             wfc_2d.get_pointer(),
-            &ncol,
+            &ldz,
             work.data(),
             &lwork,
             iwork.data(),
             ifail.data(),
-            &info);*/
-
-        // Use dsygvx_ to compute only the requested eigenvalues (range 'I').
-        int lwork2 = std::max(1, ncol * ncol);
-        std::vector<double> work_large(lwork2, 0.0);
-        int liwork = std::max(1, 5 * ncol);
-        std::vector<int> iwork_large(liwork, 0);
-
-        dsygvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &PARAM.globalv.nlocal,
-            h_tmp.c,
-            &ncol,
-            s_tmp.c,
-            &ncol,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ncol,
-            work_large.data(),
-            &lwork2,
-            iwork_large.data(),
-            ifail.data(),
             &info);
-
-        return info;
+    //	GlobalV::ofs_running<<"M="<<M<<"\t"<<"NZ="<<NZ<<std::endl;
+
+    if (info == 0) {
+        return std::make_pair(info, std::vector<int>{});
+    } else if (info < 0) {
+        return std::make_pair(info, std::vector<int>{});
+    } else if (info % 2) {
+        return std::make_pair(info, ifail);
+    } else if (info / 2 % 2) {
+        return std::make_pair(info, iclustr);
+    } else if (info / 4 % 2) {
+        return std::make_pair(info, std::vector<int>{M, NZ});
+    } else if (info / 16 % 2) {
+        return std::make_pair(info, ifail);
+    } else {
+        throw std::runtime_error("info = " + ModuleBase::GlobalFunc::TO_STRING(info) + ".\n"
+                                 + std::string(__FILE__) + " line "
+                                 + std::to_string(__LINE__));
+    }
 }
 
 template <typename T>
-int DiagoLapack<T>::zhegvx_once(const int ncol,
+std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
                                 const int nrow,
                                 const std::complex<double>* const h_mat,
                                 const std::complex<double>* const s_mat,
@@ -185,34 +200,38 @@ int DiagoLapack<T>::zhegvx_once(const int ncol,
 {
     ModuleBase::ComplexMatrix h_tmp(ncol, nrow, false);
     memcpy(h_tmp.c, h_mat, sizeof(std::complex<double>) * ncol * nrow);
-
     ModuleBase::ComplexMatrix s_tmp(ncol, nrow, false);
     memcpy(s_tmp.c, s_mat, sizeof(std::complex<double>) * ncol * nrow);
 
     const char jobz = 'V', range = 'I', uplo = 'U';
     const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
-    int M = 0, lrwork = -1, info = 0;
-    const double abstol = 0;
-
-    int lwork = (ncol + 2) * ncol;
-
+    int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0;
+    const double abstol = 0, orfac = -1;
+    //Note: pzhegvx_ has a bug
+    //      We must give vl,vu a value, although we do not use range 'V'
+    //      We must give rwork at least a memory of sizeof(double) * 3
     const double vl = 0, vu = 0;
     std::vector<std::complex<double>> work(1, 0);
-    double *rwork = new double[3 * ncol - 2];
+    std::vector<double> rwork(3, 0);
     std::vector<int> iwork(1, 0);
     std::vector<int> ifail(PARAM.globalv.nlocal, 0);
-
-    // Original Lapack caculate, obelsete
-    /*
-    zhegvx_(&itype,
+    std::vector<int> iclustr(2 * GlobalV::DSIZE);
+    std::vector<double> 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 lda = n, ldb = n, ldz = n;
+        zhegvx_(&itype,
             &jobz,
             &range,
             &uplo,
-            &PARAM.globalv.nlocal,
+            &n,
             h_tmp.c,
-            &PARAM.globalv.nlocal,
+            &lda,
             s_tmp.c,
-            &PARAM.globalv.nlocal,
+            &ldb,
             &vl,
             &vu,
             &il,
@@ -221,36 +240,38 @@ int DiagoLapack<T>::zhegvx_once(const int ncol,
             &M,
             ekb,
             wfc_2d.get_pointer(),
-            &ncol,
+            &ldz,
             work.data(),
             &lwork,
             rwork.data(),
             iwork.data(),
             ifail.data(),
             &info);
-
-    if (info)
+    if (info) {
         throw std::runtime_error("info=" + ModuleBase::GlobalFunc::TO_STRING(info) + ". "
                                  + std::string(__FILE__) + " line "
                                  + std::to_string(__LINE__));
+}
 
-    //	GlobalV::ofs_running<<"lwork="<<work[0]<<"\t"<<"lrwork="<<rwork[0]<<"\t"<<"liwork="<<iwork[0]<<std::endl;
-
-    //lwork = work[0].real();
-    //work.resize(lwork, 0);
-    //int maxlrwork = std::max(lrwork, 3);
-    //rwork.resize(maxlrwork, 0);
-    //iwork.resize(iwork[0], 0);
+    // Query returned optimal lwork in work[0]
+    lwork = static_cast<int>(work[0].real());
+    work.resize(std::max(lwork, 1), 0);
+    // rwork: choose conservative size. Use baseline 7*N plus degeneracy margin
+    lrwork = std::max(3, 7 * n + this->degeneracy_max * n);
+    rwork.resize(lrwork, 0);
+    // LAPACK integer workspace: use conservative size (5*N)
+    liwork = std::max(1, 5 * n);
+    iwork.resize(liwork, 0);
 
-    zhegvx_(&itype,
+        zhegvx_(&itype,
             &jobz,
             &range,
             &uplo,
-            &PARAM.globalv.nlocal,
+            &n,
             h_tmp.c,
-            &PARAM.globalv.nlocal,
+            &lda,
             s_tmp.c,
-            &PARAM.globalv.nlocal,
+            &ldb,
             &vl,
             &vu,
             &il,
@@ -259,52 +280,32 @@ int DiagoLapack<T>::zhegvx_once(const int ncol,
             &M,
             ekb,
             wfc_2d.get_pointer(),
-            &ncol,
+            &ldz,
             work.data(),
             &lwork,
             rwork.data(),
             iwork.data(),
             ifail.data(),
             &info);
-
-    */
-
-        // Use zhegvx_ to compute only the requested eigenvalues (range 'I').
-        // Allocate a reasonably large complex workspace and integer workspace.
-        int lwork2 = std::max(1, ncol * ncol);
-        std::complex<double> *workc = new std::complex<double>[lwork2];
-        int liwork = std::max(1, 5 * ncol);
-        std::vector<int> iwork_large(liwork, 0);
-
-        zhegvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &PARAM.globalv.nlocal,
-            h_tmp.c,
-            &ncol,
-            s_tmp.c,
-            &ncol,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ncol,
-            workc,
-            &lwork2,
-            rwork,
-            iwork_large.data(),
-            ifail.data(),
-            &info);
-
-        delete[] workc;
-        delete[] rwork;
-
-        return info;
+    //	GlobalV::ofs_running<<"M="<<M<<"\t"<<"NZ="<<NZ<<std::endl;
+
+    if (info == 0) {
+        return std::make_pair(info, std::vector<int>{});
+    } else if (info < 0) {
+        return std::make_pair(info, std::vector<int>{});
+    } else if (info % 2) {
+        return std::make_pair(info, ifail);
+    } else if (info / 2 % 2) {
+        return std::make_pair(info, iclustr);
+    } else if (info / 4 % 2) {
+        return std::make_pair(info, std::vector<int>{M, NZ});
+    } else if (info / 16 % 2) {
+        return std::make_pair(info, ifail);
+    } else {
+        throw std::runtime_error("info = " + ModuleBase::GlobalFunc::TO_STRING(info) + ".\n"
+                                 + std::string(__FILE__) + " line "
+                                 + std::to_string(__LINE__));
+    }
 }
 
 template <typename T>
@@ -317,9 +318,9 @@ void DiagoLapack<T>::dsygvx_diag(const int ncol,
 {
     while (true)
     {
-
-        int info_result = dsygvx_once(ncol, nrow, h_mat, s_mat, ekb, wfc_2d);
-        if (info_result == 0) {
+        const std::pair<int, std::vector<int>> info_vec = dsygvx_once(ncol, nrow, h_mat, s_mat, ekb, wfc_2d);
+        post_processing(info_vec.first, info_vec.second);
+        if (info_vec.first == 0) {
             break;
         }
     }
@@ -335,8 +336,9 @@ void DiagoLapack<T>::zhegvx_diag(const int ncol,
 {
     while (true)
     {
-        int info_result = zhegvx_once(ncol, nrow, h_mat, s_mat, ekb, wfc_2d);
-        if (info_result == 0) {
+        const std::pair<int, std::vector<int>> info_vec = zhegvx_once(ncol, nrow, h_mat, s_mat, ekb, wfc_2d);
+        post_processing(info_vec.first, info_vec.second);
+        if (info_vec.first == 0) {
             break;
         }
     }
@@ -354,5 +356,60 @@ void DiagoLapack<T>::post_processing(const int info, const std::vector<int>& vec
     {
         return;
     }
+    else if (info < 0)
+    {
+        const int info_negative = -info;
+        const std::string str_index
+            = (info_negative > 100)
+                  ? ModuleBase::GlobalFunc::TO_STRING(info_negative / 100) + "-th argument "
+                        + ModuleBase::GlobalFunc::TO_STRING(info_negative % 100) + "-entry is illegal.\n"
+                  : ModuleBase::GlobalFunc::TO_STRING(info_negative) + "-th argument is illegal.\n";
+        throw std::runtime_error(str_info_FILE + str_index);
+    }
+    else if (info % 2)
+    {
+        std::string str_ifail = "ifail = ";
+        for (const int i: vec) {
+            str_ifail += ModuleBase::GlobalFunc::TO_STRING(i) + " ";
+}
+        throw std::runtime_error(str_info_FILE + str_ifail);
+    }
+    else if (info / 2 % 2)
+    {
+        int degeneracy_need = 0;
+        for (int irank = 0; irank < GlobalV::DSIZE; ++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";
+        const std::string str_saved
+            = "degeneracy_saved = " + ModuleBase::GlobalFunc::TO_STRING(this->degeneracy_max) + ".\n";
+        if (degeneracy_need <= this->degeneracy_max)
+        {
+            throw std::runtime_error(str_info_FILE + str_need + str_saved);
+        }
+        else
+        {
+            GlobalV::ofs_running << str_need << str_saved;
+            this->degeneracy_max = degeneracy_need;
+            return;
+        }
+    }
+    else if (info / 4 % 2)
+    {
+        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";
+        throw std::runtime_error(str_info_FILE + str_M + str_NZ + str_NBANDS);
+    }
+    else if (info / 16 % 2)
+    {
+        const std::string str_npos = "not positive definite = " + ModuleBase::GlobalFunc::TO_STRING(vec[0]) + ".\n";
+        throw std::runtime_error(str_info_FILE + str_npos);
+    }
+    else
+    {
+        throw std::runtime_error(str_info_FILE);
+    }
 }
 } // namespace hsolver
\ No newline at end of file
diff --git a/source/source_hsolver/diago_lapack.h b/source/source_hsolver/diago_lapack.h
index 53b710ae63..bfdf78ac34 100644
--- a/source/source_hsolver/diago_lapack.h
+++ b/source/source_hsolver/diago_lapack.h
@@ -27,6 +27,10 @@ class DiagoLapack
 
   public:
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
+  #ifdef __MPI
+    // diagnolization used in parallel-k case
+    void diag_pool(hamilt::MatrixBlock<T>& h_mat, hamilt::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
+#endif
 
     void dsygvx_diag(const int ncol,
                      const int nrow,
@@ -41,18 +45,18 @@ class DiagoLapack
                      double* const ekb,
                      psi::Psi<std::complex<double>>& wfc_2d);
 
-    int dsygvx_once(const int ncol,
-                    const int nrow,
-                    const double* const h_mat,
-                    const double* const s_mat,
-                    double* const ekb,
-                    psi::Psi<double>& wfc_2d) const;
-    int zhegvx_once(const int ncol,
-                    const int nrow,
-                    const std::complex<double>* const h_mat,
-                    const std::complex<double>* const s_mat,
-                    double* const ekb,
-                    psi::Psi<std::complex<double>>& wfc_2d) const;
+    std::pair<int, std::vector<int>> dsygvx_once(const int ncol,
+            const int nrow,
+            const double* const h_mat,
+            const double* const s_mat,
+            double* const ekb,
+            psi::Psi<double>& wfc_2d) const;
+    std::pair<int, std::vector<int>> zhegvx_once(const int ncol,
+            const int nrow,
+            const std::complex<double>* const h_mat,
+            const std::complex<double>* const s_mat,
+            double* const ekb,
+            psi::Psi<std::complex<double>>& wfc_2d) const;
 
     int degeneracy_max = 12; // For reorthogonalized memory. 12 followes siesta.
 
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
index 1c4c5a5d61..413c86777a 100644
--- a/source/source_hsolver/hsolver_lcao.cpp
+++ b/source/source_hsolver/hsolver_lcao.cpp
@@ -3,10 +3,10 @@
 #ifdef __MPI
 #include "diago_scalapack.h"
 #include "source_base/module_external/scalapack_connector.h"
-#else
-#include "diago_lapack.h"
 #endif
 
+#include "diago_lapack.h"
+
 #ifdef __CUSOLVERMP
 #include "diago_cusolvermp.h"
 #endif
@@ -174,13 +174,11 @@ void HSolverLCAO<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>&
     }
 #endif
 #endif
-#ifndef __MPI
     else if (this->method == "lapack") // only for single core
     {
         DiagoLapack<T> la;
         la.diag(hm, psi, eigenvalue);
     }
-#endif
     else
     {
         ModuleBase::WARNING_QUIT("HSolverLCAO::solve", "This method is not supported for lcao basis in ABACUS!");
diff --git a/tests/performance/P100_si16_lcao/INPUT b/tests/performance/P100_si16_lcao/INPUT
index 824e290b9a..7ef394405d 100644
--- a/tests/performance/P100_si16_lcao/INPUT
+++ b/tests/performance/P100_si16_lcao/INPUT
@@ -14,6 +14,7 @@ cal_force 1
 cal_stress 1
 #Parameters (3.Basis)
 basis_type		lcao
+ks_solver lapack
 
 #Parameters (4.Smearing)
 smearing_method		gauss

From 89ab2f9f6a11128c6475b4e5c47c7fe62462d243 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 13:19:36 +0800
Subject: [PATCH 03/10] Remove changes on INPU

---
 tests/performance/P100_si16_lcao/INPUT | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/performance/P100_si16_lcao/INPUT b/tests/performance/P100_si16_lcao/INPUT
index 7ef394405d..824e290b9a 100644
--- a/tests/performance/P100_si16_lcao/INPUT
+++ b/tests/performance/P100_si16_lcao/INPUT
@@ -14,7 +14,6 @@ cal_force 1
 cal_stress 1
 #Parameters (3.Basis)
 basis_type		lcao
-ks_solver lapack
 
 #Parameters (4.Smearing)
 smearing_method		gauss

From 29b3cb48d9c4b361c932bd873520c838ce2ddb6d Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 13:30:52 +0800
Subject: [PATCH 04/10] Clean code format

---
 source/source_hsolver/diago_lapack.cpp | 216 ++++++++++++-------------
 1 file changed, 108 insertions(+), 108 deletions(-)

diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp
index f78298a488..a0271da244 100644
--- a/source/source_hsolver/diago_lapack.cpp
+++ b/source/source_hsolver/diago_lapack.cpp
@@ -105,39 +105,39 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> 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 lda = n, ldb = n, ldz = n;
-        dsygvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &n,
-            h_tmp.c,
-            &lda,
-            s_tmp.c,
-            &ldb,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ldz,
-            work.data(),
-            &lwork,
-            iwork.data(),
-            ifail.data(),
-            &info);
+    // 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 lda = n, ldb = n, ldz = n;
+    dsygvx_(&itype,
+        &jobz,
+        &range,
+        &uplo,
+        &n,
+        h_tmp.c,
+        &lda,
+        s_tmp.c,
+        &ldb,
+        &vl,
+        &vu,
+        &il,
+        &iu,
+        &abstol,
+        &M,
+        ekb,
+        wfc_2d.get_pointer(),
+        &ldz,
+        work.data(),
+        &lwork,
+        iwork.data(),
+        ifail.data(),
+        &info);
     if (info) {
         throw std::runtime_error("info = " + ModuleBase::GlobalFunc::TO_STRING(info) + ".\n"
                                  + std::string(__FILE__) + " line "
                                  + std::to_string(__LINE__));
-}
+    }
 
     // Query returned optimal lwork in work[0]
     lwork = static_cast<int>(work[0]);
@@ -146,29 +146,29 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
     liwork = std::max(1, 5 * n);
     iwork.resize(liwork, 0);
 
-        dsygvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &n,
-            h_tmp.c,
-            &lda,
-            s_tmp.c,
-            &ldb,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ldz,
-            work.data(),
-            &lwork,
-            iwork.data(),
-            ifail.data(),
-            &info);
+    dsygvx_(&itype,
+        &jobz,
+        &range,
+        &uplo,
+        &n,
+        h_tmp.c,
+        &lda,
+        s_tmp.c,
+        &ldb,
+        &vl,
+        &vu,
+        &il,
+        &iu,
+        &abstol,
+        &M,
+        ekb,
+        wfc_2d.get_pointer(),
+        &ldz,
+        work.data(),
+        &lwork,
+        iwork.data(),
+        ifail.data(),
+        &info);
     //	GlobalV::ofs_running<<"M="<<M<<"\t"<<"NZ="<<NZ<<std::endl;
 
     if (info == 0) {
@@ -218,40 +218,40 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> 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 lda = n, ldb = n, ldz = n;
-        zhegvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &n,
-            h_tmp.c,
-            &lda,
-            s_tmp.c,
-            &ldb,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ldz,
-            work.data(),
-            &lwork,
-            rwork.data(),
-            iwork.data(),
-            ifail.data(),
-            &info);
+    // 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 lda = n, ldb = n, ldz = n;
+    zhegvx_(&itype,
+        &jobz,
+        &range,
+        &uplo,
+        &n,
+        h_tmp.c,
+        &lda,
+        s_tmp.c,
+        &ldb,
+        &vl,
+        &vu,
+        &il,
+        &iu,
+        &abstol,
+        &M,
+        ekb,
+        wfc_2d.get_pointer(),
+        &ldz,
+        work.data(),
+        &lwork,
+        rwork.data(),
+        iwork.data(),
+        ifail.data(),
+        &info);
     if (info) {
         throw std::runtime_error("info=" + ModuleBase::GlobalFunc::TO_STRING(info) + ". "
                                  + std::string(__FILE__) + " line "
                                  + std::to_string(__LINE__));
-}
+    }
 
     // Query returned optimal lwork in work[0]
     lwork = static_cast<int>(work[0].real());
@@ -263,30 +263,30 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     liwork = std::max(1, 5 * n);
     iwork.resize(liwork, 0);
 
-        zhegvx_(&itype,
-            &jobz,
-            &range,
-            &uplo,
-            &n,
-            h_tmp.c,
-            &lda,
-            s_tmp.c,
-            &ldb,
-            &vl,
-            &vu,
-            &il,
-            &iu,
-            &abstol,
-            &M,
-            ekb,
-            wfc_2d.get_pointer(),
-            &ldz,
-            work.data(),
-            &lwork,
-            rwork.data(),
-            iwork.data(),
-            ifail.data(),
-            &info);
+    zhegvx_(&itype,
+        &jobz,
+        &range,
+        &uplo,
+        &n,
+        h_tmp.c,
+        &lda,
+        s_tmp.c,
+        &ldb,
+        &vl,
+        &vu,
+        &il,
+        &iu,
+        &abstol,
+        &M,
+        ekb,
+        wfc_2d.get_pointer(),
+        &ldz,
+        work.data(),
+        &lwork,
+        rwork.data(),
+        iwork.data(),
+        ifail.data(),
+        &info);
     //	GlobalV::ofs_running<<"M="<<M<<"\t"<<"NZ="<<NZ<<std::endl;
 
     if (info == 0) {
@@ -371,7 +371,7 @@ void DiagoLapack<T>::post_processing(const int info, const std::vector<int>& vec
         std::string str_ifail = "ifail = ";
         for (const int i: vec) {
             str_ifail += ModuleBase::GlobalFunc::TO_STRING(i) + " ";
-}
+        }
         throw std::runtime_error(str_info_FILE + str_ifail);
     }
     else if (info / 2 % 2)
@@ -379,7 +379,7 @@ void DiagoLapack<T>::post_processing(const int info, const std::vector<int>& vec
         int degeneracy_need = 0;
         for (int irank = 0; irank < GlobalV::DSIZE; ++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";
         const std::string str_saved
             = "degeneracy_saved = " + ModuleBase::GlobalFunc::TO_STRING(this->degeneracy_max) + ".\n";

From 8ec73c0e9a8421ab767d1ecdd9197a8195656df7 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 13:40:09 +0800
Subject: [PATCH 05/10] Add unit tiest for lapack solver

---
 source/source_hsolver/test/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt
index 72ad05e0ba..31c41752ed 100644
--- a/source/source_hsolver/test/CMakeLists.txt
+++ b/source/source_hsolver/test/CMakeLists.txt
@@ -119,7 +119,7 @@ if (ENABLE_MPI)
     ../kernels/cuda/diag_cusolver.cu
   )
   endif()
-else()
+#else()
   if(ENABLE_LCAO)
   AddTest(
         TARGET MODULE_HSOLVER_Lapack

From 7d60e3b1db49960e5ed760d04cf401796f46dc20 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 16:40:59 +0800
Subject: [PATCH 06/10] remove original diago_lapack_test and add lapack test
 to diago_lcao_test

---
 source/source_hsolver/test/CMakeLists.txt      | 8 --------
 source/source_hsolver/test/diago_lcao_test.cpp | 8 ++++++++
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt
index 31c41752ed..c1fd32b832 100644
--- a/source/source_hsolver/test/CMakeLists.txt
+++ b/source/source_hsolver/test/CMakeLists.txt
@@ -119,14 +119,6 @@ if (ENABLE_MPI)
     ../kernels/cuda/diag_cusolver.cu
   )
   endif()
-#else()
-  if(ENABLE_LCAO)
-  AddTest(
-        TARGET MODULE_HSOLVER_Lapack
-        LIBS parameter  ${math_libs} base psi device
-        SOURCES diago_lapack_test.cpp ../diago_lapack.cpp 
-      )
-  endif()
 endif()
 install(FILES H-KPoints-Si2.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 install(FILES H-GammaOnly-Si2.dat DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
index bbfaad4a1e..c58e514c26 100644
--- a/source/source_hsolver/test/diago_lcao_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -1,4 +1,5 @@
 #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"
@@ -226,6 +227,11 @@ class DiagoPrepare
                 hsolver::DiagoScalapack<T> dh;
                 dh.diag(&hmtest, psi, e_solver.data());
             }
+            else if (ks_solver == "lapack")
+            {
+                hsolver::DiagoLapack<T> la;
+                la.diag(&hmtest, psi, e_solver.data());
+            }
     #ifdef __ELPA
             else if (ks_solver == "genelpa")
             {
@@ -317,6 +323,8 @@ INSTANTIATE_TEST_SUITE_P(
 #endif
         DiagoPrepare<double>(0, 0, 1, 0, "scalapack_gvx", "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat"),
         DiagoPrepare<double>(0, 0, 32, 0, "scalapack_gvx", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat")));
+        DiagoPrepare<double>(0, 0, 1, 0, "lapack", "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat"),
+        DiagoPrepare<double>(0, 0, 32, 0, "lapack", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat")));
 
 class DiagoKPointsTest : public ::testing::TestWithParam<DiagoPrepare<std::complex<double>>>
 {

From 8c447d53c655b8d40c3ff83a9c5d358b0b62ad29 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 16:43:50 +0800
Subject: [PATCH 07/10] Fix a problem

---
 source/source_hsolver/test/diago_lcao_test.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
index c58e514c26..2668748975 100644
--- a/source/source_hsolver/test/diago_lcao_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -322,7 +322,7 @@ INSTANTIATE_TEST_SUITE_P(
         DiagoPrepare<double>(0, 0, 32, 0, "genelpa", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat"),
 #endif
         DiagoPrepare<double>(0, 0, 1, 0, "scalapack_gvx", "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat"),
-        DiagoPrepare<double>(0, 0, 32, 0, "scalapack_gvx", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat")));
+        DiagoPrepare<double>(0, 0, 32, 0, "scalapack_gvx", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat"),
         DiagoPrepare<double>(0, 0, 1, 0, "lapack", "H-GammaOnly-Si2.dat", "S-GammaOnly-Si2.dat"),
         DiagoPrepare<double>(0, 0, 32, 0, "lapack", "H-GammaOnly-Si64.dat", "S-GammaOnly-Si64.dat")));
 

From 00ccc42dcc1522f84d16dc8f2c2d77dab9365493 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 16:46:43 +0800
Subject: [PATCH 08/10] Fix Cmake error

---
 source/source_hsolver/test/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt
index c1fd32b832..41db30fa3e 100644
--- a/source/source_hsolver/test/CMakeLists.txt
+++ b/source/source_hsolver/test/CMakeLists.txt
@@ -92,7 +92,7 @@ if (ENABLE_MPI)
     AddTest(
       TARGET MODULE_HSOLVER_LCAO
       LIBS parameter  ${math_libs} ELPA::ELPA base genelpa psi device
-      SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_scalapack.cpp 
+      SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_scalapack.cpp ../diago_lapack.cpp
     )
     else()
       AddTest(

From e27d3b5e38b6a6f376d96290a173c4d3cac10215 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 17:00:05 +0800
Subject: [PATCH 09/10] Fix another cmake bug

---
 source/source_hsolver/test/CMakeLists.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt
index 41db30fa3e..217f8251b3 100644
--- a/source/source_hsolver/test/CMakeLists.txt
+++ b/source/source_hsolver/test/CMakeLists.txt
@@ -98,7 +98,7 @@ if (ENABLE_MPI)
       AddTest(
         TARGET MODULE_HSOLVER_LCAO
         LIBS parameter  ${math_libs} base psi device
-        SOURCES diago_lcao_test.cpp ../diago_scalapack.cpp 
+        SOURCES diago_lcao_test.cpp ../diago_scalapack.cpp ../diago_lapack.cpp
       )
     endif()
 

From 4eadfe3c6d3a0915b6cc75fa6d245b5315d1c201 Mon Sep 17 00:00:00 2001
From: Critsium-xy <tsfxwbbzxy@163.com>
Date: Fri, 9 Jan 2026 17:01:49 +0800
Subject: [PATCH 10/10] Final fix cmake bug

---
 source/source_hsolver/test/diago_lcao_test.cpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
index 2668748975..8952a9d7c6 100644
--- a/source/source_hsolver/test/diago_lcao_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -75,6 +75,8 @@ class DiagoPrepare
         if (ks_solver == "scalapack_gvx")
             ;
 //             dh = new hsolver::DiagoScalapack<T>;
+        else if (ks_solver == "lapack")
+            ;
 #ifdef __ELPA
         else if (ks_solver == "genelpa")
             ;