Merge 32e6d58356 into 557fa483ff

Grid/cartesian/Cartesian_base.h

@@ -91,7 +91,6 @@ public:
   ////////////////////////////////////////////////////////////////
   virtual int CheckerBoarded(int dim)=0;
   virtual int CheckerBoard(const Coordinate &site)=0;
-  virtual int CheckerDim(void){ return 0; };
   virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;
   virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0;
   virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0;

Grid/cartesian/Cartesian_red_black.h

@@ -60,7 +60,6 @@ public:
   int              _checker_dim;
   std::vector<int> _checker_board;
 
-  virtual int CheckerDim(void){ return _checker_dim; };
   virtual int CheckerBoarded(int dim){
     if( dim==_checker_dim) return 1;
     else return 0;

Grid/lattice/Lattice_transfer.h

@@ -42,21 +42,50 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
     assert((fine->_rdimensions[d] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[d]); 
   }
 }
+
  
 ////////////////////////////////////////////////////////////////////////////////////////////
 // remove and insert a half checkerboard
 ////////////////////////////////////////////////////////////////////////////////////////////
-
 template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full)
 {
-  acceleratorPickCheckerboard(cb,half,full);
+  half.Checkerboard() = cb;
+
+  autoView( half_v, half, CpuWrite);
+  autoView( full_v, full, CpuRead);
+  thread_for(ss, full.Grid()->oSites(),{
+    int cbos;
+    Coordinate coor;
+    full.Grid()->oCoorFromOindex(coor,ss);
+    cbos=half.Grid()->CheckerBoard(coor);
+
+    if (cbos==cb) {
+      int ssh=half.Grid()->oIndex(coor);
+      half_v[ssh] = full_v[ss];
+    }
+  });
 }
 template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half)
 {
-  acceleratorSetCheckerboard(full,half);
+  int cb = half.Checkerboard();
+  autoView( half_v , half, CpuRead);
+  autoView( full_v , full, CpuWrite);
+  thread_for(ss,full.Grid()->oSites(),{
+
+    Coordinate coor;
+    int cbos;
+
+    full.Grid()->oCoorFromOindex(coor,ss);
+    cbos=half.Grid()->CheckerBoard(coor);
+      
+    if (cbos==cb) {
+      int ssh=half.Grid()->oIndex(coor);
+      full_v[ss]=half_v[ssh];
+    }
+  });
 }
 
-template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int dummy=0)
+template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int checker_dim_half=0)
 {
   half.Checkerboard() = cb;
   autoView(half_v, half, AcceleratorWrite);
@@ -66,7 +95,6 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
-  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss, full.Grid()->oSites(),full.Grid()->Nsimd(),{
     
     Coordinate coor;
@@ -91,7 +119,7 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
     }
   });
 }
-template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int dummy=0)
+template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int checker_dim_half=0)
 {
   int cb = half.Checkerboard();
   autoView(half_v , half, AcceleratorRead);
@@ -101,7 +129,6 @@ template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
-  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss,full.Grid()->oSites(),full.Grid()->Nsimd(),{
 
     Coordinate coor;

Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -86,8 +86,13 @@ public:
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);
 
+#if defined(GRID_CUDA) || defined(GRID_HIP)  || defined(GRID_SYCL)
+    acceleratorSetCheckerboard(Force,ForceE);
+    acceleratorSetCheckerboard(Force,ForceO);
+#else
     setCheckerboard(Force,ForceE); 
     setCheckerboard(Force,ForceO);
+#endif
     Force=-Force;
 
     delete forcecb;
@@ -130,8 +135,13 @@ public:
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);
 
+#if defined(GRID_CUDA) || defined(GRID_HIP)  || defined(GRID_SYCL)
+    acceleratorSetCheckerboard(Force,ForceE);
+    acceleratorSetCheckerboard(Force,ForceO);
+#else
     setCheckerboard(Force,ForceE); 
     setCheckerboard(Force,ForceO);
