mohanchen · Sep 14, 2026 · Sep 7, 2026 · Sep 6, 2026 · Sep 6, 2026 · Sep 7, 2026
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -3419,6 +3419,7 @@
 - **Availability**: *[`symmetry`](#symmetry)==1 and ([`dft_functional`](#dft_functional) in [hse, hf, pbe0, scan0] or ([`basis_type`](#basis_type)==lcao and [`rpa`](#rpa)==true))*
 - **Description**: - False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization
   - True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation
+  For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting.
 - **Default**: True
 
 ### out_ri_cv

diff --git a/docs/parameters.yaml b/docs/parameters.yaml
@@ -4681,6 +4681,7 @@ parameters:
     description: |
       * False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization
       * True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation
+      For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting.
     default_value: "True"
     unit: ""
     availability: "symmetry==1 and (dft_functional in [hse, hf, pbe0, scan0] or (basis_type==lcao and rpa==true))"

diff --git a/source/source_hsolver/diago_elpa_native.cpp b/source/source_hsolver/diago_elpa_native.cpp
@@ -3,6 +3,7 @@
 #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/module_external/scalapack_connector.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_hsolver/module_genelpa/elpa_new.h"
@@ -71,6 +72,79 @@ void DiagoElpaNative<T>::diag_pool(ModuleBase::MatrixBlock<T>& h_mat,
     std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(narows * nacols);
 
+    // The complex LCAO matrices follow the LAPACK UPLO='U' convention: only
+    // their upper triangles are guaranteed to contain the Hermitian matrix.
+    // ELPA's native generalized solver consumes both triangles, so complete
+    // private Hermitian copies before the solve.
+    std::vector<T> h_work;
+    std::vector<T> s_work;
+    T* h_elpa = h_mat.p;
+    T* s_elpa = s_mat.p;
+    int decomposed_state = this->DecomposedState;
+    if (!std::is_same<T, double>::value)
+    {
+        h_work.resize(narows * nacols);
+        s_work.resize(narows * nacols);
+        const int one = 1;
+        ScalapackConnector::tranc(nFull,
+                                  nFull,
+                                  T(1.0),
+                                  h_mat.p,
+                                  one,
+                                  one,
+                                  h_mat.desc,
+                                  T(0.0),
+                                  h_work.data(),
+                                  one,
+                                  one,
+                                  h_mat.desc);
+        ScalapackConnector::tranc(nFull,
+                                  nFull,
+                                  T(1.0),
+                                  s_mat.p,
+                                  one,
+                                  one,
+                                  s_mat.desc,
+                                  T(0.0),
+                                  s_work.data(),
+                                  one,
+                                  one,
+                                  s_mat.desc);
+        const auto local_to_global = [](const int local_index,
+                                        const int block_size,
+                                        const int process_coordinate,
+                                        const int source_coordinate,
+                                        const int process_count) {
+            if (source_coordinate < 0)
+            {
+                return local_index;
+            }
+            const int process_offset
+                = (process_coordinate - source_coordinate + process_count) % process_count;
+            return ((local_index / block_size) * process_count + process_offset) * block_size
+                   + local_index % block_size;
+        };
+        for (int local_col = 0; local_col < nacols; ++local_col)
+        {
+            const int global_col
+                = local_to_global(local_col, h_mat.desc[5], mypcol, h_mat.desc[7], npcols);
+            for (int local_row = 0; local_row < narows; ++local_row)
+            {
+                const int global_row
+                    = local_to_global(local_row, h_mat.desc[4], myprow, h_mat.desc[6], nprows);
+                if (global_row <= global_col)
+                {
+                    const int local_index = local_row + local_col * narows;
+                    h_work[local_index] = h_mat.p[local_index];
+                    s_work[local_index] = s_mat.p[local_index];
+                }
+            }
+        }
+        h_elpa = h_work.data();
+        s_elpa = s_work.data();
+        decomposed_state = 0;
+    }
+
     if (elpa_init(20210430) != ELPA_OK)
     {
         fprintf(stderr, "Error: ELPA API version not supported");
@@ -114,11 +188,11 @@ void DiagoElpaNative<T>::diag_pool(ModuleBase::MatrixBlock<T>& h_mat,
 #endif
 
     elpa_generalized_eigenvectors(handle,
-                                  h_mat.p,
-                                  s_mat.p,
+                                  h_elpa,
+                                  s_elpa,
                                   eigen.data(),
                                   eigenvectors.data(),
-                                  this->DecomposedState,
+                                  decomposed_state,
                                   &success);
     elpa_deallocate(handle, &success);
     elpa_uninit(&success);

diff --git a/source/source_hsolver/kernels/cuda/diag_cusolver.cu b/source/source_hsolver/kernels/cuda/diag_cusolver.cu
@@ -9,7 +9,8 @@ Diag_Cusolver_gvd::Diag_Cusolver_gvd(){
 
     itype = CUSOLVER_EIG_TYPE_1; // A*x = (lambda)*B*x
     jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors.
-    uplo = CUBLAS_FILL_MODE_LOWER;
+    // LCAO supplies the authoritative upper triangle of the Hermitian H/S matrices.
+    uplo = CUBLAS_FILL_MODE_UPPER;
 
     d_A = NULL;
     d_B = NULL;

diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt
@@ -93,7 +93,7 @@ if (ENABLE_MPI)
     AddTest(
       TARGET MODULE_HSOLVER_LCAO
       LIBS parameter  ELPA::ELPA base genelpa psi device
-      SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_scalapack.cpp ../diago_lapack.cpp
+      SOURCES diago_lcao_test.cpp ../diago_elpa.cpp ../diago_elpa_native.cpp ../diago_scalapack.cpp ../diago_lapack.cpp
     )
     else()
       AddTest(

diff --git a/source/source_hsolver/test/diago_elpa_utils.h b/source/source_hsolver/test/diago_elpa_utils.h
@@ -167,7 +167,7 @@ void lapack_diago(double *hmatrix, double *smatrix, double *e, int &nFull)
     const char jobz = 'V'; // 'N':only calc eigenvalue, 'V': eigenvalues and eigenvectors
     const char uplo = 'U'; // Upper triangles
     int lwork = (nFull + 2) * nFull, info = 0;
-    double *ev = new double[nFull * nFull];
+    double *ev = new double[lwork];
 
     double *a = new double[nFull * nFull];
     double *b = new double[nFull * nFull];
@@ -196,7 +196,7 @@ void lapack_diago(std::complex<double> *hmatrix, std::complex<double> *smatrix,
     const char uplo = 'U'; // Upper triangles
     int lwork = (nFull + 1) * nFull, info = 0;
     double *rwork = new double[3 * nFull - 2];
-    std::complex<double> *ev = new std::complex<double>[nFull * nFull];
+    std::complex<double> *ev = new std::complex<double>[lwork];
 
     std::complex<double> *a = new std::complex<double>[nFull * nFull];
     std::complex<double> *b = new std::complex<double>[nFull * nFull];

diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
@@ -134,6 +134,21 @@ class DiagoPrepare
         return ok;
     }
 
+    void poison_lower_triangle()
+    {
+        // The distributed solver buffers are column-major. Keep the original
+        // row-major fixtures intact for the independent LAPACK reference.
+        for (int col = 0; col < nlocal; ++col)
+        {
+            for (int row = col + 1; row < nlocal; ++row)
+            {
+                const int index = row + col * nlocal;
+                this->h_local[index] = T(123.0 + row + col);
+                this->s_local[index] = T(0.0);
+            }
+        }
+    }
+
     void print_hs()
     {
         if (!PRINT_HS)
@@ -203,6 +218,10 @@ class DiagoPrepare
     {
         this->pb2d();
         this->distribute_data();
+        if (ks_solver == "cusolver")
+        {
+            this->poison_lower_triangle();
+        }
         this->print_hs();
         this->set_env();
 

diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -9,6 +9,7 @@
 #include <vector>
 #ifdef __ELPA
 #include "source_hsolver/diago_elpa.h"
+#include "source_hsolver/diago_elpa_native.h"
 #include "source_hsolver/module_genelpa/elpa_solver.h"
 #endif
 #include "source_base/module_external/scalapack_connector.h"
@@ -136,6 +137,47 @@ class DiagoPrepare
         return ok;
     }
 
+#ifdef __ELPA
+    void diago_native_with_poisoned_lower()
+    {
+        this->pb2d();
+        this->distribute_data();
+        int nprows;
+        int npcols;
+        int myprow;
+        int mypcol;
+        Cblacs_gridinfo(icontxt, &nprows, &npcols, &myprow, &mypcol);
+        for (int col = 0; col < hmtest.ncol; ++col)
+        {
+            const int global_col = (col / nb2d * npcols + mypcol) * nb2d + col % nb2d;
+            for (int row = 0; row < hmtest.nrow; ++row)
+            {
+                const int global_row = (row / nb2d * nprows + myprow) * nb2d + row % nb2d;
+                if (global_row > global_col)
+                {
+                    const int index = row + col * hmtest.nrow;
+                    h_local[index] = T(123.0 + global_row + global_col);
+                    s_local[index] = T(0.0);
+                }
+            }
+        }
+        hmtest.h_local = h_local;
+        hmtest.s_local = s_local;
+        ModuleBase::MatrixBlock<T> h_mat;
+        ModuleBase::MatrixBlock<T> s_mat;
+        hmtest.matrix(h_mat, s_mat);
+        hsolver::DiagoElpaNative<T> solver(nlocal, nbands, false);
+        solver.diag(h_mat, s_mat, psi, e_solver.data());
+        EXPECT_EQ(hmtest.h_local, h_local);
+        EXPECT_EQ(hmtest.s_local, s_local);
+        // A second solve must not reuse an in-place decomposition of the
+        // private overlap copy or alter the caller's matrix storage.
+        solver.diag(h_mat, s_mat, psi, e_solver.data());
+        EXPECT_EQ(hmtest.h_local, h_local);
+        EXPECT_EQ(hmtest.s_local, s_local);
+    }
+#endif
+
     void print_hs()
     {
         if (!PRINT_HS)
@@ -371,6 +413,26 @@ INSTANTIATE_TEST_SUITE_P(
         DiagoPrepare<std::complex<double>>(0, 0, 32, 0, "scalapack_gvx", "H-KPoints-Si64.dat", "S-KPoints-Si64.dat")));
 
 #ifdef __ELPA
+class DiagoElpaNativeUpperTest : public ::testing::TestWithParam<int>
+{
+};
+
+TEST_P(DiagoElpaNativeUpperTest, PreservesInputsAndIgnoresLowerTriangle)
+{
+    DiagoPrepare<std::complex<double>> dp(0, 0, GetParam(), 0, "genelpa",
+                                          "H-KPoints-Si2.dat", "S-KPoints-Si2.dat");
+    ASSERT_TRUE(dp.produce_HS());
+    dp.diago_native_with_poisoned_lower();
+    if (dp.myrank == 0)
+    {
+        dp.diago_lapack();
+        std::stringstream out_info;
+        EXPECT_TRUE(dp.compare_eigen(out_info)) << out_info.str();
+    }
+}
+
+INSTANTIATE_TEST_SUITE_P(BlockSizes, DiagoElpaNativeUpperTest, ::testing::Values(1, 2, 3));
+
 TEST(DiagoElpaComplexTest, UsesAuthoritativeUpperTriangle)
 {
     std::stringstream out_info;

diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp
@@ -525,7 +525,8 @@ void ReadInput::item_exx()
         item.category = "Exact Exchange (LCAO)";
         item.type = "Boolean";
         item.description = R"(* False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization
-* True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation)";
+* True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation
+For multi-k calculations, D(k) is averaged over the unitary little group of each irreducible k point before star expansion, for either setting.)";
         item.default_value = "True";
         item.unit = "";
         item.set_availability("symmetry==1 and (dft_functional in [hse, hf, pbe0, scan0] or (basis_type==lcao and rpa==true))");