mirror of
				https://github.com/paboyle/Grid.git
				synced 2025-11-04 14:04:32 +00:00 
			
		
		
		
	Jackson smoothed chebyshev and (untested) completion of force terms
for Cayley, Partial and Cont fraction dwf and overlap. have even odd and unprec forces.
This commit is contained in:
		@@ -59,6 +59,8 @@ 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) ){
 | 
					    Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
 | 
				
			||||||
      lo=_lo;
 | 
					      lo=_lo;
 | 
				
			||||||
      hi=_hi;
 | 
					      hi=_hi;
 | 
				
			||||||
@@ -77,7 +79,34 @@ namespace Grid {
 | 
				
			|||||||
	Coeffs[j] = s * 2.0/order;
 | 
						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 approx(double x) // Convenience for plotting the approximation
 | 
				
			||||||
    {
 | 
					    {
 | 
				
			||||||
      double Tn;
 | 
					      double Tn;
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -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
 | 
					    // Assemble Din
 | 
				
			||||||
    for(int s=0;s<Ls;s++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
      if ( s==0 ) {
 | 
					      if ( s==0 ) {
 | 
				
			||||||
@@ -37,11 +34,52 @@ namespace QCD {
 | 
				
			|||||||
	axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
 | 
						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);
 | 
					    DW(Din,chi,DaggerNo);
 | 
				
			||||||
    // ((b D_W + D_w hop terms +1) on s-diag
 | 
					    // ((b D_W + D_w hop terms +1) on s-diag
 | 
				
			||||||
    axpby(chi,1.0,1.0,chi,psi); 
 | 
					    axpby(chi,1.0,1.0,chi,psi); 
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					    // Call Mooee??
 | 
				
			||||||
    for(int s=0;s<Ls;s++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
      if ( s==0 ){
 | 
					      if ( s==0 ){
 | 
				
			||||||
	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
 | 
						axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
 | 
				
			||||||
@@ -67,10 +105,14 @@ namespace QCD {
 | 
				
			|||||||
    // Apply Dw
 | 
					    // Apply Dw
 | 
				
			||||||
    DW(psi,Din,DaggerYes); 
 | 
					    DW(psi,Din,DaggerYes); 
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					    Meooe5D(Din,chi);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    for(int s=0;s<Ls;s++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
 | 
					
 | 
				
			||||||
      // Collect the terms in DW
 | 
					      // Collect the terms in DW
 | 
				
			||||||
      //	Chi = bs Din[s] + cs[s] Din[s+1}
 | 
					      //	Chi = bs Din[s] + cs[s] Din[s+1}
 | 
				
			||||||
      //    Chi+= -mass*cs[s] psi[s+1}
 | 
					      //    Chi+= -mass*cs[s] psi[s+1}
 | 
				
			||||||
 | 
					      /*
 | 
				
			||||||
      if ( s==0 ) {
 | 
					      if ( s==0 ) {
 | 
				
			||||||
	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
 | 
						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);
 | 
						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_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);
 | 
						axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
 | 
				
			||||||
      }
 | 
					      }
 | 
				
			||||||
 | 
					      */
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					      // FIXME just call MooeeDag??
 | 
				
			||||||
 | 
					
 | 
				
			||||||
      // Collect the terms indept of DW
 | 
					      // Collect the terms indept of DW
 | 
				
			||||||
      if ( s==0 ){
 | 
					      if ( s==0 ){
 | 
				
			||||||
	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
 | 
						axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
 | 
				
			||||||
@@ -103,6 +149,9 @@ namespace QCD {
 | 
				
			|||||||
  {
 | 
					  {
 | 
				
			||||||
    LatticeFermion tmp(psi._grid);
 | 
					    LatticeFermion tmp(psi._grid);
 | 
				
			||||||
    // Assemble the 5d matrix
 | 
					    // 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++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
      if ( s==0 ) {
 | 
					      if ( s==0 ) {
 | 
				
			||||||
	//	tmp = bs psi[s] + cs[s] psi[s+1}
 | 
						//	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);
 | 
						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
 | 
					    // Apply 4d dslash
 | 
				
			||||||
    if ( psi.checkerboard == Odd ) {
 | 
					    if ( psi.checkerboard == Odd ) {
 | 
				
			||||||
      DhopEO(tmp,chi,DaggerNo);
 | 
					      DhopEO(tmp,chi,DaggerNo);
 | 
				
			||||||
@@ -134,6 +184,10 @@ namespace QCD {
 | 
				
			|||||||
    } else {
 | 
					    } else {
 | 
				
			||||||
      DhopOE(psi,tmp,DaggerYes);
 | 
					      DhopOE(psi,tmp,DaggerYes);
 | 
				
			||||||
    }
 | 
					    }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					    Meooe5D(tmp,chi); 
 | 
				
			||||||
 | 
					    std::cout << "Meooe Test replacement norm2 chi new = " <<norm2(chi)<<std::endl;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    // Assemble the 5d matrix
 | 
					    // Assemble the 5d matrix
 | 
				
			||||||
    for(int s=0;s<Ls;s++){
 | 
					    for(int s=0;s<Ls;s++){
 | 
				
			||||||
      if ( s==0 ) {
 | 
					      if ( s==0 ) {
 | 
				
			||||||
@@ -147,6 +201,8 @@ namespace QCD {
 | 
				
			|||||||
	axpby_ssp_pminus(chi,1.0   ,chi,-ceo[s-1],tmp,s,s-1);
 | 
						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)
 | 
					  void CayleyFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
 | 
				
			||||||
@@ -250,6 +306,50 @@ namespace QCD {
 | 
				
			|||||||
    }
 | 
					    }
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					  // 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
 | 
					  // Tanh
 | 
				
			||||||
  void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
 | 
					  void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
 | 
				
			||||||
  {
 | 
					  {
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -22,9 +22,17 @@ namespace Grid {
 | 
				
			|||||||
      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 | 
					      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 | 
				
			||||||
      virtual void   Instantiatable(void)=0;
 | 
					      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
 | 
					      // Efficient support for multigrid coarsening
 | 
				
			||||||
      virtual void  Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp);
 | 
					      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:
 | 
					      //    protected:
 | 
				
			||||||
      RealD mass;
 | 
					      RealD mass;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -170,6 +170,53 @@ namespace Grid {
 | 
				
			|||||||
      MooeeInv(psi,chi);
 | 
					      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
 | 
					    // Constructors
 | 
				
			||||||
    ContinuedFractionFermion5D::ContinuedFractionFermion5D(
 | 
					    ContinuedFractionFermion5D::ContinuedFractionFermion5D(
 | 
				
			||||||
							   LatticeGaugeField &_Umu,
 | 
												   LatticeGaugeField &_Umu,
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -21,6 +21,11 @@ namespace Grid {
 | 
				
			|||||||
      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
 | 
					      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
 | 
				
			||||||
      virtual void   MooeeInvDag (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;
 | 
				
			||||||
      virtual void   Instantiatable(void) =0;
 | 
					      virtual void   Instantiatable(void) =0;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -254,6 +254,51 @@ namespace Grid {
 | 
				
			|||||||
      MooeeInv_internal(in,out,DaggerYes);
 | 
					      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){
 | 
					    void  PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){
 | 
				
			||||||
      SetCoefficientsZolotarev(1.0/scale,zdata);
 | 
					      SetCoefficientsZolotarev(1.0/scale,zdata);
 | 
				
			||||||
    }
 | 
					    }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -28,6 +28,11 @@ namespace Grid {
 | 
				
			|||||||
      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
 | 
					      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
 | 
				
			||||||
      virtual void   MooeeInvDag (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
 | 
					      virtual void   Instantiatable(void) =0; // ensure no make-eee
 | 
				
			||||||
 | 
					
 | 
				
			||||||
      // Efficient support for multigrid coarsening
 | 
					      // Efficient support for multigrid coarsening
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -124,11 +124,14 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st,
 | 
				
			|||||||
 | 
					
 | 
				
			||||||
  WilsonCompressor compressor(dag);
 | 
					  WilsonCompressor compressor(dag);
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
  LatticeColourMatrix tmp(B._grid);
 | 
					  LatticeColourMatrix tmp(mat._grid);
 | 
				
			||||||
  LatticeFermion Btilde(B._grid);
 | 
					  LatticeFermion Btilde(B._grid);
 | 
				
			||||||
 | 
					  LatticeFermion Atilde(B._grid);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
 | 
					  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(B,comm_buf,compressor);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 | 
					  Atilde=A;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  for(int mu=0;mu<Nd;mu++){
 | 
					  for(int mu=0;mu<Nd;mu++){
 | 
				
			||||||
      
 | 
					      
 | 
				
			||||||
    ////////////////////////////////////////////////////////////////////////
 | 
					    ////////////////////////////////////////////////////////////////////////
 | 
				
			||||||
@@ -140,20 +143,28 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st,
 | 
				
			|||||||
    ////////////////////////
 | 
					    ////////////////////////
 | 
				
			||||||
    // Call the single hop
 | 
					    // Call the single hop
 | 
				
			||||||
    ////////////////////////
 | 
					    ////////////////////////
 | 
				
			||||||
 | 
					    tmp = zero;
 | 
				
			||||||
 | 
					
 | 
				
			||||||
PARALLEL_FOR_LOOP
 | 
					PARALLEL_FOR_LOOP
 | 
				
			||||||
    for(int sss=0;sss<B._grid->oSites();sss++){
 | 
					    for(int sss=0;sss<B._grid->oSites();sss++){
 | 
				
			||||||
      for(int s=0;s<Ls;s++){
 | 
					      for(int s=0;s<Ls;s++){
 | 
				
			||||||
	int sU=sss;
 | 
						int sU=sss;
 | 
				
			||||||
	int sF = s+Ls*sU;
 | 
						int sF = s+Ls*sU;
 | 
				
			||||||
	DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma);
 | 
						DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma);
 | 
				
			||||||
      }
 | 
					 | 
				
			||||||
    }
 | 
					 | 
				
			||||||
 | 
					
 | 
				
			||||||
    ////////////////////////////
 | 
					    ////////////////////////////
 | 
				
			||||||
    // spin trace outer product
 | 
					    // 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);
 | 
					    PokeIndex<LorentzIndex>(mat,tmp,mu);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  }
 | 
					  }
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -1,6 +1,17 @@
 | 
				
			|||||||
#!/bin/bash
 | 
					#!/bin/bash
 | 
				
			||||||
 | 
					
 | 
				
			||||||
export LANG=C
 | 
					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
 | 
				
			||||||
							
								
								
									
										6
									
								
								scripts/substitute-global
									
									
									
									
									
										Executable file
									
								
							
							
						
						
									
										6
									
								
								scripts/substitute-global
									
									
									
									
									
										Executable file
									
								
							@@ -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 {} \;
 | 
				
			||||||
@@ -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
 | 
					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_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_SOURCES=Test_contfrac_cg.cc
 | 
				
			||||||
Test_contfrac_cg_LDADD=-lGrid
 | 
					Test_contfrac_cg_LDADD=-lGrid
 | 
				
			||||||
 | 
					
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -285,6 +285,8 @@ public:
 | 
				
			|||||||
    //    Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
 | 
					    //    Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
 | 
				
			||||||
    Chebyshev<FineField> Cheby    (2.0,70.0,10,InverseApproximation);
 | 
					    Chebyshev<FineField> Cheby    (2.0,70.0,10,InverseApproximation);
 | 
				
			||||||
    Chebyshev<FineField> ChebyAccu(2.0,70.0,10,InverseApproximation);
 | 
					    Chebyshev<FineField> ChebyAccu(2.0,70.0,10,InverseApproximation);
 | 
				
			||||||
 | 
					    Cheby.JacksonSmooth();
 | 
				
			||||||
 | 
					    ChebyAccu.JacksonSmooth();
 | 
				
			||||||
 | 
					
 | 
				
			||||||
    _Aggregates.ProjectToSubspace  (Csrc,in);
 | 
					    _Aggregates.ProjectToSubspace  (Csrc,in);
 | 
				
			||||||
    _Aggregates.PromoteFromSubspace(Csrc,out);
 | 
					    _Aggregates.PromoteFromSubspace(Csrc,out);
 | 
				
			||||||
@@ -385,7 +387,7 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
  LatticeGaugeField Umu(UGrid); 
 | 
					  LatticeGaugeField Umu(UGrid); 
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  NerscField header;
 | 
					  NerscField header;
 | 
				
			||||||
  std::string file("./ckpoint_lat.400");
 | 
					  std::string file("./ckpoint_lat.4000");
 | 
				
			||||||
  readNerscConfiguration(Umu,header,file);
 | 
					  readNerscConfiguration(Umu,header,file);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  //  SU3::ColdConfiguration(RNG4,Umu);
 | 
					  //  SU3::ColdConfiguration(RNG4,Umu);
 | 
				
			||||||
 
 | 
				
			|||||||
@@ -24,13 +24,14 @@ int main (int argc, char ** argv)
 | 
				
			|||||||
  // simplify template declaration? Strip the lorentz from the second template
 | 
					  // simplify template declaration? Strip the lorentz from the second template
 | 
				
			||||||
  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
 | 
					  WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(5.6);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  Real mass=0.01;
 | 
					  Real mass=-0.77;
 | 
				
			||||||
  WilsonFermion FermOp(U,Fine,RBFine,mass);
 | 
					  WilsonFermion FermOp(U,Fine,RBFine,mass);
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
 | 
					  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
 | 
				
			||||||
 | 
					  ConjugateGradient<LatticeFermion>  CGmd(1.0e-6,10000);
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
  TwoFlavourPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion> 
 | 
					  TwoFlavourPseudoFermionAction<LatticeLorentzColourMatrix, LatticeColourMatrix,LatticeFermion> 
 | 
				
			||||||
    Pseudofermion(FermOp,CG,CG,Fine);
 | 
					    Pseudofermion(FermOp,CGmd,CG,Fine);
 | 
				
			||||||
 | 
					
 | 
				
			||||||
  
 | 
					  
 | 
				
			||||||
  //Collect actions
 | 
					  //Collect actions
 | 
				
			||||||
 
 | 
				
			|||||||
		Reference in New Issue
	
	Block a user