Conjugate residual added

lib/Algorithms.h

@@ -5,6 +5,7 @@
 #include <algorithms/LinearOperator.h>
 
 #include <algorithms/iterative/ConjugateGradient.h>
+#include <algorithms/iterative/ConjugateResidual.h>
 #include <algorithms/iterative/NormalEquations.h>
 #include <algorithms/iterative/SchurRedBlack.h>
 

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/Dirac.h ./qcd/LinalgUtils.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
 
-CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/Dirac.cc ./qcd/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
+CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/algorithms/LinearOperator.h

@@ -18,7 +18,8 @@ namespace Grid {
     public:
       virtual void Op     (const Field &in, Field &out) = 0; // Abstract base
       virtual void AdjOp  (const Field &in, Field &out) = 0; // Abstract base
-      virtual void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2)=0;
+      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0;
+      virtual void HermOp(const Field &in, Field &out)=0;
     };
 
 
@@ -48,9 +49,13 @@ namespace Grid {
       void AdjOp     (const Field &in, Field &out){
 	_Mat.Mdag(in,out);
       }
-      void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
+      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
 	_Mat.MdagM(in,out,n1,n2);
       }
+      void HermOp(const Field &in, Field &out){
+	RealD n1,n2;
+	HermOpAndNorm(in,out,n1,n2);
+      }
     };
 
     ////////////////////////////////////////////////////////////////////
@@ -67,7 +72,7 @@ namespace Grid {
       void AdjOp     (const Field &in, Field &out){
 	_Mat.M(in,out);
       }
-      void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
+      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
 	ComplexD dot;
 
 	_Mat.M(in,out);
@@ -78,6 +83,9 @@ namespace Grid {
 	dot = innerProduct(out,out);
 	n2=real(dot);
       }
+      void HermOp(const Field &in, Field &out){
+	_Mat.M(in,out);
+      }
     };
 
     //////////////////////////////////////////////////////////
@@ -98,6 +106,10 @@ namespace Grid {
       void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
 	MpcDagMpc(in,out,n1,n2);
       }
+      void HermOp(const Field &in, Field &out){
+	RealD n1,n2;
+	HermOpAndNorm(in,out,n1,n2);
+      }
       void Op     (const Field &in, Field &out){
 	Mpc(in,out);
       }

lib/algorithms/SparseMatrix.h

@@ -18,6 +18,7 @@ namespace Grid {
 	Field tmp (in._grid);
 	ni=M(in,tmp);
 	no=Mdag(tmp,out);
+	std::cout << "MdagM "<< ni<<" "<<no<<std::endl;
       }
     };
 

lib/algorithms/iterative/ConjugateResidual.h

@@ -0,0 +1,85 @@
+#ifndef GRID_CONJUGATE_RESIDUAL_H
+#define GRID_CONJUGATE_RESIDUAL_H
+
+namespace Grid {
+
+    /////////////////////////////////////////////////////////////
+    // Base classes for iterative processes based on operators
+    // single input vec, single output vec.
+    /////////////////////////////////////////////////////////////
+
+  template<class Field> 
+    class ConjugateResidual : public OperatorFunction<Field> {
+  public:                                                
+    RealD   Tolerance;
+    Integer MaxIterations;
+    int verbose;
+
+    ConjugateResidual(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { 
+      verbose=1;
+    };
+
+    void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi){
+
+      RealD a, b, c, d;
+      RealD cp, ssq,rsq;
+      
+      RealD rAr, rAAr, rArp;
+      RealD pAp, pAAp;
+
+      GridBase *grid = src._grid;
+      psi=zero;
+      Field r(grid),  p(grid), Ap(grid), Ar(grid);
+      
+      r=src;
+      p=src;
+
+      Linop.HermOpAndNorm(p,Ap,pAp,pAAp);
+      Linop.HermOpAndNorm(r,Ar,rAr,rAAr);
+
+      std::cout << "pAp, pAAp"<< pAp<<" "<<pAAp<<std::endl;
+      std::cout << "rAr, rAAr"<< rAr<<" "<<rAAr<<std::endl;
+
+      cp =norm2(r);
+      ssq=norm2(src);
+      rsq=Tolerance*Tolerance*ssq;
+
+      std::cout<<"ConjugateResidual: iteration " <<0<<" residual "<<cp<< " target"<< rsq<<std::endl;
+
+      for(int k=1;k<MaxIterations;k++){
+
+	a = rAr/pAAp;
+
+	axpy(psi,a,p,psi);
+
+	cp = axpy_norm(r,-a,Ap,r);
+
+	rArp=rAr;
+
+	Linop.HermOpAndNorm(r,Ar,rAr,rAAr);
+
+	b   =rAr/rArp;
+ 
+	axpy(p,b,p,r);
+	pAAp=axpy_norm(Ap,b,Ap,Ar);
+
+	std::cout<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
+
+	if(cp<rsq) {
+	  Linop.HermOp(psi,Ap);
+	  axpy(r,-1.0,src,Ap);
+	  RealD true_resid = norm2(r);
+	  std::cout<<"ConjugateResidual: Converged on iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
+	  std::cout<<"ConjugateResidual: true   residual  is "<<true_resid<<std::endl;
+	  std::cout<<"ConjugateResidual: target residual was "<<Tolerance <<std::endl;
+	  return;
+	}
+
+      }
+
+      std::cout<<"ConjugateResidual did NOT converge"<<std::endl;
+      assert(0);
+    }
+  };
+}
+#endif

