From f82702872dbf1d0ec301c999e30c6e8e5d2b7cee Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 27 Aug 2024 11:16:44 -0400
Subject: [PATCH] Normal residual

---
 Grid/algorithms/iterative/NormalEquations.h | 26 +++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/Grid/algorithms/iterative/NormalEquations.h b/Grid/algorithms/iterative/NormalEquations.h
index df90bee7..805937e1 100644
--- a/Grid/algorithms/iterative/NormalEquations.h
+++ b/Grid/algorithms/iterative/NormalEquations.h
@@ -60,6 +60,32 @@ public:
   }     
 };
 
+template<class Field> class NormalResidual : public LinearFunction<Field>{
+private:
+  SparseMatrixBase<Field> & _Matrix;
+  OperatorFunction<Field> & _HermitianSolver;
+  LinearFunction<Field>   & _Guess;
+public:
+
+  /////////////////////////////////////////////////////
+  // Wrap the usual normal equations trick
+  /////////////////////////////////////////////////////
+ NormalResidual(SparseMatrixBase<Field> &Matrix, OperatorFunction<Field> &HermitianSolver,
+		 LinearFunction<Field> &Guess) 
+   :  _Matrix(Matrix), _HermitianSolver(HermitianSolver), _Guess(Guess) {}; 
+
+  void operator() (const Field &in, Field &out){
+ 
+    Field res(in.Grid());
+    Field tmp(in.Grid());
+
+    MMdagLinearOperator<SparseMatrixBase<Field>,Field> MMdagOp(_Matrix);
+    _Guess(in,res);
+    _HermitianSolver(MMdagOp,in,res);  // M Mdag res = in ;
+    _Matrix.Mdag(res,out);             // out = Mdag res
+  }     
+};
+
 template<class Field> class HPDSolver : public LinearFunction<Field> {
 private:
   LinearOperatorBase<Field> & _Matrix;