First pass at continued fraction; solver and even odd decomposition tests pass.

Have to make ContFrac class virtual and derive end non-abstract actions for the particular
cases.

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -229,7 +229,14 @@ namespace QCD {
     }
   }
   
-  void CayleyFermion5D::SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c)
+  // Tanh
+  void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
+  {
+    SetCoefficientsZolotarev(1.0,zdata,b,c);
+
+  }
+  //Zolo
+  void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
   {
 
     ///////////////////////////////////////////////////////////
@@ -266,7 +273,7 @@ namespace QCD {
     double bmc = b-c;
     for(int i=0; i < Ls; i++){
       as[i] = 1.0;
-      omega[i] = ((double)zdata->gamma[i]); //NB reciprocal relative to Chroma NEF code
+      omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code
       bs[i] = 0.5*(bpc/omega[i] + bmc);
       cs[i] = 0.5*(bpc/omega[i] - bmc);
     }

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -20,7 +20,7 @@ namespace Grid {
       virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
-
+      virtual void   Instantiatable(void)=0;
       //    protected:
       RealD mass;
 
@@ -52,7 +52,8 @@ namespace Grid {
 		      RealD _mass,RealD _M5);
 
     protected:
-      void SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c);
+      void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
     };
 
   }

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -1,9 +1,56 @@
 #include <Grid.h>
 
 namespace Grid {
-
   namespace QCD {
 
+    void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
+    {
+      SetCoefficientsZolotarev(1.0,zdata,b,c);
+    }
+    void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
+    {
+      R=(1+this->mass)/(1-this->mass);
+
+      Beta.resize(Ls);
+      cc.resize(Ls);
+      cc_d.resize(Ls);
+      sqrt_cc.resize(Ls);
+      for(int i=0; i < Ls ; i++){
+	Beta[i] = zdata -> beta[i];
+	cc[i] = 1.0/Beta[i];
+	cc_d[i]=sqrt(cc[i]);
+      }
+    
+      cc_d[Ls-1]=1.0;
+      for(int i=0; i < Ls-1 ; i++){
+	sqrt_cc[i]= sqrt(cc[i]*cc[i+1]);
+      }    
+      sqrt_cc[Ls-2]=sqrt(cc[Ls-2]);
+
+
+      ZoloHiInv =1.0/zolo_hi;
+      double dw_diag = (4.0-M5)*ZoloHiInv;
+    
+      See.resize(Ls);
+      Aee.resize(Ls);
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	Aee[s] = sign * Beta[s] * dw_diag;
+	sign   = - sign;
+      }
+      Aee[Ls-1] += R;
+    
+      See[0] = Aee[0];
+      for(int s=1;s<Ls;s++){
+	See[s] = Aee[s] - 1.0/See[s-1];
+      }
+      for(int s=0;s<Ls;s++){
+	std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
+      }
+    }
+
+
+
     RealD  ContinuedFractionFermion5D::M           (const LatticeFermion &psi, LatticeFermion &chi)
     {
       LatticeFermion D(psi._grid);
@@ -13,13 +60,13 @@ namespace Grid {
       int sign=1;
       for(int s=0;s<Ls;s++){
 	if ( s==0 ) {
-	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*scale,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
+	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*ZoloHiInv,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
 	} else if ( s==(Ls-1) ){
 	  RealD R=(1.0+mass)/(1.0-mass);
-	  ag5xpby_ssp(chi,Beta[s]*scale,D,sqrt_cc[s-1],psi,s,s-1);
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,D,sqrt_cc[s-1],psi,s,s-1);
 	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
 	} else {
-	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*scale,D,sqrt_cc[s],psi,s,s+1);
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,D,sqrt_cc[s],psi,s,s+1);
   	  axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
 	}
 	sign=-sign; 
@@ -35,18 +82,22 @@ namespace Grid {
     }
     void   ContinuedFractionFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
     {
-      Dhop(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      // Apply 4d dslash
+      if ( psi.checkerboard == Odd ) {
+	DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      } else {
+	DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      }
       
       int sign=1;
       for(int s=0;s<Ls;s++){
 	if ( s==(Ls-1) ){
-	  ag5xpby_ssp(chi,Beta[s]*scale,chi,0.0,chi,s,s);
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,chi,0.0,chi,s,s);
 	} else {
-	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*scale,chi,0.0,chi,s,s);
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,chi,0.0,chi,s,s);
 	}
 	sign=-sign; 
-    }
+      }
-
     }
     void   ContinuedFractionFermion5D::MeooeDag    (const LatticeFermion &psi, LatticeFermion &chi)
     {
@@ -54,7 +105,7 @@ namespace Grid {
     }
     void   ContinuedFractionFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
     {
-      double dw_diag = (4.0-this->M5)*scale;
+      double dw_diag = (4.0-M5)*ZoloHiInv;
     
       int sign=1;
       for(int s=0;s<Ls;s++){
@@ -62,7 +113,7 @@ namespace Grid {
 	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*dw_diag,psi,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
 	} else if ( s==(Ls-1) ){
 	  // Drop the CC here.
-	  double R=(1+this->mass)/(1-this->mass);
+	  double R=(1+mass)/(1-mass);
 	  ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1);
 	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
 	} else {
@@ -80,7 +131,7 @@ namespace Grid {
     void   ContinuedFractionFermion5D::MooeeInv    (const LatticeFermion &psi, LatticeFermion &chi)
     {
       // Apply Linv
-      axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0);
+      axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); 
       for(int s=1;s<Ls;s++){
 	axpbg5y_ssp(chi,1.0/cc_d[s],psi,-1.0/See[s-1],chi,s,s-1);
       }
@@ -89,7 +140,7 @@ namespace Grid {
 	ag5xpby_ssp(chi,1.0/See[s],chi,0.0,chi,s,s); //only appearance of See[0]
       }
       // Apply Uinv = (Linv)^T
-      axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,this->Ls-1,this->Ls-1);
+      axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,Ls-1,Ls-1);
       for(int s=Ls-2;s>=0;s--){
 	axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1);
       }
