Jackson smoothed chebyshev and (untested) completion of force terms

for Cayley, Partial and Cont fraction dwf and overlap.
have even odd and unprec forces.

lib/algorithms/approx/Chebyshev.h

@@ -58,7 +58,9 @@ namespace Grid {
       }
     };
 
-
+    
+    // c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation".
+    //
     Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
       lo=_lo;
       hi=_hi;
@@ -77,7 +79,34 @@ namespace Grid {
 	Coeffs[j] = s * 2.0/order;
       }
     };
+    void JacksonSmooth(void){
+      double M=order;
+      double alpha = M_PI/(M+2);
+      double lmax = std::cos(alpha);
+      double sumUsq =0;
+      std::vector<double> U(M);
+      std::vector<double> a(M);
+      std::vector<double> g(M);
+      for(int n=0;n<=M;n++){
+	U[n] = std::sin((n+1)*std::acos(lmax))/std::sin(std::acos(lmax));
+	sumUsq += U[n]*U[n];
+      }      
+      sumUsq = std::sqrt(sumUsq);
 
+      for(int i=1;i<=M;i++){
+	a[i] = U[i]/sumUsq;
+      }
+      g[0] = 1.0;
+      for(int m=1;m<=M;m++){
+	g[m] = 0;
+	for(int i=0;i<=M-m;i++){
+	  g[m]+= a[i]*a[m+i];
+	}
+      }
+      for(int m=1;m<=M;m++){
+	Coeffs[m]*=g[m];
+      }
+    }
     double approx(double x) // Convenience for plotting the approximation
     {
       double Tn;

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -17,11 +17,8 @@ namespace QCD {
  {
  }
 
-  // override multiply
+  void CayleyFermion5D::Meooe5D    (const LatticeFermion &psi, LatticeFermion &Din)
-  RealD CayleyFermion5D::M    (const LatticeFermion &psi, LatticeFermion &chi)
   {
-    LatticeFermion Din(psi._grid);
-
     // Assemble Din
     for(int s=0;s<Ls;s++){
       if ( s==0 ) {
@@ -37,11 +34,52 @@ namespace QCD {
 	axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
       }
     }
+  }
+  void CayleyFermion5D::MeooeDag5D    (const LatticeFermion &psi, LatticeFermion &Din)
+  {
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	axpby_ssp_pplus (Din,bs[s],psi,cs[s+1],psi,s,s+1);
+	axpby_ssp_pminus(Din,1.0,Din,-mass*cs[Ls-1],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus (Din,bs[s],psi,-mass*cs[0],psi,s,0);
+	axpby_ssp_pminus(Din,1.0,Din,cs[s-1],psi,s,s-1);
+      } else {
+	axpby_ssp_pplus (Din,bs[s],psi,cs[s+1],psi,s,s+1);
+	axpby_ssp_pminus(Din,1.0,Din,cs[s-1],psi,s,s-1);
+      }
+    }
+  }
+
+  // override multiply
+  RealD CayleyFermion5D::M    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion Din(psi._grid);
+
+    // Assemble Din
+    /*
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	Din = bs psi[s] + cs[s] psi[s+1}
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	//      Din+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pplus (Din,1.0,Din,-mass*cs[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(Din,bs[s],psi,-mass*cs[s],psi,s,0);
+	axpby_ssp_pplus (Din,1.0,Din,cs[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
+      }
+    }
+    */
+    Meooe5D(psi,Din);
 
     DW(Din,chi,DaggerNo);
     // ((b D_W + D_w hop terms +1) on s-diag
     axpby(chi,1.0,1.0,chi,psi); 
 
+    // Call Mooee??
     for(int s=0;s<Ls;s++){
       if ( s==0 ){
 	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
@@ -67,10 +105,14 @@ namespace QCD {
     // Apply Dw
     DW(psi,Din,DaggerYes); 
 
+    Meooe5D(Din,chi);
+
     for(int s=0;s<Ls;s++){
+
       // Collect the terms in DW
       //	Chi = bs Din[s] + cs[s] Din[s+1}
       //    Chi+= -mass*cs[s] psi[s+1}
+      /*
       if ( s==0 ) {
 	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
 	axpby_ssp_pminus(chi,1.0,chi,-mass*cs[Ls-1],Din,s,Ls-1);
@@ -81,6 +123,10 @@ namespace QCD {
 	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
 	axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
       }
+      */
+
+      // FIXME just call MooeeDag??
+
       // Collect the terms indept of DW
       if ( s==0 ){
 	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
@@ -103,6 +149,9 @@ namespace QCD {
   {
     LatticeFermion tmp(psi._grid);
     // Assemble the 5d matrix
+    Meooe5D(psi,tmp); 
+    std::cout << "Meooe Test replacement norm2 tmp = " <<norm2(tmp)<<std::endl;
+
     for(int s=0;s<Ls;s++){
       if ( s==0 ) {
 	//	tmp = bs psi[s] + cs[s] psi[s+1}
@@ -117,6 +166,7 @@ namespace QCD {
 	axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
       }
     }
+    std::cout << "Meooe Test replacement norm2 tmp old = " <<norm2(tmp)<<std::endl;
     // Apply 4d dslash
     if ( psi.checkerboard == Odd ) {
       DhopEO(tmp,chi,DaggerNo);
@@ -134,6 +184,10 @@ namespace QCD {
     } else {
       DhopOE(psi,tmp,DaggerYes);
     }
+
+    Meooe5D(tmp,chi); 
+    std::cout << "Meooe Test replacement norm2 chi new = " <<norm2(chi)<<std::endl;
+
     // Assemble the 5d matrix
     for(int s=0;s<Ls;s++){
       if ( s==0 ) {
@@ -147,6 +201,8 @@ namespace QCD {
 	axpby_ssp_pminus(chi,1.0   ,chi,-ceo[s-1],tmp,s,s-1);
       }
     }
+    std::cout << "Meooe Test replacement norm2 chi old = " <<norm2(chi)<<std::endl;
+
   }
 
   void CayleyFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
@@ -249,6 +305,50 @@ namespace QCD {
       axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1);  // chi[Ls]
     }
   }
+
+  // force terms; five routines; default to Dhop on diagonal
+  void CayleyFermion5D::MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion Din(V._grid);
+
+    if ( dag == DaggerNo ) {
+      //      U d/du [D_w D5] V = U d/du DW D5 V
+      Meooe5D(V,Din);
+      DhopDeriv(mat,U,Din,dag);
+    } else {
+      //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
+      Meooe5D(U,Din);
+      DhopDeriv(mat,Din,V,dag);
+    }
+  };
+  void CayleyFermion5D::MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion Din(V._grid);
+
+    if ( dag == DaggerNo ) {
+      //      U d/du [D_w D5] V = U d/du DW D5 V
+      Meooe5D(V,Din);
+      DhopDerivOE(mat,U,Din,dag);
+    } else {
+      //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
+      Meooe5D(U,Din);
+      DhopDerivOE(mat,Din,V,dag);
+    }
+  };
+  void CayleyFermion5D::MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion Din(V._grid);
+
+    if ( dag == DaggerNo ) {
+      //      U d/du [D_w D5] V = U d/du DW D5 V
+      Meooe5D(V,Din);
+      DhopDerivEO(mat,U,Din,dag);
+    } else {
+      //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
+      Meooe5D(U,Din);
+      DhopDerivEO(mat,Din,V,dag);
+    }
+  };
   
   // Tanh
   void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -22,8 +22,16 @@ namespace Grid {
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
       virtual void   Instantiatable(void)=0;
 
+      // force terms; five routines; default to Dhop on diagonal
+      virtual void MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+
       // Efficient support for multigrid coarsening
       virtual void  Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
+      
+      void   Meooe5D       (const LatticeFermion &in, LatticeFermion &out);
+      void   MeooeDag5D    (const LatticeFermion &in, LatticeFermion &out);
 
       //    protected:
       RealD mass;

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -169,6 +169,53 @@ namespace Grid {
     {
       MooeeInv(psi,chi);
     }
+
+  // force terms; five routines; default to Dhop on diagonal
+   void ContinuedFractionFermion5D::MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int sign=1;
+    for(int s=0;s<Ls;s++){
+      if ( s==(Ls-1) ){
+	ag5xpby_ssp(D,Beta[s]*ZoloHiInv,U,0.0,U,s,s);
+      } else {
+	ag5xpby_ssp(D,cc[s]*Beta[s]*sign*ZoloHiInv,U,0.0,U,s,s);
+      }
+      sign=-sign; 
+    }
+    DhopDeriv(mat,D,V,DaggerNo); 
+  };
+   void ContinuedFractionFermion5D::MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int sign=1;
+    for(int s=0;s<Ls;s++){
+      if ( s==(Ls-1) ){
+	ag5xpby_ssp(D,Beta[s]*ZoloHiInv,U,0.0,U,s,s);
+      } else {
+	ag5xpby_ssp(D,cc[s]*Beta[s]*sign*ZoloHiInv,U,0.0,U,s,s);
+      }
+      sign=-sign; 
+    }
+    DhopDerivOE(mat,D,V,DaggerNo); 
+  };
+   void ContinuedFractionFermion5D::MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int sign=1;
+    for(int s=0;s<Ls;s++){
+      if ( s==(Ls-1) ){
+	ag5xpby_ssp(D,Beta[s]*ZoloHiInv,U,0.0,U,s,s);
+      } else {
+	ag5xpby_ssp(D,cc[s]*Beta[s]*sign*ZoloHiInv,U,0.0,U,s,s);
+      }
+      sign=-sign; 
+    }
+    DhopDerivEO(mat,D,V,DaggerNo); 
+  };
     
     // Constructors
     ContinuedFractionFermion5D::ContinuedFractionFermion5D(

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -21,6 +21,11 @@ namespace Grid {
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
+      // force terms; five routines; default to Dhop on diagonal
+      virtual void MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+
       //      virtual void   Instantiatable(void)=0;
       virtual void   Instantiatable(void) =0;
 

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -254,6 +254,51 @@ namespace Grid {
       MooeeInv_internal(in,out,DaggerYes);
     }
 
+
+  // force terms; five routines; default to Dhop on diagonal
+   void PartialFractionFermion5D::MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int nblock=(Ls-1)/2;
+    for(int b=0;b<nblock;b++){
+      int s = 2*b;
+      ag5xpby_ssp(D,-scale,U,0.0,U,s,s); 
+      ag5xpby_ssp(D, scale,U,0.0,U,s+1,s+1); 
+    }
+    ag5xpby_ssp(D,p[nblock]*scale/amax,U,0.0,U,Ls-1,Ls-1);
+
+    DhopDeriv(mat,D,V,DaggerNo); 
+  };
+   void PartialFractionFermion5D::MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int nblock=(Ls-1)/2;
+    for(int b=0;b<nblock;b++){
+      int s = 2*b;
+      ag5xpby_ssp(D,-scale,U,0.0,U,s,s); 
+      ag5xpby_ssp(D, scale,U,0.0,U,s+1,s+1); 
+    }
+    ag5xpby_ssp(D,p[nblock]*scale/amax,U,0.0,U,Ls-1,Ls-1);
+
+    DhopDerivOE(mat,D,V,DaggerNo); 
+  };
+   void PartialFractionFermion5D::MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag)
+  {
+    LatticeFermion D(V._grid);
+
+    int nblock=(Ls-1)/2;
+    for(int b=0;b<nblock;b++){
+      int s = 2*b;
+      ag5xpby_ssp(D,-scale,U,0.0,U,s,s); 
+      ag5xpby_ssp(D, scale,U,0.0,U,s+1,s+1); 
+    }
+    ag5xpby_ssp(D,p[nblock]*scale/amax,U,0.0,U,Ls-1,Ls-1);
+
+    DhopDerivEO(mat,D,V,DaggerNo); 
+  };
+
     void  PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){
       SetCoefficientsZolotarev(1.0/scale,zdata);
     }

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -28,6 +28,11 @@ namespace Grid {
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
+      // force terms; five routines; default to Dhop on diagonal
+      virtual void MDeriv  (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+      virtual void MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag);
+
       virtual void   Instantiatable(void) =0; // ensure no make-eee
 
       // Efficient support for multigrid coarsening

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -124,11 +124,14 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st,
 
   WilsonCompressor compressor(dag);
   
-  LatticeColourMatrix tmp(B._grid);
+  LatticeColourMatrix tmp(mat._grid);
   LatticeFermion Btilde(B._grid);
+  LatticeFermion Atilde(B._grid);
 
   st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
-  
+
+  Atilde=A;
+
   for(int mu=0;mu<Nd;mu++){
       
     ////////////////////////////////////////////////////////////////////////
@@ -140,20 +143,28 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st,
     ////////////////////////
     // Call the single hop
     ////////////////////////
+    tmp = zero;
+
 PARALLEL_FOR_LOOP
     for(int sss=0;sss<B._grid->oSites();sss++){
       for(int s=0;s<Ls;s++){
 	int sU=sss;
 	int sF = s+Ls*sU;
 	DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma);
-      }
-    }
 
     ////////////////////////////
     // spin trace outer product
     ////////////////////////////
-// FIXME : need to sum over fifth direction.
+
-    tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); // ordering here
+	tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(
+		   outerProduct(Btilde[sF],Atilde[sF])); // ordering here
+
+      }
+
+
+    }
+
+
     PokeIndex<LorentzIndex>(mat,tmp,mu);
 
   }

