From da6219b19dd92231cd0afe380ae4880b62bfe88d Mon Sep 17 00:00:00 2001
From: Desire NUENTSA <desire.nuentsa_wakam@inria.fr>
Date: Wed, 20 Mar 2013 16:15:18 +0100
Subject: [PATCH] Bug567 : Fix iterative solvers to immediately return when the
 initial guess is the true solution and for trivial solution

---
 Eigen/src/IterativeLinearSolvers/BiCGSTAB.h   |  5 ++
 .../ConjugateGradient.h                       | 28 ++++++----
 test/sparse_solver.h                          | 56 ++++++++++++-------
 .../Eigen/src/IterativeSolvers/GMRES.h        |  2 +-
 4 files changed, 59 insertions(+), 32 deletions(-)

diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
index 5a822e0ea..fbefb696f 100644
--- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
+++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
@@ -44,6 +44,11 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
   VectorType r0 = r;
   
   RealScalar r0_sqnorm = rhs.squaredNorm();
+  if(r0_sqnorm == 0)
+  {
+    x.setZero();
+    return true;
+  }
   Scalar rho    = 1;
   Scalar alpha  = 1;
   Scalar w      = 1;
diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
index eadf711e5..00b5647c6 100644
--- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
+++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
@@ -41,21 +41,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
   int n = mat.cols();
 
   VectorType residual = rhs - mat * x; //initial residual
-  VectorType p(n);
 
+  RealScalar rhsNorm2 = rhs.squaredNorm();
+  if(rhsNorm2 == 0) 
+  {
+    x.setZero();
+    iters = 0;
+    tol_error = 0;
+    return;
+  }
+  RealScalar threshold = tol*tol*rhsNorm2;
+  RealScalar residualNorm2 = residual.squaredNorm();
+  if (residualNorm2 < threshold)
+  {
+    iters = 0;
+    tol_error = sqrt(residualNorm2 / rhsNorm2);
+    return;
+  }
+  
+  VectorType p(n);
   p = precond.solve(residual);      //initial search direction
 
   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;
   while(i < maxIters)
   {
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 5a1be67e7..6b3c48274 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -176,25 +176,32 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver)
   // generate the problem
   Mat A, halfA;
   DenseMatrix dA;
-  int size = generate_sparse_spd_problem(solver, A, halfA, dA);
-
-  // generate the right hand sides
-  int rhsCols = internal::random<int>(1,16);
-  double density = (std::max)(8./(size*rhsCols), 0.1);
-  SpMat B(size,rhsCols);
-  DenseVector b = DenseVector::Random(size);
-  DenseMatrix dB(size,rhsCols);
-  initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
-  
   for (int i = 0; i < g_repeat; i++) {
+    int size = generate_sparse_spd_problem(solver, A, halfA, dA);
+
+    // generate the right hand sides
+    int rhsCols = internal::random<int>(1,16);
+    double density = (std::max)(8./(size*rhsCols), 0.1);
+    SpMat B(size,rhsCols);
+    DenseVector b = DenseVector::Random(size);
+    DenseMatrix dB(size,rhsCols);
+    initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
+  
     check_sparse_solving(solver, A,     b,  dA, b);
     check_sparse_solving(solver, halfA, b,  dA, b);
     check_sparse_solving(solver, A,     dB, dA, dB);
     check_sparse_solving(solver, halfA, dB, dA, dB);
     check_sparse_solving(solver, A,     B,  dA, dB);
     check_sparse_solving(solver, halfA, B,  dA, dB);
+    
+    // check only once
+    if(i==0)
+    {
+      b = DenseVector::Zero(size);
+      check_sparse_solving(solver, A, b, dA, b);
+    }
   }
-
+  
   // First, get the folder 
 #ifdef TEST_REAL_CASES  
   if (internal::is_same<Scalar, float>::value 
@@ -265,21 +272,28 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
 
   Mat A;
   DenseMatrix dA;
-  int size = generate_sparse_square_problem(solver, A, dA);
-
-  A.makeCompressed();
-  DenseVector b = DenseVector::Random(size);
-  DenseMatrix dB(size,rhsCols);
-  SpMat B(size,rhsCols);
-  double density = (std::max)(8./(size*rhsCols), 0.1);
-  initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
-  B.makeCompressed();
   for (int i = 0; i < g_repeat; i++) {
+    int size = generate_sparse_square_problem(solver, A, dA);
+
+    A.makeCompressed();
+    DenseVector b = DenseVector::Random(size);
+    DenseMatrix dB(size,rhsCols);
+    SpMat B(size,rhsCols);
+    double density = (std::max)(8./(size*rhsCols), 0.1);
+    initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
+    B.makeCompressed();
     check_sparse_solving(solver, A, b,  dA, b);
     check_sparse_solving(solver, A, dB, dA, dB);
     check_sparse_solving(solver, A, B,  dA, dB);
+    
+    // check only once
+    if(i==0)
+    {
+      b = DenseVector::Zero(size);
+      check_sparse_solving(solver, A, b, dA, b);
+    }
   }
-   
+  
   // First, get the folder 
 #ifdef TEST_REAL_CASES
   if (internal::is_same<Scalar, float>::value 
diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
index 2efb6ff92..f82472159 100644
--- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h
+++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h
@@ -349,7 +349,7 @@ public:
   void _solve(const Rhs& b, Dest& x) const
   {
     x = b;
-    if(!x.squaredNorm()) return; // Check Zero right hand side
+    if(x.squaredNorm() == 0) return; // Check Zero right hand side
     _solveWithGuess(b,x);
   }