Added a test an fixed it for red black precon Ls innermost vectorised DWF

benchmarks/Benchmark_dwf.cc

@@ -71,7 +71,6 @@ int main (int argc, char ** argv)
 
   std::cout << GridLogMessage << "Making s innermost grids"<<std::endl;
   GridCartesian         * sUGrid   = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
-  GridRedBlackCartesian * sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
   GridCartesian         * sFGrid   = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
   std::cout << GridLogMessage << "Making s innermost rb grids"<<std::endl;
   GridRedBlackCartesian * sFrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
@@ -162,13 +161,13 @@ int main (int argc, char ** argv)
     //    Dw.Report();
   }
 
-  if (0)
+  if (1)
   {
     typedef WilsonFermion5D<DomainWallRedBlack5dImplF> WilsonFermion5DF;
     LatticeFermionF ssrc(sFGrid);
     LatticeFermionF sref(sFGrid);
     LatticeFermionF sresult(sFGrid);
-    WilsonFermion5DF sDw(1,Umu,*sFGrid,*sFrbGrid,*sUGrid,*sUrbGrid,M5,params);
+    WilsonFermion5DF sDw(1,Umu,*sFGrid,*sFrbGrid,*sUGrid,M5,params);
   
     for(int x=0;x<latt4[0];x++){
     for(int y=0;y<latt4[1];y++){
@@ -196,7 +195,7 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << "mflop/s per node =  "<< flops/(t1-t0)/NP<<std::endl;
     //  sDw.Report();
   
-    if(1){
+    if(0){
       for(int i=0;i< PerformanceCounter::NumTypes(); i++ ){
 	sDw.Dhop(ssrc,sresult,0);
 	PerformanceCounter Counter(i);
@@ -226,6 +225,53 @@ int main (int argc, char ** argv)
     }}}}}
     std::cout<<" difference between normal and simd is "<<sum<<std::endl;
 
+
+    if (1) {
+
+      LatticeFermionF sr_eo(sFGrid);
+      LatticeFermionF serr(sFGrid);
+
+      LatticeFermion ssrc_e (sFrbGrid);
+      LatticeFermion ssrc_o (sFrbGrid);
+      LatticeFermion sr_e   (sFrbGrid);
+      LatticeFermion sr_o   (sFrbGrid);
+
+      pickCheckerboard(Even,ssrc_e,ssrc);
+      pickCheckerboard(Odd,ssrc_o,ssrc);
+
+      setCheckerboard(sr_eo,ssrc_o);
+      setCheckerboard(sr_eo,ssrc_e);
+      serr = sr_eo-ssrc; 
+      std::cout<<GridLogMessage << "EO src norm diff   "<< norm2(serr)<<std::endl;
+
+      sr_e = zero;
+      sr_o = zero;
+
+      double t0=usecond();
+      for(int i=0;i<ncall;i++){
+	sDw.DhopEO(ssrc_o,sr_e,DaggerNo);
+      }
+      double t1=usecond();
+
+      double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
+      double flops=(1344.0*volume*ncall)/2;
+
+      std::cout<<GridLogMessage << "sDeo mflop/s =   "<< flops/(t1-t0)<<std::endl;
+      std::cout<<GridLogMessage << "sDeo mflop/s per node   "<< flops/(t1-t0)/NP<<std::endl;
+
+      sDw.DhopEO(ssrc_o,sr_e,DaggerNo);
+      sDw.DhopOE(ssrc_e,sr_o,DaggerNo);
+      sDw.Dhop  (ssrc  ,sresult,DaggerNo);
+
+      pickCheckerboard(Even,ssrc_e,sresult);
+      pickCheckerboard(Odd ,ssrc_o,sresult);
+      ssrc_e = ssrc_e - sr_e;
+      std::cout<<GridLogMessage << "sE norm diff   "<< norm2(ssrc_e)<<std::endl;
+      ssrc_o = ssrc_o - sr_o;
+      std::cout<<GridLogMessage << "sO norm diff   "<< norm2(ssrc_o)<<std::endl;
+    }
+
+
   }
 
   if (1)

benchmarks/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Benchmark_comms Benchmark_dwf Benchmark_memory_asynch Benchmark_memory_bandwidth Benchmark_su3 Benchmark_wilson Benchmark_zmm
+bin_PROGRAMS = Benchmark_comms Benchmark_dwf Benchmark_dwf_ntpf Benchmark_memory_asynch Benchmark_memory_bandwidth Benchmark_su3 Benchmark_wilson Benchmark_zmm
 
 
 Benchmark_comms_SOURCES=Benchmark_comms.cc
@@ -10,6 +10,10 @@ Benchmark_dwf_SOURCES=Benchmark_dwf.cc
 Benchmark_dwf_LDADD=-lGrid
 
 
+Benchmark_dwf_ntpf_SOURCES=Benchmark_dwf_ntpf.cc
+Benchmark_dwf_ntpf_LDADD=-lGrid
+
+
 Benchmark_memory_asynch_SOURCES=Benchmark_memory_asynch.cc
 Benchmark_memory_asynch_LDADD=-lGrid
 

lib/cartesian/Cartesian_red_black.h

@@ -170,17 +170,13 @@ public:
 	// Use a reduced simd grid
 	_simd_layout[d] = simd_layout[d];
 	_rdimensions[d]= _ldimensions[d]/_simd_layout[d];
