DGMRES: guard against Arnoldi happy breakdown and singular Hessenberg

libeigen/eigen!2454

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

unsupported/Eigen/src/IterativeSolvers/DGMRES.h

diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
index b251ff6..25f16be 100644
--- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
@@ -316,18 +316,31 @@
       m_H(i, it) = coef;
       m_Hes(i, it) = coef;
     }
-    // Normalize the vector
+    // Normalize the vector. coef == 0 is an Arnoldi happy breakdown: the new
+    // direction lies in span(V[0..it]), so skip the division (which would
+    // poison m_V.col(it+1) with NaN) and fall through to the termination
+    // check below.
     coef = tv1.norm();
-    m_V.col(it + 1) = tv1 / coef;
+    const bool happy_breakdown = numext::is_exactly_zero(coef);
+    if (!happy_breakdown) {
+      m_V.col(it + 1) = tv1 / coef;
+    }
     m_H(it + 1, it) = coef;
     //     m_Hes(it+1,it) = coef;
 
-    // FIXME: Check for happy breakdown.
-
     // Update Hessenberg matrix with Givens rotations
     for (Index i = 1; i <= it; ++i) {
       m_H.col(it).applyOnTheLeft(i - 1, i, gr[i - 1].adjoint());
     }
+
+    // If the rotated diagonal is also zero, the reduced triangular system
+    // becomes singular and the back-substitution below would produce Inf/NaN.
+    // Stop with NumericalIssue instead of polluting x.
+    if (happy_breakdown && numext::is_exactly_zero(m_H(it, it))) {
+      m_info = NumericalIssue;
+      break;
+    }
+
     // Compute the new plane rotation
     gr[it].makeGivens(m_H(it, it), m_H(it + 1, it));
     // Apply the new rotation
@@ -339,15 +352,15 @@
     it++;
     nbIts++;
 
-    if (m_error < m_tolerance) {
-      // The method has converged
+    if (m_error < m_tolerance || happy_breakdown) {
+      // Happy breakdown: residual on the current subspace is exactly zero, so
+      // the it-dim triangular system yields the exact solution.
       m_info = Success;
       break;
     }
   }
 
-  // Compute the new coefficients by solving the least square problem
-  // FIXME: Check first if the matrix is singular (zero diagonal).
+  // Compute the new coefficients by solving the least square problem.
   DenseVector nrs(m_restart);
   nrs = m_H.topLeftCorner(it, it).template triangularView<Upper>().solve(g.head(it));
 

unsupported/test/dgmres.cpp

diff --git a/unsupported/test/dgmres.cpp b/unsupported/test/dgmres.cpp
index 940800c..04a4414 100644
--- a/unsupported/test/dgmres.cpp
+++ b/unsupported/test/dgmres.cpp
@@ -24,7 +24,46 @@
   // CALL_SUBTEST( check_sparse_square_solving(dgmres_colmajor_ssor)     );
 }
 
+// Regression: Arnoldi breakdown used to divide by zero (producing NaN in the
+// Krylov basis) and solve a singular triangular system, silently returning
+// Inf with info() == Success. Exercise both the pathological (rank-deficient
+// pivot) and benign (exact Krylov subspace) breakdown paths.
+template <typename T>
+void test_dgmres_breakdown_T() {
+  typedef SparseMatrix<T> Mat;
+  typedef Matrix<T, 2, 1> Vec;
+
+  // Nilpotent A with singular Hessenberg pivot on the first step.
+  Mat A(2, 2);
+  A.insert(0, 1) = T(1);
+  A.makeCompressed();
+  Vec b;
+  b << T(1), T(0);
+
+  DGMRES<Mat, IdentityPreconditioner> solver;
+  solver.compute(A);
+  Vec x = solver.solve(b);
+  VERIFY(x.allFinite());
+  VERIFY(solver.info() != Success);
+
+  // Diagonal A with b in an eigenspace: Arnoldi converges after one step.
+  Mat D(2, 2);
+  D.insert(0, 0) = T(2);
+  D.insert(1, 1) = T(2);
+  D.makeCompressed();
+  Vec d;
+  d << T(2), T(2);
+
+  DGMRES<Mat, DiagonalPreconditioner<T> > solver2;
+  solver2.compute(D);
+  Vec y = solver2.solve(d);
+  VERIFY_IS_EQUAL(solver2.info(), Success);
+  VERIFY_IS_APPROX(y, (Vec() << T(1), T(1)).finished());
+}
+
 EIGEN_DECLARE_TEST(dgmres) {
   CALL_SUBTEST_1(test_dgmres_T<double>());
   CALL_SUBTEST_2(test_dgmres_T<std::complex<double> >());
+  CALL_SUBTEST_3(test_dgmres_breakdown_T<double>());
+  CALL_SUBTEST_4(test_dgmres_breakdown_T<std::complex<double> >());
 }