Mshift update

Conflicts:
	systems/PVC/benchmarks/run-2tile-mpi.sh
	systems/PVC/config-command

Grid/algorithms/iterative/ConjugateGradientMultiShiftCleanup.h

@@ -0,0 +1,373 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/ConjugateGradientMultiShift.h
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Christopher Kelly <ckelly@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+//CK 2020: A variant of the multi-shift conjugate gradient with the matrix multiplication in single precision. 
+//The residual is stored in single precision, but the search directions and solution are stored in double precision. 
+//Every update_freq iterations the residual is corrected in double precision. 
+//For safety the a final regular CG is applied to clean up if necessary
+
+//PB Pure single, then double fixup
+
+template<class FieldD, class FieldF,
+	 typename std::enable_if< getPrecision<FieldD>::value == 2, int>::type = 0,
+	 typename std::enable_if< getPrecision<FieldF>::value == 1, int>::type = 0> 
+class ConjugateGradientMultiShiftMixedPrecCleanup : public OperatorMultiFunction<FieldD>,
+					     public OperatorFunction<FieldD>
+{
+public:                                                
+
+  using OperatorFunction<FieldD>::operator();
+
+  RealD   Tolerance;
+  Integer MaxIterationsMshift;
+  Integer MaxIterations;
+  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  std::vector<int> IterationsToCompleteShift;  // Iterations for this shift
+  int verbose;
+  MultiShiftFunction shifts;
+  std::vector<RealD> TrueResidualShift;
+
+  int ReliableUpdateFreq; //number of iterations between reliable updates
+
+  GridBase* SinglePrecGrid; //Grid for single-precision fields
+  LinearOperatorBase<FieldF> &Linop_f; //single precision
+
+  ConjugateGradientMultiShiftMixedPrecCleanup(Integer maxit, const MultiShiftFunction &_shifts,
+				       GridBase* _SinglePrecGrid, LinearOperatorBase<FieldF> &_Linop_f,
+				       int _ReliableUpdateFreq) : 
+    MaxIterationsMshift(maxit),  shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq),
+    MaxIterations(20000)
+  { 
+    verbose=1;
+    IterationsToCompleteShift.resize(_shifts.order);
+    TrueResidualShift.resize(_shifts.order);
+  }
+
+  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, FieldD &psi)
+  {
+    GridBase *grid = src.Grid();
+    int nshift = shifts.order;
+    std::vector<FieldD> results(nshift,grid);
+    (*this)(Linop,src,results,psi);
+  }
+  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, std::vector<FieldD> &results, FieldD &psi)
+  {
+    int nshift = shifts.order;
+
+    (*this)(Linop,src,results);
+  
+    psi = shifts.norm*src;
+    for(int i=0;i<nshift;i++){
+      psi = psi + shifts.residues[i]*results[i];
+    }
+
+    return;
+  }
+
+  void operator() (LinearOperatorBase<FieldD> &Linop_d, const FieldD &src_d, std::vector<FieldD> &psi_d)
+  { 
+    GRID_TRACE("ConjugateGradientMultiShiftMixedPrecCleanup");
+    GridBase *DoublePrecGrid = src_d.Grid();
+
+    ////////////////////////////////////////////////////////////////////////
+    // Convenience references to the info stored in "MultiShiftFunction"
+    ////////////////////////////////////////////////////////////////////////
+    int nshift = shifts.order;
+
+    std::vector<RealD> &mass(shifts.poles); // Make references to array in "shifts"
+    std::vector<RealD> &mresidual(shifts.tolerances);
+    std::vector<RealD> alpha(nshift,1.0);
+
+    //Double precision search directions
+    FieldD p_d(DoublePrecGrid);
+    std::vector<FieldF> ps_f (nshift, SinglePrecGrid);// Search directions (single precision)
+    std::vector<FieldF> psi_f(nshift, SinglePrecGrid);// solutions (single precision)
+
+    FieldD tmp_d(DoublePrecGrid);
+    FieldD r_d(DoublePrecGrid);
+    FieldF r_f(SinglePrecGrid);
+    FieldD mmp_d(DoublePrecGrid);
+
+    assert(psi_d.size()==nshift);
+    assert(mass.size()==nshift);
+    assert(mresidual.size()==nshift);
+  
+    // dynamic sized arrays on stack; 2d is a pain with vector
+    RealD  bs[nshift];
+    RealD  rsq[nshift];
+    RealD  rsqf[nshift];
+    RealD  z[nshift][2];
+    int     converged[nshift];
+  
+    const int       primary =0;
+  
+    //Primary shift fields CG iteration
+    RealD a,b,c,d;
+    RealD cp,bp,qq; //prev
+  
+    // Matrix mult fields
+    FieldF p_f(SinglePrecGrid);
+    FieldF mmp_f(SinglePrecGrid);
+
+    // Check lightest mass
+    for(int s=0;s<nshift;s++){
+      assert( mass[s]>= mass[primary] );
+      converged[s]=0;
+    }
+  
+    // Wire guess to zero
+    // Residuals "r" are src
+    // First search direction "p" is also src
+    cp = norm2(src_d);
+
+    // Handle trivial case of zero src.
+    if( cp == 0. ){
+      for(int s=0;s<nshift;s++){
+	psi_d[s] = Zero();
+	psi_f[s] = Zero();
+	IterationsToCompleteShift[s] = 1;
+	TrueResidualShift[s] = 0.;
+      }
+      return;
+    }
+
+    for(int s=0;s<nshift;s++){
+      rsq[s] = cp * mresidual[s] * mresidual[s];
+      rsqf[s] =rsq[s];
+      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrecCleanup: shift "<< s <<" target resid "<<rsq[s]<<std::endl;
+      //      ps_d[s] = src_d;
+      precisionChangeFast(ps_f[s],src_d);
+    }
+    // r and p for primary
+    p_d = src_d; //primary copy --- make this a reference to ps_d to save axpys
+    r_d = p_d;
+    
+    //MdagM+m[0]
+    precisionChangeFast(p_f,p_d);
+    Linop_f.HermOpAndNorm(p_f,mmp_f,d,qq); // mmp = MdagM p        d=real(dot(p, mmp)),  qq=norm2(mmp)
+    precisionChangeFast(tmp_d,mmp_f);
+    Linop_d.HermOpAndNorm(p_d,mmp_d,d,qq); // mmp = MdagM p        d=real(dot(p, mmp)),  qq=norm2(mmp)
+    tmp_d = tmp_d - mmp_d;
+    std::cout << " Testing operators match "<<norm2(mmp_d)<<" f "<<norm2(mmp_f)<<" diff "<< norm2(tmp_d)<<std::endl;
+    //    assert(norm2(tmp_d)< 1.0e-4);
+
+    axpy(mmp_d,mass[0],p_d,mmp_d);
+    RealD rn = norm2(p_d);
+    d += rn*mass[0];
+
+    b = -cp /d;
+  
+    // Set up the various shift variables
+    int       iz=0;
+    z[0][1-iz] = 1.0;
+    z[0][iz]   = 1.0;
+    bs[0]      = b;
+    for(int s=1;s<nshift;s++){
+      z[s][1-iz] = 1.0;
+      z[s][iz]   = 1.0/( 1.0 - b*(mass[s]-mass[0]));
+      bs[s]      = b*z[s][iz]; 
+    }
+  
+    // r += b[0] A.p[0]
+    // c= norm(r)
+    c=axpy_norm(r_d,b,mmp_d,r_d);
+  
+    for(int s=0;s<nshift;s++) {
+      axpby(psi_d[s],0.,-bs[s]*alpha[s],src_d,src_d);
+      precisionChangeFast(psi_f[s],psi_d[s]);
+    }
+  
+    ///////////////////////////////////////
+    // Timers
+    ///////////////////////////////////////
+    GridStopWatch AXPYTimer, ShiftTimer, QRTimer, MatrixTimer, SolverTimer, PrecChangeTimer, CleanupTimer;
+
+    SolverTimer.Start();
+  
+    // Iteration loop
+    int k;
+  
+    for (k=1;k<=MaxIterationsMshift;k++){    
+
+      a = c /cp;
+      AXPYTimer.Start();
+      axpy(p_d,a,p_d,r_d); 
+      AXPYTimer.Stop();
+
+      PrecChangeTimer.Start();
+      precisionChangeFast(r_f, r_d);
+      PrecChangeTimer.Stop();
+
+      AXPYTimer.Start();
+      for(int s=0;s<nshift;s++){
+	if ( ! converged[s] ) { 
+	  if (s==0){
+	    axpy(ps_f[s],a,ps_f[s],r_f);
+	  } else{
+	    RealD as =a *z[s][iz]*bs[s] /(z[s][1-iz]*b);
+	    axpby(ps_f[s],z[s][iz],as,r_f,ps_f[s]);
+	  }
+	}
+      }
+      AXPYTimer.Stop();
+
+      cp=c;
+      PrecChangeTimer.Start();
+      precisionChangeFast(p_f, p_d); //get back single prec search direction for linop
+      PrecChangeTimer.Stop();
+      MatrixTimer.Start();  
+      Linop_f.HermOp(p_f,mmp_f);
+      MatrixTimer.Stop();  
+      PrecChangeTimer.Start();
+      precisionChangeFast(mmp_d, mmp_f); // From Float to Double
+      PrecChangeTimer.Stop();
+
+      d=real(innerProduct(p_d,mmp_d));    
+      axpy(mmp_d,mass[0],p_d,mmp_d);
+      RealD rn = norm2(p_d);
+      d += rn*mass[0];
+    
+      bp=b;
+      b=-cp/d;
+
+      // Toggle the recurrence history
+      bs[0] = b;
+      iz = 1-iz;
+      ShiftTimer.Start();
+      for(int s=1;s<nshift;s++){
+	if((!converged[s])){
+	  RealD z0 = z[s][1-iz];
+	  RealD z1 = z[s][iz];
+	  z[s][iz] = z0*z1*bp
+	    / (b*a*(z1-z0) + z1*bp*(1- (mass[s]-mass[0])*b)); 
+	  bs[s] = b*z[s][iz]/z0; // NB sign  rel to Mike
+	}
+      }
+      ShiftTimer.Stop();
+
+      //Update single precision solutions
+      AXPYTimer.Start();
+      for(int s=0;s<nshift;s++){
+	int ss = s;
+	if( (!converged[s]) ) { 
+	  axpy(psi_f[ss],-bs[s]*alpha[s],ps_f[s],psi_f[ss]);
+	}
+      }
+      c = axpy_norm(r_d,b,mmp_d,r_d);
+      AXPYTimer.Stop();
+    
+      // Convergence checks
+      int all_converged = 1;
+      for(int s=0;s<nshift;s++){
+      
+	if ( (!converged[s]) ){
+	  IterationsToCompleteShift[s] = k;
+	
+	  RealD css  = c * z[s][iz]* z[s][iz];
+	
+	  if(css<rsqf[s]){
+	    if ( ! converged[s] )
+	      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrecCleanup k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
+	    converged[s]=1;
+	  } else {
+	    all_converged=0;
+	  }
+
+	}
+      }
+
+      if ( all_converged || k == MaxIterationsMshift-1){
+
+	SolverTimer.Stop();
+
+	for(int s=0;s<nshift;s++){
+	  precisionChangeFast(psi_d[s],psi_f[s]);
+	}
+
+	
+	if ( all_converged ){
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrecCleanup: All shifts have converged iteration "<<k<<std::endl;
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrecCleanup: Checking solutions"<<std::endl;
+	} else {
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrecCleanup: Not all shifts have converged iteration "<<k<<std::endl;
+	}
+	
+	// Check answers 
+	for(int s=0; s < nshift; s++) { 
+	  Linop_d.HermOpAndNorm(psi_d[s],mmp_d,d,qq);
+	  axpy(tmp_d,mass[s],psi_d[s],mmp_d);
+	  axpy(r_d,-alpha[s],src_d,tmp_d);
+	  RealD rn = norm2(r_d);
+	  RealD cn = norm2(src_d);
+	  TrueResidualShift[s] = std::sqrt(rn/cn);
+	  std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrecCleanup: shift["<<s<<"] true residual "<< TrueResidualShift[s] << " target " << mresidual[s] << std::endl;
+
+	  //If we have not reached the desired tolerance, do a (mixed precision) CG cleanup
+	  if(rn >= rsq[s]){
+	    CleanupTimer.Start();
+	    std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrecCleanup: performing cleanup step for shift " << s << std::endl;
+
+	    //Setup linear operators for final cleanup
+	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldD> Linop_shift_d(Linop_d, mass[s]);
+	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldF> Linop_shift_f(Linop_f, mass[s]);
+					       
+	    MixedPrecisionConjugateGradient<FieldD,FieldF> cg(mresidual[s], MaxIterations, MaxIterations, SinglePrecGrid, Linop_shift_f, Linop_shift_d); 
+	    cg(src_d, psi_d[s]);
+	    
+	    TrueResidualShift[s] = cg.TrueResidual;
+	    CleanupTimer.Stop();
+	  }
+	}
+
+	std::cout << GridLogMessage << "ConjugateGradientMultiShiftMixedPrecCleanup: Time Breakdown for body"<<std::endl;
+	std::cout << GridLogMessage << "\tSolver    " << SolverTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tAXPY    " << AXPYTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tMatrix    " << MatrixTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tShift    " << ShiftTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\t\tPrecision Change " << PrecChangeTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\tFinal Cleanup " << CleanupTimer.Elapsed()     <<std::endl;
+	std::cout << GridLogMessage << "\tSolver+Cleanup " << SolverTimer.Elapsed() + CleanupTimer.Elapsed() << std::endl;
+
+	IterationsToComplete = k;	
+
+	return;
+      }
+   
+    }
+    std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;
+    assert(0);
+  }
+
+};
+NAMESPACE_END(Grid);
+

Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h