+#endif
     Force=-Force;
 
     delete forcecb;

Grid/qcd/smearing/GaugeConfigurationMasked.h

@@ -32,9 +32,7 @@ private:
   //  Smear_Stout<Gimpl> *StoutSmearing;
   //  std::vector<GaugeField> SmearedSet;
   
-  GridRedBlackCartesian * UrbGrid; // keep a copy of the redblack grid for life of object
   std::vector<LatticeLorentzComplex> masks;
-  std::vector<int> cbs;
 
   typedef typename SU3Adjoint::AMatrix AdjMatrix;
   typedef typename SU3Adjoint::LatticeAdjMatrix  AdjMatrixField;
@@ -149,25 +147,6 @@ private:
     }
     pokeLorentz(Fdet, Fdet_pol, nu);
   }
-
-  void Compute_MpInvJx_dNxxdSy(int cb,
-			       const GaugeLinkField &PlaqL,
-			       const GaugeLinkField &PlaqR,
-			       AdjMatrixField MpInvJx,
-			       AdjVectorField &Fdet2 )
-  {
-    GaugeLinkField PlaqLeo(UrbGrid);
-    GaugeLinkField PlaqReo(UrbGrid);
-    AdjMatrixField MpInvJxeo(UrbGrid);
-    AdjVectorField Fdet2eo(UrbGrid);
-    pickCheckerboard(cb,PlaqLeo,PlaqL);
-    pickCheckerboard(cb,PlaqReo,PlaqR);
-    pickCheckerboard(cb,MpInvJxeo,MpInvJx);
-    Fdet2eo.Checkerboard()=cb;
-    Compute_MpInvJx_dNxxdSy(PlaqLeo,PlaqReo,MpInvJxeo,Fdet2eo);
-    setCheckerboard(Fdet2,Fdet2eo);
-  }
-  
   void Compute_MpInvJx_dNxxdSy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
   {
     GaugeLinkField UtaU(PlaqL.Grid());
@@ -299,9 +278,8 @@ public:
     ////////////////////////////////////////////////////////////////////////////////
     // Mask the gauge field
     ////////////////////////////////////////////////////////////////////////////////
-    int cb = cbs[smr];
     auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask
-    
+
     Umsk = U;
     ApplyMask(Umsk,smr);
     Utmp = peekLorentz(Umsk,mu);
@@ -464,7 +442,7 @@ public:
     AdjMatrixField MpInvJx_nu(grid);
     MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor
 
-    Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx,FdetV);
+    Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
     Fdet2_mu=FdetV;
     Fdet1_mu=Zero();
     
@@ -521,7 +499,7 @@ public:
 
 	time=-usecond();
 	PlaqR=(-1.0)*PlaqR;
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
 	Fdet2_nu = FdetV;
 	time+=usecond();
 	std::cout << GridLogMessage << "Compute_MpInvJx_dNxxSy (occurs 6x) took "<<time<< " us"<<std::endl;
@@ -542,7 +520,7 @@ public:
 	
 
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 	
 	///////////////// -ve nu /////////////////
@@ -561,7 +539,7 @@ public:
 	Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 	
 	// x==
@@ -582,7 +560,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
 	MpInvJx_nu = Cshift(MpInvJx_nu,nu,1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 
 	/////////////////////////////////////////////////////////////////////
@@ -611,7 +589,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,-1);
 
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
 
 	//  __
@@ -631,7 +609,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
 
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
 	
       }
@@ -953,10 +931,6 @@ private:
 public:
 
   /* Standard constructor */
-  virtual ~SmearedConfigurationMasked()
-  {
-    delete UrbGrid;
-  }
   SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout)
     : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
   {
@@ -965,6 +939,7 @@ public:
     // was resized in base class
     assert(this->SmearedSet.size()==Nsmear);
     
+    GridRedBlackCartesian * UrbGrid;
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);
     LatticeComplex one(_UGrid); one = ComplexD(1.0,0.0);
     LatticeComplex tmp(_UGrid);
