Working towards solvers

lib/algorithms/LinearOperator.h

@@ -5,6 +5,22 @@
 
 namespace Grid {
 
+  /////////////////////////////////////////////////////////////////////////////////////////////
+  // Interface defining what I expect of a general sparse matrix, such as a Fermion action
+  /////////////////////////////////////////////////////////////////////////////////////////////
+    template<class Field> class SparseMatrixBase {
+    public:
+      // Full checkerboar operations
+      virtual void M    (const Field &in, Field &out);
+      virtual void Mdag (const Field &in, Field &out);
+      virtual RealD MdagM(const Field &in, Field &out);
+      
+      // half checkerboard operaions
+      virtual void Mpc      (const Field &in, Field &out);
+      virtual void MpcDag   (const Field &in, Field &out);
+      virtual RealD MpcDagMpc(const Field &in, Field &out);
+    };
+
   /////////////////////////////////////////////////////////////////////////////////////////////
   // LinearOperators Take a something and return a something.
   /////////////////////////////////////////////////////////////////////////////////////////////
@@ -27,25 +43,13 @@ namespace Grid {
   /////////////////////////////////////////////////////////////////////////////////////////////
     template<class Field> class HermitianOperatorBase : public LinearOperatorBase<Field> {
     public:
+      virtual RealD OpAndNorm(const Field &in, Field &out);
       void AdjOp(const Field &in, Field &out) {
-	return Op(in,out);
+	Op(in,out);
+      };
+      void Op(const Field &in, Field &out) {
+	OpAndNorm(in,out);
       };
-    };
-
-  /////////////////////////////////////////////////////////////////////////////////////////////
-  // Interface defining what I expect of a general sparse matrix
-  /////////////////////////////////////////////////////////////////////////////////////////////
-    template<class Field> class SparseMatrixBase {
-    public:
-      // Full checkerboar operations
-      virtual void M    (const Field &in, Field &out);
-      virtual void Mdag (const Field &in, Field &out);
-      virtual void MdagM(const Field &in, Field &out);
-      
-      // half checkerboard operaions
-      virtual void Mpc      (const Field &in, Field &out);
-      virtual void MpcDag   (const Field &in, Field &out);
-      virtual void MpcDagMpc(const Field &in, Field &out);
     };
 
   /////////////////////////////////////////////////////////////////////////////////////////////
@@ -79,10 +83,10 @@ namespace Grid {
     public:
       NonHermitianRedBlackOperator(SparseMatrix &Mat): _Mat(Mat){};
       void Op     (const Field &in, Field &out){
-	this->Mpc(in,out);
+	_Mat.Mpc(in,out);
       }
       void AdjOp     (const Field &in, Field &out){ //
-	this->MpcDag(in,out);
+	_Mat.MpcDag(in,out);
       }
     };
 
@@ -94,12 +98,8 @@ namespace Grid {
       SparseMatrix &_Mat;
     public:
       HermitianOperator(SparseMatrix &Mat): _Mat(Mat) {};
-
+      RealD OpAndNorm(const Field &in, Field &out){
-      void Op     (const Field &in, Field &out){
+	return _Mat.MdagM(in,out);
-	this->Mpc(in,out);
-      }
-      void AdjOp     (const Field &in, Field &out){ //
-	this->MpcDag(in,out);
       }
     };
 
@@ -111,8 +111,8 @@ namespace Grid {
       SparseMatrix &_Mat;
     public:
       HermitianRedBlackOperator(SparseMatrix &Mat): _Mat(Mat) {};
-      void Op     (const Field &in, Field &out){
+      RealD OpAndNorm(const Field &in, Field &out){
-	this->MpcDagMpc(in,out);
+	return _Mat.MpcDagMpc(in,out);
       }
     };
 
@@ -141,6 +141,7 @@ namespace Grid {
     public:
       RealD   Tolerance;
       Integer MaxIterations;
+      IterativeProcess(RealD tol,Integer maxit) : Tolerance(tol),MaxIterations(maxit) {};
       virtual void operator() (LinearOperatorBase<Field> *Linop,const Field &in, Field &out) = 0;
     };
 

lib/algorithms/approx/PolynomialApprox.h

@@ -1,6 +1,9 @@
 #ifndef GRID_POLYNOMIAL_APPROX_H
 #define GRID_POLYNOMIAL_APPROX_H
 
+#include<Grid.h>
+#include<algorithms/LinearOperator.h>
+
 namespace Grid {
 
   ////////////////////////////////////////////////////////////////////////////////////////////

lib/algorithms/iterative/ConjugateGradient.h

@@ -0,0 +1,88 @@
+#ifndef GRID_CONJUGATE_GRADIENT_H
+#define GRID_CONJUGATE_GRADIENT_H
+
+#include<Grid.h>
+#include<algorithms/LinearOperator.h>
+
+namespace Grid {
+
+  template<class Field> class ConjugateGradient : public IterativeProcess<Field> {
+  public:
+
+    ConjugateGradient(RealD tol,Integer maxit): IterativeProces(tol,maxit) {};
+
+    void operator() (HermitianOperatorBase<Field> *Linop,const Field &src, Field &psi){
+
+      RealD residual = Tolerance;
+      RealD cp,c,a,d,b,ssq,qq,b_pred;
+      
+      Field   p(src);
+      Field mmp(src);
+      Field   r(src);
+      
+      //Initial residual computation & set up
+      RealD guess = norm2(psi);
+      
+      Linop.OpAndNorm(psi,mmp,d,b);
+      
+      r= src-mmp;
+      p= r;
+      
+      a  =norm2(p);
+      cp =a;
+      ssq=norm2(src);
+      
+      std::cout <<setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
+      std::cout <<setprecision(4)<< "ConjugateGradient:   src "<<ssq  <<std::endl;
+      std::cout <<setprecision(4)<< "ConjugateGradient:    mp "<<d    <<std::endl;
+      std::cout <<setprecision(4)<< "ConjugateGradient:   mmp "<<b    <<std::endl;
+      std::cout <<setprecision(4)<< "ConjugateGradient:     r "<<cp   <<std::endl;
+      std::cout <<setprecision(4)<< "ConjugateGradient:     p "<<a    <<std::endl;
+      
+      RealD rsq =  residual* residual*ssq;
+      
+      //Check if guess is really REALLY good :)
+      if ( cp <= rsq ) {
+	return 0;
+      }
+      
+      std::cout << setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
+      
+      int k;
+      for (k=1;k<=max_iter;k++){
+	
+	c=cp;
+	
+	Linop.OpAndNorm(p,mmp,d,qq);
+	  
+	a      = c/d;
+	b_pred = a*(a*qq-d)/c;
+	
+	cp = axpy_norm(r,mmp,r,-a);
+	b = cp/c;
+	
+	// Fuse these loops ; should be really easy
+	// Update psi; New (conjugate/M-orthogonal) search direction
+	psi= a*p+psi;
+	p  = p*b+r;
+	  
+	// Stopping condition
+	if ( cp <= rsq ) { 
+	  
+	  Linop.Op(p,mmp);
+	  p=mmp-src;
+	  
+	  RealD mmpnorm = sqrt(norm2(mmp));
+	  RealD psinorm = sqrt(norm2(psi));
+	  RealD srcnorm = sqrt(norm2(src));
+	  RealD resnorm = sqrt(norm2(p));
+	  RealD true_residual = resnorm/srcnorm;
+	  std::cout<<"ConjugateGradient: true   residual  is "<<true_residual<<" sol "<<psinorm<<" src "<<srcnorm<<std::endl;
+	  std::cout<<"ConjugateGradient: target residual was "<<residual<<std::endl;
+	  return k;
+	}
+      }
+    }
+  };
+}
+#endif