@@ -81,6 +81,7 @@ public:
   using OperatorFunction<FieldD>::operator();
 
   RealD   Tolerance;
+  Integer MaxIterationsMshift;
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
   std::vector<int> IterationsToCompleteShift;  // Iterations for this shift
@@ -95,9 +96,9 @@ public:
 
   ConjugateGradientMultiShiftMixedPrec(Integer maxit, const MultiShiftFunction &_shifts,
 				       GridBase* _SinglePrecGrid, LinearOperatorBase<FieldF> &_Linop_f,
-				       int _ReliableUpdateFreq
-				       ) : 
-    MaxIterations(maxit),  shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq)
+				       int _ReliableUpdateFreq) : 
+    MaxIterationsMshift(maxit),  shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq),
+    MaxIterations(20000)
   { 
     verbose=1;
     IterationsToCompleteShift.resize(_shifts.order);
@@ -247,7 +248,7 @@ public:
     // Iteration loop
     int k;
   
-    for (k=1;k<=MaxIterations;k++){    
+    for (k=1;k<=MaxIterationsMshift;k++){    
 
       a = c /cp;
       AXPYTimer.Start();
@@ -353,12 +354,17 @@ public:
 	}
       }
 
-      if ( all_converged ){
+      if ( all_converged || k == MaxIterationsMshift-1){
 
 	SolverTimer.Stop();
-	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: All shifts have converged iteration "<<k<<std::endl;
-	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: Checking solutions"<<std::endl;
-      
+
+	if ( all_converged ){
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: All shifts have converged iteration "<<k<<std::endl;
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: Checking solutions"<<std::endl;
+	} else {
+	  std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: Not all shifts have converged iteration "<<k<<std::endl;
+	}
+	
 	// Check answers 
 	for(int s=0; s < nshift; s++) { 
 	  Linop_d.HermOpAndNorm(psi_d[s],mmp_d,d,qq);
@@ -399,12 +405,10 @@ public:
 
 	return;
       }
-
    
     }
-    // ugly hack
     std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;
-    //  assert(0);
+    assert(0);
   }
 
 };

