diff --git a/Grid/algorithms/iterative/BlockConjugateGradient.h b/Grid/algorithms/iterative/BlockConjugateGradient.h
index 0bd7cf51..14707d7b 100644
--- a/Grid/algorithms/iterative/BlockConjugateGradient.h
+++ b/Grid/algorithms/iterative/BlockConjugateGradient.h
@@ -52,6 +52,7 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
   Integer PrintInterval; //GridLogMessages or Iterative
+  RealD TrueResidual;
   
   BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true)
     : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100)
@@ -306,7 +307,8 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
 
       Linop.HermOp(X, AD);
       AD = AD-B;
-      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(norm2(AD)/norm2(B)) <<std::endl;
+      TrueResidual = std::sqrt(norm2(AD)/norm2(B));
+      std::cout << GridLogMessage <<"\t True residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
@@ -442,7 +444,8 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
 
       Linop.HermOp(Psi, AP);
       AP = AP-Src;
-      std::cout <<GridLogMessage << "\tTrue residual is " << std::sqrt(norm2(AP)/norm2(Src)) <<std::endl;
+      TrueResidual = std::sqrt(norm2(AP)/norm2(Src));
+      std::cout <<GridLogMessage << "\tTrue residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
@@ -668,7 +671,8 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
 
       for(int b=0;b<Nblock;b++) Linop.HermOp(X[b], AD[b]);
       for(int b=0;b<Nblock;b++) AD[b] = AD[b]-B[b];
-      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(normv(AD)/normv(B)) <<std::endl;
+      TrueResidual = std::sqrt(normv(AD)/normv(B));
+      std::cout << GridLogMessage <<"\t True residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
diff --git a/Grid/algorithms/iterative/ConjugateGradient.h b/Grid/algorithms/iterative/ConjugateGradient.h
index 398f578f..160cb185 100644
--- a/Grid/algorithms/iterative/ConjugateGradient.h
+++ b/Grid/algorithms/iterative/ConjugateGradient.h
@@ -49,6 +49,7 @@ public:
   RealD Tolerance;
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  RealD TrueResidual;
   
   ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true)
     : Tolerance(tol),
@@ -82,6 +83,14 @@ public:
     cp = a;
     ssq = norm2(src);
 
+    // Handle trivial case of zero src
+    if (ssq == 0.){
+      psi = 0;
+      IterationsToComplete = 1;
+      TrueResidual = 0.;
+      return;
+    }
+
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient: guess " << guess << std::endl;
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient:   src " << ssq << std::endl;
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient:    mp " << d << std::endl;
@@ -93,6 +102,7 @@ public:
 
     // Check if guess is really REALLY good :)
     if (cp <= rsq) {
+      TrueResidual = sqrt(a/ssq);
       std::cout << GridLogMessage << "ConjugateGradient guess is converged already " << std::endl;
       IterationsToComplete = 0;	
       return;
@@ -142,7 +152,7 @@ public:
       LinalgTimer.Stop();
 
       std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k
-                << " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
+                << " residual " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
 
       // Stopping condition
       if (cp <= rsq) {
@@ -170,10 +180,17 @@ public:
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 
 	IterationsToComplete = k;	
+	TrueResidual = true_residual;
 
         return;
       }
     }
+    // Failed. Calculate true residual before giving up                                                         
+    Linop.HermOpAndNorm(psi, mmp, d, qq);
+    p = mmp - src;
+
+    TrueResidual = sqrt(norm2(p)/ssq);
+
     std::cout << GridLogMessage << "ConjugateGradient did NOT converge "<<k<<" / "<< MaxIterations<< std::endl;
 
     if (ErrorOnNoConverge) assert(0);
diff --git a/Grid/algorithms/iterative/ConjugateGradientMultiShift.h b/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
index 7e8d01e2..b8b7e056 100644
--- a/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
+++ b/Grid/algorithms/iterative/ConjugateGradientMultiShift.h
@@ -46,15 +46,19 @@ public:
 
   RealD   Tolerance;
   Integer MaxIterations;
-    Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  std::vector<int> IterationsToComplete;
+  //    Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
   int verbose;
   MultiShiftFunction shifts;
+  std::vector<RealD> TrueResiduals;
 
   ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : 
     MaxIterations(maxit),
     shifts(_shifts)
   { 
     verbose=1;
+    IterationsToComplete.resize(_shifts.order);
+    TrueResiduals.resize(_shifts.order);
   }
 
   void operator() (LinearOperatorBase<Field> &Linop, const Field &src, Field &psi)
@@ -125,6 +129,17 @@ public:
     // Residuals "r" are src
     // First search direction "p" is also src
     cp = norm2(src);
+
+    // Handle trivial case of zero src.
+    if( cp == 0. ){
+      for(int s=0;s<nshift;s++){
+	psi[s] = 0;
+	IterationsToComplete[s] = 1;
+	TrueResiduals[s] = 0.;
+      }
+      return;
+    }
+
     for(int s=0;s<nshift;s++){
       rsq[s] = cp * mresidual[s] * mresidual[s];
       std::cout<<GridLogMessage<<"ConjugateGradientMultiShift: shift "<<s
@@ -270,6 +285,7 @@ public:
       for(int s=0;s<nshift;s++){
       
 	if ( (!converged[s]) ){
+	  IterationsToComplete[s] = k;
 	
 	  RealD css  = c * z[s][iz]* z[s][iz];
 	
@@ -299,7 +315,8 @@ public:
 	  axpy(r,-alpha[s],src,tmp);
 	  RealD rn = norm2(r);
 	  RealD cn = norm2(src);
-	  std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<<std::sqrt(rn/cn)<<std::endl;
+	  TrueResiduals[s] = std::sqrt(rn/cn);
+	  std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<< TrueResiduals[s] <<std::endl;
 	}
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
@@ -308,7 +325,7 @@ public:
       std::cout << GridLogMessage << "\tMarix    " << MatrixTimer.Elapsed()     <<std::endl;
       std::cout << GridLogMessage << "\tShift    " << ShiftTimer.Elapsed()     <<std::endl;
 
-      IterationsToComplete = k;	
+      //      IterationsToComplete = k;	
 
 	return;
       }