Merge branch 'release/0.10.0'

Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h

@@ -38,15 +38,19 @@ NAMESPACE_BEGIN(Grid);
     // cf. GeneralEvenOddRational.h for details
     /////////////////////////////////////////////////////////////////////////////////////////////////////////////
       
-    template<class ImplD, class ImplF>
+    template<class ImplD, class ImplF, class ImplD2>
     class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> {
     private:
+      typedef typename ImplD2::FermionField FermionFieldD2;
       typedef typename ImplD::FermionField FermionFieldD;
       typedef typename ImplF::FermionField FermionFieldF;
 
       FermionOperator<ImplD> & NumOpD;
       FermionOperator<ImplD> & DenOpD;
 
+      FermionOperator<ImplD2> & NumOpD2;
+      FermionOperator<ImplD2> & DenOpD2;
+     
       FermionOperator<ImplF> & NumOpF;
       FermionOperator<ImplF> & DenOpF;
 
@@ -60,31 +64,40 @@ NAMESPACE_BEGIN(Grid);
 	ConjugateGradientMultiShift<FermionFieldD> msCG(MaxIter, approx);
 	msCG(schurOp,in, out);
 #else
-	SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
+	SchurDifferentiableOperator<ImplD2> schurOpD2(numerator ? NumOpD2 : DenOpD2);
 	SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
-	FermionFieldD inD(NumOpD.FermionRedBlackGrid());
+	FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
-	FermionFieldD outD(NumOpD.FermionRedBlackGrid());
+	FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
 
 	// Action better with higher precision?
-	ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
+	ConjugateGradientMultiShiftMixedPrec<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
-	msCG(schurOpD, in, out);
+	precisionChange(inD2,in);
+	std::cout << "msCG single solve "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
+	msCG(schurOpD2, inD2, outD2);
+	precisionChange(out,outD2);
 #endif
       }
       virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &out){
-	SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
+	SchurDifferentiableOperator<ImplD2> schurOpD2(numerator ? NumOpD2 : DenOpD2);
 	SchurDifferentiableOperator<ImplF>  schurOpF (numerator ? NumOpF  : DenOpF);
 
-	FermionFieldD inD(NumOpD.FermionRedBlackGrid());
+	FermionFieldD2 inD2(NumOpD2.FermionRedBlackGrid());
-	FermionFieldD outD(NumOpD.FermionRedBlackGrid());
+	FermionFieldD2 outD2(NumOpD2.FermionRedBlackGrid());
-	std::vector<FermionFieldD> out_elemsD(out_elems.size(),NumOpD.FermionRedBlackGrid());
+	std::vector<FermionFieldD2> out_elemsD2(out_elems.size(),NumOpD2.FermionRedBlackGrid());
-	ConjugateGradientMultiShiftMixedPrecCleanup<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
+	ConjugateGradientMultiShiftMixedPrecCleanup<FermionFieldD2, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
-	msCG(schurOpD, in, out_elems, out);
+	precisionChange(inD2,in);
+	std::cout << "msCG in "<<norm2(inD2)<<" " <<norm2(in)<<std::endl;
+	msCG(schurOpD2, inD2, out_elemsD2, outD2);
+	precisionChange(out,outD2);
+	for(int i=0;i<out_elems.size();i++){
+	  precisionChange(out_elems[i],out_elemsD2[i]);
+	}
       }
       //Allow derived classes to override the gauge import
       virtual void ImportGauge(const typename ImplD::GaugeField &Ud){
 
 	typename ImplF::GaugeField Uf(NumOpF.GaugeGrid());
-	typename ImplD::GaugeField Ud2(NumOpD.GaugeGrid());
+	typename ImplD2::GaugeField Ud2(NumOpD2.GaugeGrid());
 	precisionChange(Uf, Ud);
 	precisionChange(Ud2, Ud);
 
@@ -96,18 +109,20 @@ NAMESPACE_BEGIN(Grid);
 	NumOpF.ImportGauge(Uf);
 	DenOpF.ImportGauge(Uf);
 
-	NumOpD.ImportGauge(Ud2);
+	NumOpD2.ImportGauge(Ud2);
-	DenOpD.ImportGauge(Ud2);
+	DenOpD2.ImportGauge(Ud2);
       }
       
     public:
       GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD>  &_NumOpD, FermionOperator<ImplD>  &_DenOpD, 
 							      FermionOperator<ImplF>  &_NumOpF, FermionOperator<ImplF>  &_DenOpF, 
+							      FermionOperator<ImplD2>  &_NumOpD2, FermionOperator<ImplD2>  &_DenOpD2, 
 							      const RationalActionParams & p, Integer _ReliableUpdateFreq
 							      ) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p),
 								  ReliableUpdateFreq(_ReliableUpdateFreq),
 								  NumOpD(_NumOpD), DenOpD(_DenOpD),