Grid/qcd/action/fermion/WilsonCompressor.h

@@ -36,7 +36,7 @@ NAMESPACE_BEGIN(Grid);
 // Wilson compressor will need FaceGather policies for:
 // Periodic, Dirichlet, and partial Dirichlet for DWF
 ///////////////////////////////////////////////////////////////
-const int dwf_compressor_depth=1;
+const int dwf_compressor_depth=2;
 #define DWF_COMPRESS
 class FaceGatherPartialDWF
 {

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h

@@ -127,6 +127,8 @@ NAMESPACE_BEGIN(Grid);
 	  ApproxNegPowerAction.tolerances[i]    = action_tolerance[i];
 	  ApproxHalfPowerAction.tolerances[i]   = action_tolerance[i];
 	  ApproxNegHalfPowerAction.tolerances[i]= action_tolerance[i];
+	}
+	for(int i=0;i<ApproxPowerMD.tolerances.size();i++){
 	  ApproxPowerMD.tolerances[i]       = md_tolerance[i];
 	  ApproxNegPowerMD.tolerances[i]    = md_tolerance[i];
 	  ApproxHalfPowerMD.tolerances[i]   = md_tolerance[i];

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h

@@ -29,6 +29,8 @@
 #ifndef QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
 #define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
 
+#include <Grid/algorithms/iterative/ConjugateGradientMultiShiftCleanup.h>
+
 NAMESPACE_BEGIN(Grid);
 
     /////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -58,7 +60,7 @@ NAMESPACE_BEGIN(Grid);
       //Allow derived classes to override the multishift CG
       virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){
 #if 0
-	SchurDifferentiableOperator<ImplD> schurOp(numerator ? NumOp : DenOp);
+	SchurDifferentiableOperator<ImplD> schurOp(numerator ? NumOpD : DenOpD);
 	ConjugateGradientMultiShift<FermionFieldD> msCG(MaxIter, approx);
 	msCG(schurOp,in, out);
 #else
@@ -66,7 +68,8 @@ NAMESPACE_BEGIN(Grid);
 	SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
 	FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
 	FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
-	
+
+	// Action better with higher precision?
 	ConjugateGradientMultiShiftMixedPrec<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
 	precisionChange(inD2,in);
 	std::cout << "msCG single solve "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
@@ -76,12 +79,12 @@ NAMESPACE_BEGIN(Grid);
       }
       virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &out){
 	SchurDifferentiableOperator<ImplD2> schurOpD2(numerator ? NumOpD2 : DenOpD2);
-	SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
+	SchurDifferentiableOperator<ImplF>  schurOpF (numerator ? NumOpF  : DenOpF);
 
 	FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
 	FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
 	std::vector<FermionFieldD2> out_elemsD2(out_elems.size(),NumOpD2.FermionRedBlackGrid());
-	ConjugateGradientMultiShiftMixedPrec<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
+	ConjugateGradientMultiShiftMixedPrecCleanup<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
 	precisionChange(inD2,in);
 	std::cout << "msCG in "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
 	msCG(schurOpD2, inD2, out_elemsD2, outD2);

Grid/qcd/hmc/integrators/Integrator.h

@@ -284,6 +284,15 @@ public:
 		  << as[level].actions.at(actionID)->deriv_us*1.0e-6<<" s"<< std::endl;
       }
     }
