fix: use real-symmetric fast path in solve_hermitian_gevp

When the Hermitian Gram matrices are essentially real (max |Im| < 1e-14
* max |Re|), solve the generalized eigenvalue problem directly with
dsygvd on the real part.  This avoids the 2n x 2n block-real reduction
for the common case of real wavefunctions / real operators.

Only when the matrices carry significant imaginary parts (complex
wavefunctions, e.g. k-point DFT) does the code fall back to the
full block-real Hermitian solve.

Should resolve the SEGFAULT seen in CI, which was likely caused by
the block-real path when applied to numerically degenerate real
matrices.

Author: cheerly-pku

source/source_hsolver/diago_ppcg.cpp

@@ -420,97 +420,158 @@ void DiagoPPCG<T, Device>::solve_hermitian_gevp(
     std::vector<T>& k_herm, std::vector<T>& m_herm,
     int n, std::vector<Real>& eval, std::vector<T>& evec) const
 {
-    const int n2 = 2 * n;
-    std::vector<Real> k_blk(n2 * n2, static_cast<Real>(0));
-    std::vector<Real> m_blk(n2 * n2, static_cast<Real>(0));
+    // -----------------------------------------------------------------------
+    // Determine whether the matrices carry a significant imaginary part.
+    // When wavefunctions are real (or the operator is real and the basis
+    // has been kept real) the Hermitian Gram matrices are strictly real,
+    // and we can use the fast real-symmetric dsygvd path directly.
+    // -----------------------------------------------------------------------
+    const int n2 = n * n;
+    Real max_imag = 0, max_abs = 0;
+    for (int i = 0; i < n2; ++i)
+    {
+        max_imag = std::max(max_imag, std::abs(std::imag(k_herm[i])));
+        max_imag = std::max(max_imag, std::abs(std::imag(m_herm[i])));
+        max_abs  = std::max(max_abs,  std::abs(std::real(k_herm[i])));
+        max_abs  = std::max(max_abs,  std::abs(std::real(m_herm[i])));
+    }
+    const bool is_real = (max_imag <= max_abs * static_cast<Real>(1e-14));
 
-    // Build the 2n×2n real block matrices.
-    for (int j = 0; j < n; ++j)
+    if (is_real)
     {
-        for (int i = 0; i < n; ++i)
+        // ===== Fast path: real-symmetric matrices =====
+        std::vector<Real> k_real(n2), m_real(n2);
+        for (int i = 0; i < n2; ++i)
         {
-            const T kij = k_herm[i + j * n];
-            const T mij = m_herm[i + j * n];
-            const Real kr = std::real(kij);
-            const Real ki = std::imag(kij);
-            const Real mr = std::real(mij);
-            const Real mi = std::imag(mij);
-
-            // Upper-left  n×n block: Re(K), Re(M)
-            k_blk[i     + j     * n2] = kr;
-            m_blk[i     + j     * n2] = mr;
-            // Upper-right n×n block: -Im(K), -Im(M)
-            k_blk[i     + (j+n) * n2] = -ki;
-            m_blk[i     + (j+n) * n2] = -mi;
-            // Lower-left  n×n block: Im(K),  Im(M)
-            k_blk[i + n + j     * n2] = ki;
-            m_blk[i + n + j     * n2] = mi;
-            // Lower-right n×n block: Re(K), Re(M)
-            k_blk[i + n + (j+n) * n2] = kr;
-            m_blk[i + n + (j+n) * n2] = mr;
+            k_real[i] = std::real(k_herm[i]);
+            m_real[i] = std::real(m_herm[i]);
         }
-    }
 
-    // Symmetrise explicitly.
-    for (int j = 0; j < n2; ++j)
-        for (int i = 0; i < j; ++i)
+        const std::vector<Real> m0 = m_real;
+        const Real shifts[] = {static_cast<Real>(1e-10),
+                                static_cast<Real>(1e-8),
+                                static_cast<Real>(1e-6)};
+        bool solved = false;
+        for (int attempt = 0; attempt < 3; ++attempt)
         {
-            k_blk[j + i * n2] = k_blk[i + j * n2];
-            m_blk[j + i * n2] = m_blk[i + j * n2];
+            try
+            {
+                Lapack<Real>::sygvd(n, k_real.data(), m_real.data(),
+                                    eval.data());
+                solved = true;
+                break;
+            }
+            catch (const std::runtime_error&)
+            {
+                m_real = m0;
+                for (int i = 0; i < n; ++i)
+                    m_real[i + i * n] += shifts[attempt];
+            }
         }
 
-    // M-save/restore for retry (same pattern as solve_small_generalized).
-    const std::vector<Real> m0 = m_blk;
-    std::vector<Real> eval2(n2, static_cast<Real>(0));
-    const Real shifts[] = {static_cast<Real>(1e-10),
-                           static_cast<Real>(1e-8),
-                           static_cast<Real>(1e-6)};
-    bool solved = false;
-    for (int attempt = 0; attempt < 3; ++attempt)
-    {
-        try
+        if (solved)
         {
-            Lapack<Real>::sygvd(n2, k_blk.data(), m_blk.data(), eval2.data());
-            solved = true;
-            break;
+            // k_real now holds the real eigenvectors column-by-column.
+            evec.assign(n2, T(0));
+            for (int j = 0; j < n; ++j)
+                for (int i = 0; i < n; ++i)
+                    evec[i + j * n] = T(k_real[i + j * n], 0);
         }
-        catch (const std::runtime_error&)
+        else
         {
-            m_blk = m0;
-            for (int i = 0; i < n2; ++i)
-                m_blk[i + i * n2] += shifts[attempt];
+            evec.assign(n2, T(0));
+            for (int i = 0; i < n; ++i)
+                evec[i + i * n] = T(1, 0);
+            eval.assign(n, static_cast<Real>(0));
         }
+        return;
     }
 
-    if (solved)
+    // ===== General path: complex-Hermitian matrices =====
+    // Convert to a 2n×2n real block system and solve with dsygvd.
     {
-        // Extract the n lowest eigenvalues (each pair gives two nearly-equal
-        // eigenvalues; we take the first of each pair).
-        eval.assign(n, static_cast<Real>(0));
-        for (int i = 0; i < n; ++i)
-            eval[i] = eval2[2 * i];
+        const int n2blk = 2 * n;
+        std::vector<Real> k_blk(n2blk * n2blk, static_cast<Real>(0));
+        std::vector<Real> m_blk(n2blk * n2blk, static_cast<Real>(0));
 
-        // Extract complex eigenvectors from the block eigenvectors.
-        // k_blk now holds the real eigenvectors column-by-column (n2×n2).
-        evec.assign(n * n, T(0));
         for (int j = 0; j < n; ++j)
         {
-            const int jj = 2 * j;   // use the first of each degenerate pair
             for (int i = 0; i < n; ++i)
             {
-                const Real vr = k_blk[i     + jj * n2];
-                const Real vi = k_blk[i + n + jj * n2];
-                evec[i + j * n] = T(vr, vi);
+                const T kij = k_herm[i + j * n];
+                const T mij = m_herm[i + j * n];
+                const Real kr = std::real(kij);
+                const Real ki = std::imag(kij);
+                const Real mr = std::real(mij);
+                const Real mi = std::imag(mij);
+
+                k_blk[i     + j         * n2blk] = kr;
+                m_blk[i     + j         * n2blk] = mr;
+                k_blk[i     + (j+n)     * n2blk] = -ki;
+                m_blk[i     + (j+n)     * n2blk] = -mi;
+                k_blk[i + n + j         * n2blk] = ki;
+                m_blk[i + n + j         * n2blk] = mi;
+                k_blk[i + n + (j+n)     * n2blk] = kr;
+                m_blk[i + n + (j+n)     * n2blk] = mr;
             }
         }
-    }
-    else
-    {
-        // All attempts failed — set eigenvectors to identity (no update).
-        std::fill(evec.begin(), evec.end(), T(0));
-        for (int i = 0; i < n; ++i)
-            evec[i + i * n] = T(1, 0);
-        std::fill(eval.begin(), eval.end(), static_cast<Real>(0));
+
+        // Symmetrise.
+        for (int j = 0; j < n2blk; ++j)
+            for (int i = 0; i < j; ++i)
+            {
+                k_blk[j + i * n2blk] = k_blk[i + j * n2blk];
+                m_blk[j + i * n2blk] = m_blk[i + j * n2blk];
+            }
+
+        const std::vector<Real> m0 = m_blk;
+        std::vector<Real> eval2(n2blk, static_cast<Real>(0));
+        const Real shifts[] = {static_cast<Real>(1e-10),
+                                static_cast<Real>(1e-8),
+                                static_cast<Real>(1e-6)};
+        bool solved = false;
+        for (int attempt = 0; attempt < 3; ++attempt)
+        {
+            try
+            {
+                Lapack<Real>::sygvd(n2blk, k_blk.data(), m_blk.data(),
+                                    eval2.data());
+                solved = true;
+                break;
+            }
+            catch (const std::runtime_error&)
+            {
+                m_blk = m0;
+                for (int i = 0; i < n2blk; ++i)
+                    m_blk[i + i * n2blk] += shifts[attempt];
+            }
+        }
+
+        if (solved)
+        {
+            eval.assign(n, static_cast<Real>(0));
+            for (int i = 0; i < n; ++i)
+                eval[i] = eval2[2 * i];
+
+            evec.assign(n2, T(0));
+            for (int j = 0; j < n; ++j)
+            {
+                const int jj = 2 * j;
+                for (int i = 0; i < n; ++i)
+                {
+                    const Real vr = k_blk[i     + jj * n2blk];
+                    const Real vi = k_blk[i + n + jj * n2blk];
+                    evec[i + j * n] = T(vr, vi);
+                }
+            }
+        }
+        else
+        {
+            evec.assign(n2, T(0));
+            for (int i = 0; i < n; ++i)
+                evec[i + i * n] = T(1, 0);
+            eval.assign(n, static_cast<Real>(0));
+        }
     }
 }