scripts/substitute

@@ -1,6 +1,17 @@
 #!/bin/bash
 
 export LANG=C
-find . -name "*.cc" -print -exec sed -e "s/${1}/${2}/g" -i .bak {} \;
+
-find . -name "*.h" -print -exec sed -e "s/${1}/${2}/g" -i .bak {} \;
+WAS=$1
-find . -name "*.bak" -exec rm -f {} \;
+shift
+IS=$1
+shift
+
+while (( "$#" )); do
+
+echo $1
+sed -e "s/${WAS}/${IS}/g" -i .bak $1
+
+shift
+
+done

scripts/substitute-global

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+export LANG=C
+find . -name "*.cc" -print -exec sed -e "s/${1}/${2}/g" -i .bak {} \;
+find . -name "*.h" -print -exec sed -e "s/${1}/${2}/g" -i .bak {} \;
+find . -name "*.bak" -exec rm -f {} \;

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec 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_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update 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_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby 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_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update 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_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@@ -30,6 +30,10 @@ Test_cf_cr_unprec_SOURCES=Test_cf_cr_unprec.cc
 Test_cf_cr_unprec_LDADD=-lGrid
 
 
+Test_cheby_SOURCES=Test_cheby.cc
+Test_cheby_LDADD=-lGrid
+
+
 Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc
 Test_contfrac_cg_LDADD=-lGrid
 