@@ -112,6 +163,10 @@ namespace Grid {
 		      FourDimGrid, FourDimRedBlackGrid,M5),
       mass(_mass)
     {
+      assert((Ls&0x1)==1); // Odd Ls required
+      int nrational=Ls-1;// Even rational order
+      zdata = Approx::grid_higham(1.0,nrational);// eps is ignored for higham
+      SetCoefficientsTanh(zdata,1.0,0.0);
     }
 
   }

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -21,20 +21,8 @@ namespace Grid {
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
-    private:
+      //      virtual void   Instantiatable(void)=0;
-
+      virtual void   Instantiatable(void) {};
-      Approx::zolotarev_data *zdata;
-
-      // Cont frac
-      RealD mass;
-      RealD R;
-      RealD scale;
-      std::vector<double> Beta;
-      std::vector<double> cc;;
-      std::vector<double> cc_d;;
-      std::vector<double> sqrt_cc;
-      std::vector<double> See;
-      std::vector<double> Aee;
 
       // Constructors
       ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
@@ -44,6 +32,24 @@ namespace Grid {
 				 GridRedBlackCartesian &FourDimRedBlackGrid,
 				 RealD _mass,RealD M5);
 
+    protected:
+
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
+      void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c);
+
+      Approx::zolotarev_data *zdata;
+
+      // Cont frac
+      RealD mass;
+      RealD R;
+      RealD ZoloHiInv;
+      std::vector<double> Beta;
+      std::vector<double> cc;;
+      std::vector<double> cc_d;;
+      std::vector<double> sqrt_cc;
+      std::vector<double> See;
+      std::vector<double> Aee;
+
     };
 
 

lib/qcd/action/fermion/DomainWallFermion.h