-								  NumOpF(_NumOpF), DenOpF(_DenOpF)
+								  NumOpF(_NumOpF), DenOpF(_DenOpF),
+								  NumOpD2(_NumOpD2), DenOpD2(_DenOpD2)
       {}
       
       virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";}

Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h

@@ -67,9 +67,9 @@ NAMESPACE_BEGIN(Grid);
       virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";}      
     };
 
-    template<class Impl,class ImplF>
+    template<class Impl,class ImplF,class ImplD2>
     class OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction
-      : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF> {
+      : public GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF,ImplD2> {
     public:
       typedef OneFlavourRationalParams Params;
     private:
@@ -91,9 +91,11 @@ NAMESPACE_BEGIN(Grid);
 								 FermionOperator<Impl>  &_DenOp, 
 								 FermionOperator<ImplF>  &_NumOpF, 
 								 FermionOperator<ImplF>  &_DenOpF, 
+								 FermionOperator<ImplD2>  &_NumOpD2, 
+								 FermionOperator<ImplD2>  &_DenOpD2, 
 								 const Params & p, Integer ReliableUpdateFreq
 							) : 
-	GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF>(_NumOp, _DenOp,_NumOpF, _DenOpF, transcribe(p),ReliableUpdateFreq){}
+	GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<Impl,ImplF,ImplD2>(_NumOp, _DenOp,_NumOpF, _DenOpF,_NumOpD2, _DenOpD2, transcribe(p),ReliableUpdateFreq){}
 
       virtual std::string action_name(){return "OneFlavourEvenOddRatioRationalPseudoFermionAction";}      
     };

HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc

@@ -451,7 +451,7 @@ int main(int argc, char **argv) {
   
 #define MIXED_PRECISION
 #ifdef MIXED_PRECISION
-  std::vector<GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF> *> Bdys;
+  std::vector<GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy> *> Bdys;
 #else
   std::vector<GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> *> Bdys;
 #endif
@@ -526,13 +526,15 @@ 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
-      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF>(
+      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy>(
 			   *Numerators[h],*Denominators[h],
 			   *NumeratorsF[h],*DenominatorsF[h],
+			   *Numerators[h],*Denominators[h],
 			   OFRp, SP_iters) );
-      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF>(
+      Bdys.push_back( new GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicy>(
 			   *Numerators[h],*Denominators[h],
 			   *NumeratorsF[h],*DenominatorsF[h],
+			   *Numerators[h],*Denominators[h],
 			   OFRp, SP_iters) );
 #else
       Bdys.push_back( new GeneralEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],OFRp));

HMC/Mobius2p1f_EOFA_96I_hmc.cc

@@ -164,6 +164,11 @@ 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;
 
   //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -267,6 +272,7 @@ 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
@@ -388,13 +394,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> *> Bdys;
+  std::vector<OneFlavourEvenOddRatioRationalMixedPrecPseudoFermionAction<FermionImplPolicy,FermionImplPolicyF,FermionImplPolicyD2> *> Bdys;
 #else
   std::vector<OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> *> Bdys;
 #endif

systems/Crusher/config-command

@@ -1,7 +1,7 @@
 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 \

systems/Crusher/sourceme.sh

@@ -5,8 +5,8 @@ module load emacs
 #module load gperftools
 module load PrgEnv-gnu
 module load rocm/5.3.0
-#module load cray-mpich/8.1.16
+module load cray-mpich/8.1.16
-module load cray-mpich/8.1.17
+#module load cray-mpich/8.1.17
 module load gmp
 module load cray-fftw
 module load craype-accel-amd-gfx90a

tests/Test_cayley_even_odd_vec.cc

@@ -73,12 +73,12 @@ int main (int argc, char ** argv)
   RealD M5  =1.8;
 
   std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
-  std::cout<<GridLogMessage <<"DomainWallFermion vectorised test"<<std::endl;
+  std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
   std::cout<<GridLogMessage<<"**************************************************************"<<std::endl;
   std::vector<Complex> boundary = {1,1,1,-1};
   DomainWallFermionD::ImplParams Params(boundary);
-  Coordinate Dirichlet({0,8,8,16,32});
+  //  Coordinate Dirichlet({0,8,8,16,32});
-  Params.dirichlet=Dirichlet;
+  //  Params.dirichlet=Dirichlet;
 
   DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,Params);
   TestWhat<DomainWallFermionD>(Ddwf,FGrid,FrbGrid,UGrid,mass,M5,&RNG4,&RNG5);