tests/Test_dwf_hdcr.cc

@@ -285,6 +285,8 @@ public:
     //    Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
     Chebyshev<FineField> Cheby    (2.0,70.0,10,InverseApproximation);
     Chebyshev<FineField> ChebyAccu(2.0,70.0,10,InverseApproximation);
+    Cheby.JacksonSmooth();
+    ChebyAccu.JacksonSmooth();
 
     _Aggregates.ProjectToSubspace  (Csrc,in);
     _Aggregates.PromoteFromSubspace(Csrc,out);
@@ -385,7 +387,7 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(UGrid); 
 
   NerscField header;
-  std::string file("./ckpoint_lat.400");
+  std::string file("./ckpoint_lat.4000");
   readNerscConfiguration(Umu,header,file);
 
   //  SU3::ColdConfiguration(RNG4,Umu);

tests/Test_hmc_WilsonFermionGauge.cc

@@ -24,13 +24,14 @@ int main (int argc, char ** argv)
   // simplify template declaration? Strip the lorentz from the second template
   WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
 
-  Real mass=0.01;
+  Real mass=-0.77;
   WilsonFermion FermOp(U,Fine,RBFine,mass);
   
   ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
+  ConjugateGradient<LatticeFermion>  CGmd(1.0e-6,10000);
   
   TwoFlavourPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion> 
-    Pseudofermion(FermOp,CG,CG,Fine);
+    Pseudofermion(FermOp,CGmd,CG,Fine);
 
   
   //Collect actions