Harden PPCG generalized eigensolver fallback

Author: cheerly-pku

source/source_hsolver/ppcg/diago_ppcg_lapack.hpp

@@ -48,6 +48,17 @@ struct HermitianLapack
     {
         std::vector<Scalar> eigvec(n * n, Scalar(0));
         container::kernels::lapack_hegvd<Scalar, Device>()(n, n, a, b, w, eigvec.data());
+        for (int j = 0; j < n; ++j)
+        {
+            if (!std::isfinite(w[j]))
+                throw std::runtime_error("PPCG: hegvd returned non-finite eigenvalue.");
+
+            Real nrm2 = 0;
+            for (int i = 0; i < n; ++i)
+                nrm2 += static_cast<Real>(std::norm(eigvec[i + j * n]));
+            if (nrm2 <= static_cast<Real>(1e-30))
+                throw std::runtime_error("PPCG: hegvd returned a zero eigenvector.");
+        }
         std::copy(eigvec.begin(), eigvec.end(), a);
     }
 

source/source_hsolver/ppcg/diago_ppcg_subspace.hpp

@@ -168,8 +168,12 @@ void DiagoPPCG<T, Device>::solve_small_generalized(
     // All attempts failed — set eigenvectors to identity (no update).
     std::fill(subspace.k.begin(), subspace.k.end(), T(0));
     for (int i = 0; i < dim; ++i)
+    {
         subspace.k[i + i * dim] = T(1);
-    std::fill(subspace.eval.begin(), subspace.eval.end(), static_cast<Real>(0));
+        subspace.eval[i] = static_cast<Real>(std::real(k0[i + i * dim]))
+                         / std::max(static_cast<Real>(std::real(m0[i + i * dim])),
+                                    static_cast<Real>(1e-30));
+    }
 }
 
 // ---------------------------------------------------------------------------