+	assert(_rdimensions[d]>0);
 
 	// all elements of a simd vector must have same checkerboard.
-	//
+	// If Ls vectorised, this must still be the case; e.g. dwf rb5d
-	// If Ls vectorised, this must still be the case; e.g. 
-	// layout == 8 , require _rdimensions[d] >= 2;
-	//	std::cout << GridLogMessage << " dimension "<<d<< " layout "<< _simd_layout[d] <<" ldimension "<<_ldimensions[d]<<std::endl;
 	if ( _simd_layout[d]>1 ) {
-	  if ( _rdimensions[d]&0x1 ) { 
+	  if ( d != _checker_dim ) { 
-	    std::cout << GridLogMessage << " dimension "<<d<< " bad layout "<<
+	    assert( (_rdimensions[d]&0x1) == 0 );
-	      _simd_layout[d] << " " << _ldimensions[d] << std::endl;
-	    //	    assert(0);
 	  }
 	}
 

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -103,26 +103,24 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
 }  
 
 template<class Impl>
-WilsonFermion5D<Impl>::WilsonFermion5D(int simd, GaugeField &_Umu,
+WilsonFermion5D<Impl>::WilsonFermion5D(int simd,GaugeField &_Umu,
 				       GridCartesian         &FiveDimGrid,
 				       GridRedBlackCartesian &FiveDimRedBlackGrid,
 				       GridCartesian         &FourDimGrid,
-				       GridRedBlackCartesian &FourDimRedBlackGrid,
 				       RealD _M5,const ImplParams &p) :
   Kernels(p),
   _FiveDimGrid        (&FiveDimGrid),
   _FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
   _FourDimGrid        (&FourDimGrid),
-  _FourDimRedBlackGrid(&FourDimRedBlackGrid),
   Stencil    (_FiveDimGrid,npoint,Even,directions,displacements),
   StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
   StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
   M5(_M5),
   Umu(_FourDimGrid),
-  UmuEven(_FourDimRedBlackGrid),
+  UmuEven(_FourDimGrid),
-  UmuOdd (_FourDimRedBlackGrid),
+  UmuOdd (_FourDimGrid),
   Lebesgue(_FourDimGrid),
-  LebesgueEvenOdd(_FourDimRedBlackGrid)
+  LebesgueEvenOdd(_FourDimGrid)
 {
   int nsimd = Simd::Nsimd();
 
@@ -131,7 +129,6 @@ WilsonFermion5D<Impl>::WilsonFermion5D(int simd, GaugeField &_Umu,
   assert(FiveDimRedBlackGrid._ndimension==5);
   assert(FiveDimRedBlackGrid._checker_dim==0); // Checkerboard the s-direction
   assert(FourDimGrid._ndimension==4);
-  assert(FourDimRedBlackGrid._ndimension==4);
 
   // Dimension zero of the five-d is the Ls direction
   Ls=FiveDimGrid._fdimensions[0];
@@ -144,15 +141,10 @@ WilsonFermion5D<Impl>::WilsonFermion5D(int simd, GaugeField &_Umu,
 
   // Other dimensions must match the decomposition of the four-D fields 
   for(int d=0;d<4;d++){
-    assert(FourDimRedBlackGrid._fdimensions[d]  ==FourDimGrid._fdimensions[d]);
     assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]);
-
-    assert(FourDimRedBlackGrid._processors[d]   ==FourDimGrid._processors[d]);
     assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]);
 
     assert(FourDimGrid._simd_layout[d]=1);
-    assert(FourDimRedBlackGrid._simd_layout[d]  ==1);
-    assert(FourDimRedBlackGrid._simd_layout[d]  ==1);
     assert(FiveDimRedBlackGrid._simd_layout[d+1]==1);
 
     assert(FiveDimGrid._fdimensions[d+1]        ==FourDimGrid._fdimensions[d]);
@@ -160,8 +152,13 @@ WilsonFermion5D<Impl>::WilsonFermion5D(int simd, GaugeField &_Umu,
     assert(FiveDimGrid._simd_layout[d+1]        ==FourDimGrid._simd_layout[d]);
   }
 
-  // Allocate the required comms buffer
+  {
-  ImportGauge(_Umu);
+    GaugeField HUmu(_Umu._grid);
+    HUmu = _Umu*(-0.5);
+    Impl::DoubleStore(GaugeGrid(),Umu,HUmu);
+    UmuEven=Umu;// Really want a reference.
+    UmuOdd =Umu;
+  }
 }  
 
 

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -134,7 +134,6 @@ namespace Grid {
 		      GridCartesian         &FiveDimGrid,
 		      GridRedBlackCartesian &FiveDimRedBlackGrid,
 		      GridCartesian         &FourDimGrid,
-		      GridRedBlackCartesian &FourDimRedBlackGrid,
 		      double _M5,const ImplParams &p= ImplParams());
 
       // DoubleStore

tests/Test_dwf_rb5d.cc

@@ -97,7 +97,7 @@ int main (int argc, char ** argv)
   RealD M5  =1.8;
   typename WilsonFermion5DR::ImplParams params;
 
-  WilsonFermion5DR Dw(1,Umu,*FGrid,*FrbGrid,*sUGrid,*sUrbGrid,M5,params);
+  WilsonFermion5DR Dw(1,Umu,*FGrid,*FrbGrid,*sUGrid,M5,params);
 
   Dw.Dhop(src,result,0);