From 4107b371e365a8f2d2544af9808fa90c05a81731 Mon Sep 17 00:00:00 2001
From: Desire NUENTSA <desire.nuentsa_wakam@inria.fr>
Date: Wed, 20 Mar 2013 11:22:45 +0100
Subject: [PATCH] Handle zero right hand side in CG and GMRES

---
 Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 6 ++++++
 unsupported/Eigen/src/IterativeSolvers/GMRES.h       | 3 ++-
 2 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
index f64f2534d..eadf711e5 100644
--- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
+++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
@@ -48,6 +48,12 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
   VectorType z(n), tmp(n);
   RealScalar absNew = internal::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
   RealScalar rhsNorm2 = rhs.squaredNorm();
+   // Check Zero right hand side
+  if(!rhsNorm2) 
+  {
+    x.setZero();
+    return;
+  }
   RealScalar residualNorm2 = 0;
   RealScalar threshold = tol*tol*rhsNorm2;
   int i = 0;
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 34e67db82..2efb6ff92 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -348,7 +348,8 @@ public:
   template<typename Rhs,typename Dest>
   void _solve(const Rhs& b, Dest& x) const
   {
-    x.setZero();
+    x = b;
+    if(!x.squaredNorm()) return; // Check Zero right hand side
     _solveWithGuess(b,x);
   }