@@ -11,6 +11,7 @@ namespace Grid {
     {
     public:
 
+      virtual void   Instantiatable(void) {};
       // Constructors
       DomainWallFermion(LatticeGaugeField &_Umu,
 			GridCartesian         &FiveDimGrid,
@@ -33,7 +34,7 @@ namespace Grid {
 	
 	std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
 	// Call base setter
-	this->CayleyFermion5D::SetCoefficients(1.0,zdata,1.0,0.0);
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);
  
       }
 

lib/qcd/action/fermion/MobiusFermion.h

@@ -11,6 +11,7 @@ namespace Grid {
     {
     public:
 
+      virtual void   Instantiatable(void) {};
       // Constructors
       MobiusFermion(LatticeGaugeField &_Umu,
 		    GridCartesian         &FiveDimGrid,
@@ -34,7 +35,7 @@ namespace Grid {
 	assert(zdata->n==this->Ls);
 	
 	// Call base setter
-	this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c);
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c);
  
       }
 

lib/qcd/action/fermion/MobiusZolotarevFermion.h

@@ -11,6 +11,7 @@ namespace Grid {
     {
     public:
 
+      virtual void   Instantiatable(void) {};
       // Constructors
        MobiusZolotarevFermion(LatticeGaugeField &_Umu,
 			      GridCartesian         &FiveDimGrid,
@@ -34,10 +35,9 @@ namespace Grid {
 	assert(zdata->n==this->Ls);
 
 	std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
-	std::cout << "MobiusZolotarevFermion : note there is a degeneracy between (b+c) and Zolo param hi"<<std::endl;
 	
 	// Call base setter
-	this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c);
+	this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);
  
       }
 

tests/Make.inc

@@ -1,5 +1,21 @@
 
-bin_PROGRAMS = Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_many_cg Test_many_evenodd Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_evenodd
+bin_PROGRAMS = Test_cayley_cg Test_cayley_even_odd Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_even_odd
+
+
+Test_cayley_cg_SOURCES=Test_cayley_cg.cc
+Test_cayley_cg_LDADD=-lGrid
+
+
+Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
+Test_cayley_even_odd_LDADD=-lGrid
+
+
+Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc
+Test_contfrac_cg_LDADD=-lGrid
+
+
+Test_contfrac_even_odd_SOURCES=Test_contfrac_even_odd.cc
+Test_contfrac_even_odd_LDADD=-lGrid
 
 
 Test_cshift_SOURCES=Test_cshift.cc
@@ -34,14 +50,6 @@ Test_main_SOURCES=Test_main.cc
 Test_main_LDADD=-lGrid
 
 
-Test_many_cg_SOURCES=Test_many_cg.cc
-Test_many_cg_LDADD=-lGrid
-
-
-Test_many_evenodd_SOURCES=Test_many_evenodd.cc
-Test_many_evenodd_LDADD=-lGrid
-
-
 Test_nersc_io_SOURCES=Test_nersc_io.cc
 Test_nersc_io_LDADD=-lGrid
 
@@ -78,6 +86,6 @@ Test_wilson_cg_unprec_SOURCES=Test_wilson_cg_unprec.cc
 Test_wilson_cg_unprec_LDADD=-lGrid
 
 
-Test_wilson_evenodd_SOURCES=Test_wilson_evenodd.cc
+Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc
-Test_wilson_evenodd_LDADD=-lGrid
+Test_wilson_even_odd_LDADD=-lGrid
 

tests/Sqrt.gnu

@@ -1,2 +0,0 @@
-f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866));
-f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422));

tests/Test_contfrac_cg.cc