+    std::cout << GridLogMessage << "--------------------------- "<<std::endl;
+    std::cout << GridLogMessage << " Dslash counts "<<std::endl;
+    std::cout << GridLogMessage << "------------------------- "<<std::endl;
+    uint64_t full, partial, dirichlet;
+    DslashGetCounts(dirichlet,partial,full);
+    std::cout << GridLogMessage << " Full BCs               : "<<full<<std::endl;
+    std::cout << GridLogMessage << " Partial dirichlet BCs  : "<<partial<<std::endl;
+    std::cout << GridLogMessage << " Dirichlet BCs          : "<<dirichlet<<std::endl;
+
     std::cout << GridLogMessage << "--------------------------- "<<std::endl;
     std::cout << GridLogMessage << " Force average size "<<std::endl;
     std::cout << GridLogMessage << "------------------------- "<<std::endl;

Grid/stencil/Stencil.cc

@@ -29,6 +29,27 @@
 
 NAMESPACE_BEGIN(Grid);
 
+uint64_t DslashFullCount;
+uint64_t DslashPartialCount;
+uint64_t DslashDirichletCount;
+
+void DslashResetCounts(void)
+{
+  DslashFullCount=0;
+  DslashPartialCount=0;
+  DslashDirichletCount=0;
+}
+void DslashGetCounts(uint64_t &dirichlet,uint64_t &partial,uint64_t &full)
+{
+  dirichlet = DslashDirichletCount;
+  partial   = DslashPartialCount;
+  full      = DslashFullCount;
+}
+void DslashLogFull(void)     { DslashFullCount++;}
+void DslashLogPartial(void)  { DslashPartialCount++;}
+void DslashLogDirichlet(void){ DslashDirichletCount++;}
+
+
 void Gather_plane_table_compute (GridBase *grid,int dimension,int plane,int cbmask,
 				 int off,std::vector<std::pair<int,int> > & table)
 {

Grid/stencil/Stencil.h

@@ -120,6 +120,12 @@ void Gather_plane_exchange_table(commVector<std::pair<int,int> >& table,
 }
 */
 
+void DslashResetCounts(void);
+void DslashGetCounts(uint64_t &dirichlet,uint64_t &partial,uint64_t &full);
+void DslashLogFull(void);
+void DslashLogPartial(void);
+void DslashLogDirichlet(void);
+
 struct StencilEntry {
 #ifdef GRID_CUDA
   uint64_t _byte_offset;       // 8 bytes
@@ -312,6 +318,7 @@ public:
 
   int face_table_computed;
   int partialDirichlet;
+  int fullDirichlet;
   std::vector<commVector<std::pair<int,int> > > face_table ;
   Vector<int> surface_list;
 
@@ -442,6 +449,9 @@ public:
   void CommunicateComplete(std::vector<std::vector<CommsRequest_t> > &reqs)
   {
     _grid->StencilSendToRecvFromComplete(MpiReqs,0);
+    if   ( this->partialDirichlet ) DslashLogPartial();
+    else if ( this->fullDirichlet ) DslashLogDirichlet();
+    else DslashLogFull();
   }
   ////////////////////////////////////////////////////////////////////////
   // Blocking send and receive. Either sequential or parallel.
@@ -770,6 +780,10 @@ public:
     if ( p.dirichlet.size() ==0 ) p.dirichlet.resize(grid->Nd(),0);
     partialDirichlet = p.partialDirichlet;
     DirichletBlock(p.dirichlet); // comms send/recv set up
+    fullDirichlet=0;
+    for(int d=0;d<p.dirichlet.size();d++){
+      if (p.dirichlet[d]) fullDirichlet=1;
+    }
 
     _unified_buffer_size=0;
     surface_list.resize(0);

HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc

@@ -164,11 +164,6 @@ int main(int argc, char **argv) {
   typedef MobiusEOFAFermionF FermionEOFAActionF;
   typedef typename FermionActionF::FermionField FermionFieldF;
 
-  typedef WilsonImplD2 FermionImplPolicyD2;
-  typedef MobiusFermionD2 FermionActionD2;
-  typedef MobiusEOFAFermionD2 FermionEOFAActionD2;
-  typedef typename FermionActionD2::FermionField FermionFieldD2;
-
   typedef Grid::XmlReader       Serialiser;
 
   //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -232,31 +227,34 @@ int main(int argc, char **argv) {
   //  std::vector<Real> hasenbusch({ light_mass, 0.005, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); // Updated
   //  std::vector<Real> hasenbusch({ light_mass, 0.0145, 0.045, 0.108, 0.25, 0.51 , 0.75 , pv_mass });
 
-  OneFlavourRationalParams OFRp; // Up/down
-  OFRp.lo       = 4.0e-5;
+  int SP_iters=10000;
+  
+  RationalActionParams OFRp; // Up/down
+  OFRp.lo       = 6.0e-5;
   OFRp.hi       = 90.0;
-  OFRp.MaxIter  = 60000;
-  OFRp.tolerance= 1.0e-5;
-  OFRp.mdtolerance= 1.0e-3;
+  OFRp.inv_pow  = 2;
+  OFRp.MaxIter  = SP_iters; // get most shifts by 2000, stop sharing space
+  OFRp.action_tolerance= 1.0e-8;
+  OFRp.action_degree   = 18;
+  OFRp.md_tolerance= 1.0e-5;
+  OFRp.md_degree   = 14;
   //  OFRp.degree   = 20; converges
   //  OFRp.degree   = 16;
-  OFRp.degree   = 18;
   OFRp.precision= 80;
   OFRp.BoundsCheckFreq=0;
   std::vector<RealD> ActionTolByPole({
-      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-7,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8
     });
   std::vector<RealD> MDTolByPole({
-      1.0e-5,5.0e-6,1.0e-6,1.0e-7, // soften convergence more more
+      1.6e-5,5.0e-6,1.0e-6,3.0e-7, // soften convergence more more
       //      1.0e-6,3.0e-7,1.0e-7,1.0e-7,
       //      3.0e-6,1.0e-6,1.0e-7,1.0e-7, // soften convergence
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
-      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8
     });
 
@@ -265,10 +263,8 @@ int main(int argc, char **argv) {
 
   typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
   typedef SchurDiagMooeeOperator<FermionAction ,FermionField > LinearOperatorD;
-  typedef SchurDiagMooeeOperator<FermionActionD2,FermionFieldD2 > LinearOperatorD2;
   typedef SchurDiagMooeeOperator<FermionEOFAActionF,FermionFieldF> LinearOperatorEOFAF;
   typedef SchurDiagMooeeOperator<FermionEOFAAction ,FermionField > LinearOperatorEOFAD;
-  typedef SchurDiagMooeeOperator<FermionEOFAActionD2,FermionFieldD2 > LinearOperatorEOFAD2;
   typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusFermionD,MobiusFermionF,LinearOperatorD,LinearOperatorF> MxPCG;
   typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusEOFAFermionD,MobiusEOFAFermionF,LinearOperatorEOFAD,LinearOperatorEOFAF> MxPCG_EOFA;
 
@@ -321,7 +317,6 @@ int main(int argc, char **argv) {
   // temporarily need a gauge field
   LatticeGaugeFieldD  U(GridPtr); U=Zero();
   LatticeGaugeFieldF  UF(GridPtrF); UF=Zero();
-  LatticeGaugeFieldD2 UD2(GridPtrF); UD2=Zero();
 
   std::cout << GridLogMessage << " Running the HMC "<< std::endl;
   TheHMC.ReadCommandLine(argc,argv);  // params on CML or from param file
@@ -340,6 +335,7 @@ int main(int argc, char **argv) {
   ParamsDirF.dirichlet=Dirichlet;
   ParamsDir.partialDirichlet=1;
   ParamsDirF.partialDirichlet=1;
+  std::cout << GridLogMessage<< "Partial Dirichlet depth is "<<dwf_compressor_depth<<std::endl;
 
   //  double StoppingCondition = 1e-14;
   //  double MDStoppingCondition = 1e-9;
@@ -424,7 +420,7 @@ int main(int argc, char **argv) {
 	 ActionCGL, ActionCGR,
 	 DerivativeCGL, DerivativeCGR,
 	 SFRp, true);
-  //  Level2.push_back(&EOFA);
+  Level2.push_back(&EOFA);
 
   ////////////////////////////////////
   // up down action
@@ -449,17 +445,15 @@ int main(int argc, char **argv) {
   std::vector<FermionAction *> Denominators;
   std::vector<FermionActionF *> NumeratorsF;
   std::vector<FermionActionF *> DenominatorsF;
-  std::vector<FermionActionD2 *> NumeratorsD2;
-  std::vector<FermionActionD2 *> DenominatorsD2;
   std::vector<TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy> *> Quotients;
   std::vector<MxPCG *> ActionMPCG;
   std::vector<MxPCG *> MPCG;
   
 #define MIXED_PRECISION
 #ifdef MIXED_PRECISION
-  std::vector<OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicyD2> *> Bdys;
+  std::vector<GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy> *> Bdys;
 #else
-  std::vector<OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> *> Bdys;
+  std::vector<GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> *> Bdys;
 #endif
 
   typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
@@ -532,31 +526,19 @@ int main(int argc, char **argv) {
       Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],CG));
     } else {
 #ifdef MIXED_PRECISION
-      // Use the D2 data types and make them use same grid as single
-      FermionActionD2::ImplParams ParamsDenD2(boundary);
-      FermionActionD2::ImplParams ParamsNumD2(boundary);
-
-      ParamsDenD2.dirichlet = ParamsDen.dirichlet;
-      ParamsDenD2.partialDirichlet = ParamsDen.partialDirichlet;
-      DenominatorsD2.push_back(new FermionActionD2(UD2,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_den[h],M5,b,c, ParamsDenD2));
-
-      ParamsNumD2.dirichlet = ParamsNum.dirichlet;
-      ParamsNumD2.partialDirichlet = ParamsNum.partialDirichlet;
-      NumeratorsD2.push_back  (new FermionActionD2(UD2,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_num[h],M5,b,c, ParamsNumD2));
-    
-      Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicyD2>(
+      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy>(
 			   *Numerators[h],*Denominators[h],
 			   *NumeratorsF[h],*DenominatorsF[h],
-			   *NumeratorsD2[h],*DenominatorsD2[h],
-			   OFRp, 400) );
-      Bdys.push_back( new OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicyD2>(
+			   *Numerators[h],*Denominators[h],
+			   OFRp, SP_iters) );
+      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy>(
 			   *Numerators[h],*Denominators[h],
 			   *NumeratorsF[h],*DenominatorsF[h],
-			   *NumeratorsD2[h],*DenominatorsD2[h],
-			   OFRp, 400) );
+			   *Numerators[h],*Denominators[h],
+			   OFRp, SP_iters) );
 #else
-      Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
-      Bdys.push_back( new OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
+      Bdys.push_back( new GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
+      Bdys.push_back( new GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));
 #endif
     }
   }

HMC/Mobius2p1f_EOFA_96I_hmc.cc

@@ -183,7 +183,7 @@ int main(int argc, char **argv) {
   // 4/2 => 0.6 dH
   // 3/3 => 0.8 dH .. depth 3, slower
   //MD.MDsteps =  4;
-  MD.MDsteps =  3;
+  MD.MDsteps =  12;
   MD.trajL   = 0.5;
 
   HMCparameters HMCparams;
@@ -200,8 +200,8 @@ int main(int argc, char **argv) {
   TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
 
   CheckpointerParameters CPparams;
-  CPparams.config_prefix = "ckpoint_DDHMC_lat";
-  CPparams.rng_prefix    = "ckpoint_DDHMC_rng";
+  CPparams.config_prefix = "ckpoint_HMC_lat";
+  CPparams.rng_prefix    = "ckpoint_HMC_rng";
   CPparams.saveInterval  = 1;
   CPparams.format        = "IEEE64BIG";
   TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@@ -228,7 +228,7 @@ int main(int argc, char **argv) {
   Real pv_mass      = 1.0;
   //  std::vector<Real> hasenbusch({ 0.01, 0.045, 0.108, 0.25, 0.51 , pv_mass });
   //  std::vector<Real> hasenbusch({ light_mass, 0.01, 0.045, 0.108, 0.25, 0.51 , pv_mass });
-  std::vector<Real> hasenbusch({ light_mass, 0.005, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); // Updated
+  std::vector<Real> hasenbusch({ 0.005, 0.0145, 0.045, 0.108, 0.25, 0.51 , pv_mass }); // Updated
   //  std::vector<Real> hasenbusch({ light_mass, 0.0145, 0.045, 0.108, 0.25, 0.51 , 0.75 , pv_mass });
 
   auto GridPtr   = TheHMC.Resources.GetCartesian();
@@ -299,8 +299,8 @@ int main(int argc, char **argv) {
   ////////////////////////////////////
   // Collect actions
   ////////////////////////////////////
-  ActionLevel<HMCWrapper::Field> Level1(1);
-  ActionLevel<HMCWrapper::Field> Level2(3);
+  //  ActionLevel<HMCWrapper::Field> Level1(1);
+  ActionLevel<HMCWrapper::Field> Level2(1);
   ActionLevel<HMCWrapper::Field> Level3(15);
 
   ////////////////////////////////////
@@ -369,7 +369,7 @@ int main(int argc, char **argv) {
 	 ActionCGL, ActionCGR,
 	 DerivativeCGL, DerivativeCGR,
 	 SFRp, true);
-  //  Level2.push_back(&EOFA);
+  Level2.push_back(&EOFA);
 
   ////////////////////////////////////
   // up down action
@@ -477,7 +477,7 @@ int main(int argc, char **argv) {
   // Gauge action
   /////////////////////////////////////////////////////////////
   Level3.push_back(&GaugeAction);
-  TheHMC.TheAction.push_back(Level1);
+  //  TheHMC.TheAction.push_back(Level1);
   TheHMC.TheAction.push_back(Level2);
   TheHMC.TheAction.push_back(Level3);
   std::cout << GridLogMessage << " Action complete "<< std::endl;

TODO

@@ -1,3 +1,12 @@
+- - Slice sum optimisation & A2A - atomic addition
+- - Also faster non-atomic reduction
+- - Remaining PRs
+- - DDHMC
+  - - MixedPrec is the action eval, high precision
+  - - MixedPrecCleanup is the force eval, low precision
+
+=================
+=================
 Lattice_basis.h -- > HIP and SYCL GPU code
 
 
@@ -8,6 +17,7 @@ DDHMC
 -- Multishift Mixed Precision - DONE
 -- Pole dependent residual  - DONE
 
+
 =======
 -- comms threads issue??
 -- Part done: Staggered kernel performance on GPU

systems/Crusher/config-command

@@ -1,12 +1,13 @@
 CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 ../../configure --enable-comms=mpi-auto \
 --with-lime=$CLIME \
---enable-unified=no \
+--enable-unified=yes \
 --enable-shm=nvlink \
 --enable-tracing=timer \
 --enable-accelerator=hip \
 --enable-gen-simd-width=64 \
 --enable-simd=GPU \
+--disable-accelerator-cshift \
 --with-gmp=$OLCF_GMP_ROOT \
 --with-fftw=$FFTW_DIR/.. \
 --with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \

systems/PVC/benchmarks/run-2tile-mpi.sh

@@ -23,12 +23,7 @@ export SYCL_PI_LEVEL_ZERO_DEVICE_SCOPE_EVENTS=1
 export SYCL_PI_LEVEL_ZERO_USE_IMMEDIATE_COMMANDLISTS=1
 export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE=0
 
-for i in 0 
-do
-mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.2 --grid 32.32.32.64 --accelerator-threads $NT  --shm-mpi 1  --device-mem 32768
-mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT  --shm-mpi 1  --device-mem 32768
-done
-#mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_halo --mpi 1.1.1.2 --grid 32.32.32.64 --accelerator-threads $NT  --shm-mpi 1 > halo.2tile.1x2.log
-#mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_halo --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT  --shm-mpi 1 > halo.2tile.2x1.log
+#mpiexec -launcher ssh -n 1 -host localhost  ./wrap.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.1 --grid 32.32.32.32 --accelerator-threads $NT --comms-sequential --shm-mpi 0 > 1tile.log
 
+mpiexec -launcher ssh -n 2 -host localhost  ./wrap.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-sequential --shm-mpi 0
 

systems/PVC/config-command

@@ -14,4 +14,3 @@ INSTALL=/nfs/site/home/paboylx/prereqs/
 	LDFLAGS="-fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L$INSTALL/lib" \
 	CXXFLAGS="-fsycl-unnamed-lambda -fsycl -no-fma -I$INSTALL/include -Wno-tautological-compare"
 
-

tests/forces/Test_bdy.cc

@@ -0,0 +1,305 @@
+/*
+
+  2f Full det MdagM 10^6 force ~ 1.3e7
+rid : Message : 1767.283471 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Grid : Message : 1767.283476 s : S1 : 1.52885e+09
+Grid : Message : 1767.283480 s : S2 : 1.52886e+09
+Grid : Message : 1767.283482 s : dS : 8877.34
+Grid : Message : 1767.283483 s : dSpred : 8877.7
+Grid : Message : 1767.283484 s : diff : -0.360484
+Grid : Message : 1767.283485 s : *********************************************************
+
+  2f Full det MpcdagMpc 10^6 force ~ 1.8e6
+Grid : Message : 2399.576962 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Grid : Message : 2399.576968 s : S1 : 1.52885e+09
+Grid : Message : 2399.576972 s : S2 : 1.52886e+09
+Grid : Message : 2399.576974 s : dS : 9728.49
+Grid : Message : 2399.576975 s : dSpred : 9726.58
+Grid : Message : 2399.576976 s : diff : 1.90683
+Grid : Message : 2399.576977 s : *********************************************************
+
+  2f bdy MdagM 1500 force Force ~ 2800
+Grid : Message : 4622.385061 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Grid : Message : 4622.385067 s : S1 : 1.52885e+09
+Grid : Message : 4622.385071 s : S2 : 1.52885e+09
+Grid : Message : 4622.385072 s : dS : 25.4944
+Grid : Message : 4622.385073 s : dSpred : 25.4672
+Grid : Message : 4622.385074 s : diff : 0.0271414
+Grid : Message : 4622.385075 s : *********************************************************
+
+  2f bdy MpcdagMpc 10^6 force   Force ~ 2200
+Grid : Message : 4622.385061 s : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Grid : Message : 4622.385067 s : S1 : 1.52885e+09
+Grid : Message : 4622.385071 s : S2 : 1.52885e+09
+Grid : Message : 4622.385072 s : dS : 25.4944
+Grid : Message : 4622.385073 s : dSpred : 25.4672
+Grid : Message : 4622.385074 s : diff : 0.0271414
+Grid : Message : 4622.385075 s : *********************************************************
+  
+  1f Bdy Det
+  Optimisation log:  looser rational AND MD tolerances sloppy
+MobiusForce.221179 -- same as HMC. dS is mispredicted Forece  ~2.8
+Grid : Message : 6582.258991 s : dS : 0.024478
+Grid : Message : 6582.258992 s : dSpred : 0.00791876
+Grid : Message : 6582.258994 s : diff : 0.0165592
+
+MobiusForce.221193 -- tight rational AND MD tolerances to 1e-8 ~ 2.8 same
+Grid : Message : 1964.939209 s : S1 : 7.64404e+08
+Grid : Message : 1964.939213 s : S2 : 7.64404e+08
+Grid : Message : 1964.939215 s : dS : -0.00775838 <--- too loose even on action
+Grid : Message : 1964.939216 s : dSpred : -0.00416793 
+Grid : Message : 1964.939217 s : diff : -0.00359045
+
+MobiusForce.221394 -- looser rational, MD tol 1e-8 ~ 2.8 same
+Grid : Message : 1198.346720 s : S1 : 764404649.48886
+Grid : Message : 1198.346760 s : S2 : 764404649.5133
+Grid : Message : 1198.346780 s : dS : 0.024440884590149
+Grid : Message : 1198.346800 s : dSpred : 0.0079145154465184
+Grid : Message : 1198.346810 s : diff : 0.016526369143631
+
+MobiusForce.221394 -- tight rational, MD tol sloppy Force ~ 2.8
+Grid : Message : 2376.921950 s : S1 : 764404436.44069
+Grid : Message : 2376.921954 s : S2 : 764404436.43299
+Grid : Message : 2376.921956 s : dS : -0.0076971054077148
+Grid : Message : 2376.921958 s : dSpred : -0.0041610472282526
+Grid : Message : 2376.921959 s : diff : -0.0035360581794623
+
+*/
+
+//
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_double_ratio.cc
+
+    Copyright (C) 2022
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+typedef MobiusFermionD FermionAction;
+typedef WilsonImplD FimplD;
+typedef WilsonImplD FermionImplPolicy;
+
+template<class Gimpl>
+void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
+{
+  GridBase *UGrid = U.Grid();
+
+  std::vector<int> seeds({1,2,3,5});
+  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
+
+  LatticeColourMatrix Pmu(UGrid); 
+  LatticeGaugeField P(UGrid); 
+  LatticeGaugeField UdSdU(UGrid); 
+
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  
+  RealD eps=0.005;
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Refresh "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  
+  Gimpl::generate_momenta(P,sRNG,RNG4);
+  Filter.applyFilter(P);
+
+#if  0
+  FieldMetaData header;
+  std::string file("./ckpoint_lat.2000");
+  NerscIO::readConfiguration(U,header,file);
+#else
+  U = 1.0;
+#endif
+  action.refresh(U,sRNG,RNG4);
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+
+  RealD S1 = action.S(U);
+
+  Gimpl::update_field(P,U,eps);
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Derivative "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  action.deriv(U,UdSdU);
+  UdSdU = Ta(UdSdU);
+  Filter.applyFilter(UdSdU);
+
+  DumpSliceNorm("Force",UdSdU,Nd-1);
+  
+  Gimpl::update_field(P,U,eps);
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  
+  RealD S2 = action.S(U);
+
+  // Use the derivative
+  LatticeComplex dS(UGrid); dS = Zero();
+  for(int mu=0;mu<Nd;mu++){
+    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
+    Pmu= PeekIndex<LorentzIndex>(P,mu);
+    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
+  }
+  ComplexD dSpred    = sum(dS);
+  RealD diff =  S2-S1-dSpred.real();
+
+  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
+  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
+  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
+  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
+  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
+  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
+  //  assert(diff<1.0);
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::cout << std::setprecision(14);
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate mpi_layout  = GridDefaultMpi();
+  Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  Coordinate shm;
+  GlobalSharedMemory::GetShmDims(mpi_layout,shm);
+
+  const int Ls=12;
+  const int Nt = latt_size[3];
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+ 
+  ////////////////////////////////////////////////////////////////
+  // Domain decomposed operator
+  ////////////////////////////////////////////////////////////////
+  Coordinate CommDim(Nd);
+  for(int d=0;d<Nd;d++) CommDim[d]= (mpi_layout[d]/shm[d])>1 ? 1 : 0;
+
+  Coordinate NonDirichlet(Nd+1,0);
+  Coordinate Dirichlet(Nd+1,0);
+  Dirichlet[1] = CommDim[0]*latt_size[0]/mpi_layout[0] * shm[0];
+  Dirichlet[2] = CommDim[1]*latt_size[1]/mpi_layout[1] * shm[1];
+  Dirichlet[3] = CommDim[2]*latt_size[2]/mpi_layout[2] * shm[2];
+  Dirichlet[4] = CommDim[3]*latt_size[3]/mpi_layout[3] * shm[3];
+
+  Coordinate Block4(Nd);
+  Block4[0] = Dirichlet[1];
+  Block4[1] = Dirichlet[2];
+  Block4[2] = Dirichlet[3];
+  Block4[3] = Dirichlet[4];
+
+  std::vector<Complex> boundary = {1,1,1,-1};
+  FermionAction::ImplParams Params(boundary);
+  FermionAction::ImplParams ParamsDir(boundary);
+  Params.dirichlet=NonDirichlet;
+  ParamsDir.dirichlet=Dirichlet;
+  ParamsDir.partialDirichlet=1;
+
+  ///////////////////// Gauge Field and Gauge Forces ////////////////////////////
+  LatticeGaugeField U(UGrid);
+
+  RealD beta=6.0;
+  WilsonGaugeActionR PlaqAction(beta);
+  IwasakiGaugeActionR RectAction(beta);
+
+  MomentumFilterNone<LatticeGaugeField> FilterNone;
+  ForceTest<GimplTypesR>(PlaqAction,U,FilterNone);
+  ForceTest<GimplTypesR>(RectAction,U,FilterNone);
+
+  ////////////////////////////////////
+  // Action
+  ////////////////////////////////////
+  RealD mass=0.00078; 
+  RealD pvmass=1.0; 
+  RealD M5=1.8; 
+  RealD b=1.5;
+  RealD c=0.5;
+  
+  // Double versions
+  FermionAction DdwfPeriodic(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,Params);
+  FermionAction PVPeriodic  (U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pvmass,M5,b,c,Params);
+  FermionAction DdwfDirichlet(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,ParamsDir);
+
+  double StoppingCondition = 1.0e-8;
+  double MaxCGIterations = 50000;
+  ConjugateGradient<LatticeFermion>  CG(StoppingCondition,MaxCGIterations);
+  
+  //////////////////// Two Flavour Determinant Ratio ///////////////////////////////
+  TwoFlavourRatioPseudoFermionAction<FimplD> Nf2(PVPeriodic, DdwfPeriodic,CG,CG);
+  //  ForceTest<GimplTypesR>(Nf2,U,FilterNone);
+
+  //////////////////// Two Flavour Determinant force test Even Odd ///////////////////////////////
+  TwoFlavourEvenOddRatioPseudoFermionAction<FimplD> Nf2eo(PVPeriodic, DdwfPeriodic,CG,CG);
+  //  ForceTest<GimplTypesR>(Nf2eo,U,FilterNone);
+
+  //////////////////// Domain forces ////////////////////
+  int Width=4;
+  DDHMCFilter<WilsonImplD::Field> DDHMCFilter(Block4,Width);
+  
+  //////////////////// Two flavour boundary det  ////////////////////
+  TwoFlavourRatioPseudoFermionAction<FimplD> BdyNf2(DdwfDirichlet, DdwfPeriodic,CG,CG);
+  //  ForceTest<GimplTypesR>(BdyNf2,U,DDHMCFilter);
+
+  //////////////////// Two flavour eo boundary det  ////////////////////
+  TwoFlavourEvenOddRatioPseudoFermionAction<FimplD> BdyNf2eo(DdwfDirichlet, DdwfPeriodic,CG,CG);
+  //  ForceTest<GimplTypesR>(BdyNf2eo,U,DDHMCFilter);
+
+  //////////////////// One flavour boundary det  ////////////////////
+  OneFlavourRationalParams OFRp; // Up/down
+  OFRp.lo       = 4.0e-5;
+  OFRp.hi       = 90.0;
+  OFRp.MaxIter  = 60000;
+  OFRp.tolerance= 1.0e-8;
+  OFRp.mdtolerance= 1.0e-6;
+  OFRp.degree   = 18;
+  OFRp.precision= 80;
+  OFRp.BoundsCheckFreq=0;
+  std::vector<RealD> ActionTolByPole({
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8
+    });
+  std::vector<RealD> MDTolByPole({
+      1.0e-6,3.0e-7,1.0e-7,1.0e-7,  // Orig sloppy
+      //      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8
+    });
+  OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> BdySqrt(DdwfDirichlet,DdwfPeriodic,OFRp);
+  ForceTest<GimplTypesR>(BdySqrt,U,DDHMCFilter);
+
+  Grid_finalize();
+}

tests/forces/Test_double_ratio.cc

@@ -476,6 +476,20 @@ int main (int argc, char ** argv)
   //  ForceTest<GimplTypesR>(BdyNf2eo,U,DDHMCFilter);
 
   //////////////////// One flavour boundary det  ////////////////////
+  RationalActionParams OFRp; // Up/down
+  OFRp.lo       = 6.0e-5;
+  OFRp.hi       = 90.0;
+  OFRp.inv_pow  = 2;
+  OFRp.MaxIter  = SP_iters; // get most shifts by 2000, stop sharing space
+  OFRp.action_tolerance= 1.0e-8;
+  OFRp.action_degree   = 18;
+  OFRp.md_tolerance= 1.0e-5;
+  OFRp.md_degree   = 14;
+  //  OFRp.degree   = 20; converges
+  //  OFRp.degree   = 16;
+  OFRp.precision= 80;
+  OFRp.BoundsCheckFreq=0;
+  /*
   OneFlavourRationalParams OFRp; // Up/down
   OFRp.lo       = 4.0e-5;
   OFRp.hi       = 90.0;
@@ -485,6 +499,23 @@ int main (int argc, char ** argv)
   OFRp.degree   = 18;
   OFRp.precision= 80;
   OFRp.BoundsCheckFreq=0;
+  */
+  std::vector<RealD> ActionTolByPole({
+      1.0e-7,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8
+    });
+  std::vector<RealD> MDTolByPole({
+      1.6e-5,5.0e-6,1.0e-6,3.0e-7, // soften convergence more more
+      //      1.0e-6,3.0e-7,1.0e-7,1.0e-7,
+      //      3.0e-6,1.0e-6,1.0e-7,1.0e-7, // soften convergence
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
+      1.0e-8,1.0e-8
+    });
+  /*
   std::vector<RealD> ActionTolByPole({
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
@@ -499,9 +530,9 @@ int main (int argc, char ** argv)
       //      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8,1.0e-8,1.0e-8,
-      1.0e-8,1.0e-8,1.0e-8,1.0e-8,
       1.0e-8,1.0e-8
     });
+  */
   OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> BdySqrt(DdwfDirichlet,DdwfPeriodic,OFRp);
   BdySqrt.SetTolerances(ActionTolByPole,MDTolByPole);
   ForceTest<GimplTypesR>(BdySqrt,U,DDHMCFilter);