DGMRES: fix null rhs, fix restart, fix m_isDeflInitialized for multiple solve

1 files changed, 18 insertions, 4 deletions

diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
index 2aea53dcb..9ae766995 100644
--- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h
@@ -241,7 +241,18 @@ template<typename Rhs, typename Dest>
 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
               const Preconditioner& precond) const
 {
+  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
+
+  RealScalar normRhs = rhs.norm();
+  if(normRhs <= considerAsZero) 
+  {
+    x.setZero();
+    m_error = 0;
+    return;
+  }
+
   //Initialization
+  m_isDeflInitialized = false;
   Index n = mat.rows(); 
   DenseVector r0(n); 
   Index nbIts = 0; 
@@ -249,10 +260,11 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh
   m_Hes.resize(m_restart, m_restart);
   m_V.resize(n,m_restart+1);
   //Initial residual vector and initial norm
-  x = precond.solve(x);
+  if(x.squaredNorm()==0) 
+    x = precond.solve(rhs);
   r0 = rhs - mat * x; 
   RealScalar beta = r0.norm(); 
-  RealScalar normRhs = rhs.norm();
+  
   m_error = beta/normRhs; 
   if(m_error < m_tolerance)
     m_info = Success; 
@@ -265,8 +277,10 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh
     dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 
     
     // Compute the new residual vector for the restart 
-    if (nbIts < m_iterations && m_info == NoConvergence)
-      r0 = rhs - mat * x; 
+    if (nbIts < m_iterations && m_info == NoConvergence) {
+      r0 = rhs - mat * x;
+      beta = r0.norm();
+    }
   }
 }