@@ -0,0 +1,147 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+
+template<class What> 
+void  TestCGinversions(What & Ddwf, 
+		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		       RealD mass, RealD M5,
+		       GridParallelRNG *RNG4,
+		       GridParallelRNG *RNG5);
+template<class What> 
+void  TestCGschur(What & Ddwf, 
+		  GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		  GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		  RealD mass, RealD M5,
+		  GridParallelRNG *RNG4,
+		  GridParallelRNG *RNG5);
+
+template<class What> 
+void  TestCGunprec(What & Ddwf, 
+		   GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		   GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		   RealD mass, RealD M5,
+		   GridParallelRNG *RNG4,
+		   GridParallelRNG *RNG5);
+
+template<class What> 
+void  TestCGprec(What & Ddwf, 
+		 GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		 GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		 RealD mass, RealD M5,
+		 GridParallelRNG *RNG4,
+		 GridParallelRNG *RNG5);
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  const int Ls=9;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  std::cout <<"ContinuedFractionFermion test"<<std::endl;
+  ContinuedFractionFermion5D Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  TestCGinversions<ContinuedFractionFermion5D>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
+  Grid_finalize();
+}
+template<class What> 
+void  TestCGinversions(What & Ddwf, 
+		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		       RealD mass, RealD M5,
+		       GridParallelRNG *RNG4,
+		       GridParallelRNG *RNG5)
+{
+  std::cout << "Testing unpreconditioned inverter"<<std::endl;
+  TestCGunprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
+  std::cout << "Testing red black preconditioned inverter"<<std::endl;
+  TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
+  std::cout << "Testing red black Schur inverter"<<std::endl;
+  TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
+}
+
+template<class What> 
+void  TestCGunprec(What & Ddwf, 
+		   GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		   GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		   RealD mass, RealD M5,
+		   GridParallelRNG *RNG4,
+		   GridParallelRNG *RNG5)
+{
+  LatticeFermion src   (FGrid); random(*RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+
+  HermitianOperator<What,LatticeFermion> HermOp(Ddwf);
+  ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+  CG(HermOp,src,result);
+
+}
+template<class What> 
+void  TestCGprec(What & Ddwf, 
+		 GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		 GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		 RealD mass, RealD M5,
+		 GridParallelRNG *RNG4,
+		 GridParallelRNG *RNG5)
+{
+  LatticeFermion src   (FGrid); random(*RNG5,src);
+  LatticeFermion    src_o(FrbGrid);
+  LatticeFermion result_o(FrbGrid);
+  pickCheckerboard(Odd,src_o,src);
+  result_o=zero;
+
+  HermitianCheckerBoardedOperator<What,LatticeFermion> HermOpEO(Ddwf);
+  ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+  CG(HermOpEO,src_o,result_o);
+}
+
+
+template<class What> 
+void  TestCGschur(What & Ddwf, 
+		   GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+		   GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		   RealD mass, RealD M5,
+		   GridParallelRNG *RNG4,
+		   GridParallelRNG *RNG5)
+{
+  LatticeFermion src   (FGrid); random(*RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+
+  ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+  SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
+  SchurSolver(Ddwf,src,result);
+}

tests/Test_contfrac_even_odd.cc

@@ -0,0 +1,218 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+
+template<class What> 
+void  TestWhat(What & Ddwf,
+	       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+	       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+	       RealD mass, RealD M5,
+	       GridParallelRNG *RNG4,   GridParallelRNG *RNG5);
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  const int Ls=9;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  std::cout <<"ContinuedFractionFermion test"<<std::endl;
+  ContinuedFractionFermion5D Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  TestWhat<ContinuedFractionFermion5D>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
+  Grid_finalize();
+}
+
+template<class What> 
+void  TestWhat(What & Ddwf, 
+	       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+	       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+	       RealD mass, RealD M5,
+	       GridParallelRNG *RNG4,
+	       GridParallelRNG *RNG5)
+{
+
+  LatticeFermion src   (FGrid); random(*RNG5,src);
+  LatticeFermion phi   (FGrid); random(*RNG5,phi);
+  LatticeFermion chi   (FGrid); random(*RNG5,chi);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid);    ref=zero;
+  LatticeFermion    tmp(FGrid);    tmp=zero;
+  LatticeFermion    err(FGrid);    tmp=zero;
+
+  LatticeFermion src_e (FrbGrid);
+  LatticeFermion src_o (FrbGrid);
+  LatticeFermion r_e   (FrbGrid);
+  LatticeFermion r_o   (FrbGrid);
+  LatticeFermion r_eo  (FGrid);
+  LatticeFermion r_eeoo(FGrid);
+
+  std::cout<<"=========================================================="<<std::endl;
+  std::cout<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
+  std::cout<<"=========================================================="<<std::endl;
+
+  pickCheckerboard(Even,src_e,src);
+  pickCheckerboard( Odd,src_o,src);
+
+  Ddwf.Meooe(src_e,r_o);  std::cout<<"Applied Meo "<<norm2(r_o)<<std::endl;
+  Ddwf.Meooe(src_o,r_e);  std::cout<<"Applied Moe "<<norm2(r_e)<<std::endl;
+  setCheckerboard(r_eo,r_o);
+  setCheckerboard(r_eo,r_e);
+
+  Ddwf.Mooee(src_e,r_e);  std::cout<<"Applied Mee"<<norm2(r_e)<<std::endl;
+  Ddwf.Mooee(src_o,r_o);  std::cout<<"Applied Moo"<<norm2(r_o)<<std::endl;
+  setCheckerboard(r_eeoo,r_e);
+  setCheckerboard(r_eeoo,r_o);
+
+  r_eo=r_eo+r_eeoo;
+  Ddwf.M(src,ref);  
+
+  //  std::cout << r_eo<<std::endl;
+  //  std::cout << ref <<std::endl;
+
+  err= ref - r_eo;
+  std::cout << "EO norm diff   "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
+    
+  LatticeComplex cerr(FGrid);
+  cerr = localInnerProduct(err,err);
+  //  std::cout << cerr<<std::endl;
+
+  std::cout<<"=============================================================="<<std::endl;
+  std::cout<<"= Test Ddagger is the dagger of D by requiring                "<<std::endl;
+  std::cout<<"=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  "<<std::endl;
+  std::cout<<"=============================================================="<<std::endl;
+  
+  LatticeFermion chi_e   (FrbGrid);
+  LatticeFermion chi_o   (FrbGrid);
+
+  LatticeFermion dchi_e  (FrbGrid);
+  LatticeFermion dchi_o  (FrbGrid);
+
+  LatticeFermion phi_e   (FrbGrid);
+  LatticeFermion phi_o   (FrbGrid);
+
+  LatticeFermion dphi_e  (FrbGrid);
+  LatticeFermion dphi_o  (FrbGrid);
+
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+
+  Ddwf.Meooe(chi_e,dchi_o);
+  Ddwf.Meooe(chi_o,dchi_e);
+  Ddwf.MeooeDag(phi_e,dphi_o);
+  Ddwf.MeooeDag(phi_o,dphi_e);
+
+  ComplexD pDce = innerProduct(phi_e,dchi_e);
+  ComplexD pDco = innerProduct(phi_o,dchi_o);
+  ComplexD cDpe = innerProduct(chi_e,dphi_e);
+  ComplexD cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
+  std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
+
+  std::cout<<"=============================================================="<<std::endl;
+  std::cout<<"= Test MeeInv Mee = 1                                         "<<std::endl;
+  std::cout<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.Mooee(chi_e,src_e);
+  Ddwf.MooeeInv(src_e,phi_e);
+
+  Ddwf.Mooee(chi_o,src_o);
+  Ddwf.MooeeInv(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<"=============================================================="<<std::endl;
+  std::cout<<"= Test MeeInvDag MeeDag = 1                                   "<<std::endl;
+  std::cout<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.MooeeDag(chi_e,src_e);
+  Ddwf.MooeeInvDag(src_e,phi_e);
+
+  Ddwf.MooeeDag(chi_o,src_o);
+  Ddwf.MooeeInvDag(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<"=============================================================="<<std::endl;
+  std::cout<<"= Test MpcDagMpc is Hermitian              "<<std::endl;
+  std::cout<<"=============================================================="<<std::endl;
+  
+  random(*RNG5,phi);
+  random(*RNG5,chi);
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+  RealD t1,t2;
+
+  Ddwf.MpcDagMpc(chi_e,dchi_e,t1,t2);
+  Ddwf.MpcDagMpc(chi_o,dchi_o,t1,t2);
+
+  Ddwf.MpcDagMpc(phi_e,dphi_e,t1,t2);
+  Ddwf.MpcDagMpc(phi_o,dphi_o,t1,t2);
+
+  pDce = innerProduct(phi_e,dchi_e);
+  pDco = innerProduct(phi_o,dchi_o);
+  cDpe = innerProduct(chi_e,dphi_e);
+  cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
+  std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
+  
+}