@@ -972,11 +947,10 @@ public:
     for (unsigned int i = 0; i < this->smearingLevels; ++i) {
 
       masks.push_back(*(new LatticeLorentzComplex(_UGrid)));
+
       int mu= (i/2) %Nd;
       int cb= (i%2);
       LatticeComplex tmpcb(UrbGrid);
-
-      cbs.push_back(cb);
 	
       masks[i]=Zero();
       ////////////////////
@@ -988,6 +962,7 @@ public:
       PokeIndex<LorentzIndex>(masks[i],tmp, mu);
 	
     }
+    delete UrbGrid;
   }
   
   virtual void smeared_force(GaugeField &SigmaTilde) 

Grid/qcd/utils/GaugeGroup.h

@@ -418,32 +418,32 @@ static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in,
   int hNNm1= NNm1/2;
   RealD sqrt_2 = sqrt(2.0);
   Complex ci(0.0,1.0);
-
+  for(int su2Index=0;su2Index<hNNm1;su2Index++){
-  const int nsimd=  Matrix::Nsimd();
+    int i1, i2;
-  accelerator_for(ss,grid->oSites(),nsimd,{
+    su2SubGroupIndex(i1, i2, su2Index);
-      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+    int ax = su2Index*2;
-	int i1, i2;
+    int ay = su2Index*2+1;
-	su2SubGroupIndex(i1, i2, su2Index);
+    accelerator_for(ss,grid->oSites(),1,{
-	int ax = su2Index*2;
-	int ay = su2Index*2+1;
 	// in is traceless ANTI-hermitian whereas Grid generators are Hermitian.
 	// trace( Ta x Ci in)
 	// Bet I need to move to real part with mult by -i
-	coalescedWrite(out_v[ss]()()(ax,b),0.5*(real(in_v(ss)()()(i2,i1)) - real(in_v(ss)()()(i1,i2))));
+	out_v[ss]()()(ax,b) = 0.5*(real(in_v[ss]()()(i2,i1)) - real(in_v[ss]()()(i1,i2)));
-	coalescedWrite(out_v[ss]()()(ay,b),0.5*(imag(in_v(ss)()()(i1,i2)) + imag(in_v(ss)()()(i2,i1))));
+	out_v[ss]()()(ay,b) = 0.5*(imag(in_v[ss]()()(i1,i2)) + imag(in_v[ss]()()(i2,i1)));
-      }
+      });
-      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+  }
-	int k = diagIndex + 1; // diagIndex starts from 0
+  for(int diagIndex=0;diagIndex<N-1;diagIndex++){
-	int a = NNm1+diagIndex;
+    int k = diagIndex + 1; // diagIndex starts from 0
-	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+    int a = NNm1+diagIndex;
-	auto tmp = in_v(ss)()()(0,0);
+    RealD scale = 1.0/sqrt(2.0*k*(k+1));
+    accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
+	auto tmp = in_v[ss]()()(0,0);
 	for(int i=1;i<k;i++){
-	  tmp=tmp+in_v(ss)()()(i,i);
+	  tmp=tmp+in_v[ss]()()(i,i);
 	}
-	tmp = tmp - in_v(ss)()()(k,k)*k;
+	tmp = tmp - in_v[ss]()()(k,k)*k;
-	coalescedWrite(out_v[ss]()()(a,b),imag(tmp) * scale);
+	out_v[ss]()()(a,b) =imag(tmp) * scale;
-      }
+      });
-    });
+    }
 }
 
   

Grid/qcd/utils/SUn.impl.h

@@ -118,7 +118,7 @@ static void generatorDiagonal(int diagIndex, iGroupMatrix<cplx> &ta) {
 ////////////////////////////////////////////////////////////////////////
 // Map a su2 subgroup number to the pair of rows that are non zero
 ////////////////////////////////////////////////////////////////////////
-static accelerator_inline void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
+static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
   assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2));
 
   int spare = su2_index;