lib/qcd/QCD.h

@@ -307,10 +307,10 @@ namespace QCD {
 }   //namespace QCD
 } // Grid
 
-#include <qcd/SpaceTimeGrid.h>
+#include <qcd/utils/SpaceTimeGrid.h>
-#include <qcd/Dirac.h>
+#include <qcd/spin/Dirac.h>
-#include <qcd/TwoSpinor.h>
+#include <qcd/spin/TwoSpinor.h>
-#include <qcd/LinalgUtils.h>
+#include <qcd/utils/LinalgUtils.h>
 #include <qcd/utils/CovariantCshift.h>
 #include <qcd/utils/WilsonLoops.h>
 #include <qcd/action/Actions.h>

lib/qcd/action/Actions.h

@@ -65,36 +65,9 @@
 #include <qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h>
 #include <qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h>
 
-
+///////////////////////////////////////////////////////////////////////////////
-    // Chroma interface defining FermionAction
+// G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code
-    /*
+///////////////////////////////////////////////////////////////////////////////
-     template<typename T, typename P, typename Q>  class FermAct4D : public FermionAction<T,P,Q>
+#include <qcd/action/fermion/g5HermitianLinop.h>
-     virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
-     virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
-     virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
-     class DiffFermAct4D : public FermAct4D<T,P,Q>
-     virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-    */
-
-
-    // Chroma interface defining GaugeAction
-    /*
-      template<typename P, typename Q>   class GaugeAction
-  virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
-  virtual GaugeState<P,Q>* createState(const Q& q) const
-  virtual const GaugeBC<P,Q>& getGaugeBC() const
-  virtual const Set& getSet(void) const = 0;
-  virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const 
-  virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
-
-  class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
-  typedef multi1d<LatticeColorMatrix>  P;
-    */
 
 #endif

lib/qcd/action/fermion/g5HermitianLinop.h

@@ -0,0 +1,38 @@
+#ifndef G5_HERMITIAN_LINOP
+#define G5_HERMITIAN_LINOP
+namespace Grid {
+  namespace QCD {
+////////////////////////////////////////////////////////////////////
+// Wrap an already herm matrix
+////////////////////////////////////////////////////////////////////
+template<class Matrix,class Field>
+class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+public:
+  Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
+  void Op     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+
+    HermOp(in,out);
+    
+    ComplexD dot;
+    dot= innerProduct(in,out);
+    n1=real(dot);
+    
+    dot = innerProduct(out,out);
+    n2=real(dot);
+  }
+  void HermOp(const Field &in, Field &out){
+    Field tmp(in._grid);
+    _Mat.M(in,tmp);
+    G5R5(out,tmp);
+  }
+};
+
+}}
+#endif

lib/qcd/utils/LinalgUtils.h

@@ -109,5 +109,24 @@ PARALLEL_FOR_LOOP
     vstream(z._odata[ss+s],tmp);
   }
 }
+
+template<class vobj> 
+void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
+{
+  GridBase *grid=x._grid;
+  z.checkerboard = x.checkerboard;
+  conformable(x,z);
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    for(int s=0;s<Ls;s++){
+      int sp = Ls-1-s;
+      multGamma5(tmp(),x._odata[ss+s]());
+      vstream(z._odata[ss+sp],tmp);
+    }
+  }
+}
+
 }}
 #endif 

tests/Make.inc

@@ -14,6 +14,10 @@ Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
 Test_cayley_even_odd_LDADD=-lGrid
 
 
+Test_cf_cr_unprec_SOURCES=Test_cf_cr_unprec.cc
+Test_cf_cr_unprec_LDADD=-lGrid
+
+
 Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc
 Test_contfrac_cg_LDADD=-lGrid
 
@@ -42,6 +46,10 @@ Test_dwf_cg_unprec_SOURCES=Test_dwf_cg_unprec.cc
 Test_dwf_cg_unprec_LDADD=-lGrid
 
 
+Test_dwf_cr_unprec_SOURCES=Test_dwf_cr_unprec.cc
+Test_dwf_cr_unprec_LDADD=-lGrid
+
+
 Test_dwf_even_odd_SOURCES=Test_dwf_even_odd.cc
 Test_dwf_even_odd_LDADD=-lGrid
 
@@ -90,6 +98,10 @@ Test_wilson_cg_unprec_SOURCES=Test_wilson_cg_unprec.cc
 Test_wilson_cg_unprec_LDADD=-lGrid
 
 
+Test_wilson_cr_unprec_SOURCES=Test_wilson_cr_unprec.cc
+Test_wilson_cr_unprec_LDADD=-lGrid
+
+
 Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc
 Test_wilson_even_odd_LDADD=-lGrid