Some unary ops and coarse grid support

lib/Grid.h

@@ -55,7 +55,6 @@
 #include <Cartesian.h> // subdir aggregate
 #include <Tensors.h>   // subdir aggregate
 #include <Lattice.h>   // subdir aggregate
-#include <Comparison.h>
 #include <Cshift.h>    // subdir aggregate
 #include <Stencil.h>   // subdir aggregate
 #include <Algorithms.h>// subdir aggregate

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h
 
 CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/algorithms/CoarsenedMatrix.h

@@ -6,12 +6,42 @@
 namespace Grid {
 
   class Geometry {
+    //    int dimension;
   public:
     int npoint;
-    int dimension;
     std::vector<int> directions   ;
     std::vector<int> displacements;
+    std::vector<int> opdirs;
 
+    // FIXME -- don't like xposing the operator directions
+    // as different to the geometrical dirs
+    // Also don't like special casing five dim.. should pass an object in template
+  Geometry(int _d)  {
+  
+      int base = (_d==5) ? 1:0;
+
+      // make coarse grid stencil for 4d , not 5d
+      if ( _d==5 ) _d=4;
+
+      npoint = 2*_d+1;
+      directions.resize(npoint);
+      displacements.resize(npoint);
+      opdirs.resize(npoint);
+      for(int d=0;d<_d;d++){
+	directions[2*d  ] = d+base;
+	directions[2*d+1] = d+base;
+	opdirs[2*d  ] = d;
+	opdirs[2*d+1] = d;
+	displacements[2*d  ] = +1;
+	displacements[2*d+1] = -1;
+      }
+      directions   [2*_d]=0;
+      displacements[2*_d]=0;
+      opdirs       [2*_d]=0;
+    }
+  
+    /*
+      // Original cleaner code
     Geometry(int _d) : dimension(_d), npoint(2*_d+1), directions(npoint), displacements(npoint) {
       for(int d=0;d<dimension;d++){
 	directions[2*d  ] = d;
@@ -22,12 +52,12 @@ namespace Grid {
       directions   [2*dimension]=0;
       displacements[2*dimension]=0;
     }
-    
     std::vector<int> GetDelta(int point) {
       std::vector<int> delta(dimension,0);
       delta[directions[point]] = displacements[point];
       return delta;
     };
+    */    
 
   };
   
@@ -50,7 +80,10 @@ namespace Grid {
     Geometry         geom;
     GridBase *       _grid; 
     CartesianStencil Stencil; 
+
     std::vector<CoarseMatrix> A;
+    std::vector<CoarseMatrix> Aslow;
+
     std::vector<siteVector,alignedAllocator<siteVector> >   comm_buf;
       
     ///////////////////////
@@ -63,7 +96,7 @@ namespace Grid {
       SimpleCompressor<siteVector> compressor;
       Stencil.HaloExchange(in,comm_buf,compressor);
 
-      //PARALLEL_FOR_LOOP
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<Grid()->oSites();ss++){
         siteVector res = zero;
 	siteVector tmp;
@@ -101,43 +134,112 @@ namespace Grid {
       _grid(&CoarseGrid),
       geom(CoarseGrid._ndimension),
       Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements),
-      A(geom.npoint,&CoarseGrid)
+	A(geom.npoint,&CoarseGrid),
+	Aslow(geom.npoint,&CoarseGrid)
     {
       comm_buf.resize(Stencil._unified_buffer_size);
     };
 
     void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,std::vector<Lattice<Fobj> > & subspace){
 
-      FineField  phi(FineGrid);
+      FineField iblock(FineGrid); // contributions from within this block
-      FineField Mphi(FineGrid);
+      FineField oblock(FineGrid); // contributions from outwith this block
-      CoarseVector Proj(Grid()); 
+
+      FineField     phi(FineGrid);
+      FineField     tmp(FineGrid);
+      FineField     zz(FineGrid); zz=zero;
+      FineField    Mphi(FineGrid);
+
+      Lattice<iScalar<vInteger> > coor(FineGrid);
+
+      CoarseVector iProj(Grid()); 
+      CoarseVector oProj(Grid()); 
       CoarseScalar InnerProd(Grid()); 
 
       // Orthogonalise the subblocks over the basis
       blockOrthogonalise(InnerProd,subspace);
+      blockProject(iProj,subspace[0],subspace);
 
       // Compute the matrix elements of linop between this orthonormal
       // set of vectors.
+      int self_stencil=-1;
+      for(int p=0;p<geom.npoint;p++){ 
+	A[p]=zero;
+	if( geom.displacements[p]==0){
+	  self_stencil=p;
+	}
+      }
+      assert(self_stencil!=-1);
+
       for(int i=0;i<nbasis;i++){
 	phi=subspace[i];
 	for(int p=0;p<geom.npoint;p++){ 
-	  int dir = geom.directions[p];
+
+	  int dir   = geom.directions[p];
+	  int opdir = geom.opdirs[p];
 	  int disp= geom.displacements[p];
 
-	  if ( disp==0 )linop.OpDiag(phi,Mphi);
+	  int block=(FineGrid->_rdimensions[dir])/(Grid()->_rdimensions[dir]);
-	  else linop.OpDir(phi,Mphi,dir,disp); 
 
-	  blockProject(Proj,Mphi,subspace);
+	  LatticeCoordinate(coor,dir);
 
-	  for(int ss=0;ss<Grid()->oSites();ss++){
+	  if ( disp==0 ){
-	    for(int j=0;j<nbasis;j++){
+	    linop.OpDiag(phi,Mphi);
-	      A[p]._odata[ss](j,i) = Proj._odata[ss](j);
+	  }
-	    }
+	  else  {
+	    linop.OpDir(phi,Mphi,opdir,disp); 
 	  }
 
-	}	
+	  ////////////////////////////////////////////////////////////////////////
+	  // Pick out contributions coming from this cell and neighbour cell
+	  ////////////////////////////////////////////////////////////////////////
+	  if ( disp==0 ) {
+	    iblock = Mphi;
+	    oblock = zero;
+	  } else if ( disp==1 ) {
+	    oblock = where(mod(coor,block)==(block-1),Mphi,zz);
+	    iblock = where(mod(coor,block)!=(block-1),Mphi,zz);
+	  } else if ( disp==-1 ) {
+	    oblock = where(mod(coor,block)==0,Mphi,zz);
+	    iblock = where(mod(coor,block)!=0,Mphi,zz);
+	  } else {
+	    assert(0);
+	  }
+
+	  blockProject(iProj,iblock,subspace);
+	  blockProject(oProj,oblock,subspace);
+	  for(int ss=0;ss<Grid()->oSites();ss++){
+	    for(int j=0;j<nbasis;j++){
+	      if( disp!= 0 ) {
+		A[p]._odata[ss](j,i) = oProj._odata[ss](j);
+	      }
+	      A[self_stencil]._odata[ss](j,i) =	A[self_stencil]._odata[ss](j,i) + iProj._odata[ss](j);
+	    }
+	  }
+	}
       }
+
+#if 0
+      ///////////////////////////
+      // test code worth preserving in if block
+      ///////////////////////////
+      std::cout<< " Computed matrix elements "<< self_stencil <<std::endl;
+      for(int p=0;p<geom.npoint;p++){
+	std::cout<< "A["<<p<<"]" << std::endl;
+	std::cout<< A[p] << std::endl;
+      }
+      std::cout<< " picking by block0 "<< self_stencil <<std::endl;
+
+      phi=subspace[0];
+      std::vector<int> bc(FineGrid->_ndimension,0);
+
+      blockPick(Grid(),phi,tmp,bc);      // Pick out a block
+      linop.Op(tmp,Mphi);                // Apply big dop
+      blockProject(iProj,Mphi,subspace); // project it and print it
+      std::cout<< " Computed matrix elements from block zero only "<<std::endl;
+      std::cout<< iProj <<std::endl;
       std::cout<<"Computed Coarse Operator"<<std::endl;
+#endif
     }
     
   };

lib/lattice/Lattice_base.h

@@ -297,11 +297,13 @@ PARALLEL_FOR_LOOP
 #include <lattice/Lattice_reduction.h>
 #include <lattice/Lattice_peekpoke.h>
 #include <lattice/Lattice_reality.h>
+#include <lattice/Lattice_comparison_utils.h>
+#include <lattice/Lattice_comparison.h>
 #include <lattice/Lattice_coordinate.h>
+#include <lattice/Lattice_where.h>
 #include <lattice/Lattice_rng.h>
 #include <lattice/Lattice_unary.h>
 #include <lattice/Lattice_transfer.h>
 
 
-
 #endif

lib/lattice/Lattice_comparison_utils.h

@@ -141,7 +141,5 @@ namespace Grid {
     }
 }
 
-#include <lattice/Lattice_comparison.h>
-#include <lattice/Lattice_where.h>
 
 #endif

lib/lattice/Lattice_coordinate.h

@@ -3,20 +3,8 @@
 
 namespace Grid {
 
-  /*
-  depbase=`echo Grid_main.o | sed 's|[^/]*$|.deps/&|;s|\.o$||'`;\
-        icpc -DHAVE_CONFIG_H -I. -I../lib    -I../lib -mmic -O3 -std=c++11 -fopenmp -MT Grid_main.o -MD -MP -MF $depbase.Tpo -c -o Grid_main.o Grid_main.cc &&\
-        mv -f $depbase.Tpo $depbase.Po
-	  ../lib/lattice/Grid_lattice_coordinate.h(25): error: no suitable user-defined conversion from "vector_type" to "const Grid::iScalar<Grid::iScalar<Grid::iScalar<Grid::vInteger>>>" exists
-    l._odata[o]=vI;
-                    ^
-          detected during instantiation of "void Grid::LatticeCoordinate(Grid::Lattice<iobj> &, int) [with iobj=Grid::QCD::vTInteger]" at line 283 of "Grid_main.cc"
-
-	    compilation aborted for Grid_main.cc (code 2)
-*/
     template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
     {
-      typedef typename iobj::scalar_object scalar_object;
       typedef typename iobj::scalar_type scalar_type;
       typedef typename iobj::vector_type vector_type;
 
@@ -33,7 +21,7 @@ namespace Grid {
 	  mergebuf[i]=(Integer)gcoor[mu];
 	}
 	merge<vector_type,scalar_type>(vI,mergebuf);
-	l._odata[o]._internal._internal._internal=vI;
+	l._odata[o]=vI;
       }
     };
 

lib/lattice/Lattice_transfer.h

@@ -166,9 +166,10 @@ template<class vobj,class CComplex>
 inline void blockNormalise(Lattice<CComplex> &ip,Lattice<vobj> &fineX)
 {
   GridBase *coarse = ip._grid;
+  Lattice<vobj> zz(fineX._grid); zz=zero;
   blockInnerProduct(ip,fineX,fineX);
   ip = rsqrt(ip);
-  blockZAXPY(fineX,ip,fineX,fineX);
+  blockZAXPY(fineX,ip,fineX,zz);
 }
 // useful in multigrid project;
 // Generic name : Coarsen?
@@ -205,6 +206,26 @@ inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
   return;
 }
 
+template<class vobj>
+inline void blockPick(GridBase *coarse,const Lattice<vobj> &unpicked,Lattice<vobj> &picked,std::vector<int> coor)
+{
+  GridBase * fine = unpicked._grid;
+
+  Lattice<vobj> zz(fine);
+  Lattice<iScalar<vInteger> > fcoor(fine);
+
+  zz = zero;
+
+  picked = unpicked;
+  for(int d=0;d<fine->_ndimension;d++){
+    LatticeCoordinate(fcoor,d);
+    int block= fine->_rdimensions[d] / coarse->_rdimensions[d];
+    int lo   = (coor[d])*block;
+    int hi   = (coor[d]+1)*block;
+    picked = where( (fcoor<hi) , picked, zz);
+    picked = where( (fcoor>=lo), picked, zz);
+  }
+}
 
 template<class vobj,class CComplex>
 inline void blockOrthogonalise(Lattice<CComplex> &ip,std::vector<Lattice<vobj> > &Basis)

lib/lattice/Lattice_unary.h

@@ -26,7 +26,48 @@ PARALLEL_FOR_LOOP
     }
     return ret;
   }
+  template<class obj> Lattice<obj> sin(const Lattice<obj> &rhs){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=sin(rhs._odata[ss]);
+    }
+    return ret;
+  }
+  template<class obj> Lattice<obj> cos(const Lattice<obj> &rhs){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=cos(rhs._odata[ss]);
+    }
+    return ret;
+  }
 
 
+  template<class obj> Lattice<obj> pow(const Lattice<obj> &rhs,RealD y){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=pow(rhs._odata[ss],y);
+    }
+    return ret;
+  }
+  template<class obj> Lattice<obj> mod(const Lattice<obj> &rhs,Integer y){
+    Lattice<obj> ret(rhs._grid);
+    ret.checkerboard = rhs.checkerboard;
+    conformable(ret,rhs);
+PARALLEL_FOR_LOOP
+    for(int ss=0;ss<rhs._grid->oSites();ss++){
+      ret._odata[ss]=mod(rhs._odata[ss],y);
+    }
+    return ret;
+  }
+
 }
 #endif

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -85,6 +85,7 @@ void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGau
 void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp)
 {
   assert( (disp==1)||(disp==-1) );
+  assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t;
 
   WilsonCompressor compressor(DaggerNo);
   Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
@@ -93,6 +94,9 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int
 
   int dirdisp = dir+skip*4;
 
+  assert(dirdisp<=7);
+  assert(dirdisp>=0);
+
 PARALLEL_FOR_LOOP
   for(int ss=0;ss<Umu._grid->oSites();ss++){
     for(int s=0;s<Ls;s++){

lib/simd/Grid_vector_unops.h

@@ -35,6 +35,14 @@ namespace Grid {
     }
   };
 
+  template<class scalar> struct ModIntFunctor {
+    Integer y;
+  ModIntFunctor(Integer _y) : y(_y) {};
+    scalar operator()(const scalar &a)  const {
+      return Integer(a)%y;
+    }
+  };
+
   template < class S, class V > 
   inline Grid_simd<S,V> sqrt(const Grid_simd<S,V> &r) {
     return SimdApply(SqrtRealFunctor<S>(),r);
@@ -55,6 +63,10 @@ namespace Grid {
   inline Grid_simd<S,V> pow(const Grid_simd<S,V> &r,double y) {
     return SimdApply(PowRealFunctor<S>(y),r);
   }
+  template < class S, class V > 
+  inline Grid_simd<S,V> mod(const Grid_simd<S,V> &r,Integer y) {
+    return SimdApply(ModIntFunctor<S>(y),r);
+  }
 
 }
 #endif

lib/tensors/Tensor_class.h

@@ -90,10 +90,10 @@ public:
   operator ComplexD () const { return(TensorRemove(_internal)); };
   operator RealD () const { return(real(TensorRemove(_internal))); }
   
-  // convert from a something to a scalar
+  // convert from a something to a scalar via constructor of something arg
-  template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iScalar<vtype>
+  template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline iScalar<vtype> operator = (T arg)
     { 
-      _internal = vtype(arg);
+      _internal = arg;
       return *this;
     }
 
@@ -316,7 +316,8 @@ public:
 	stream<<o._internal[i][j];
 	if (i<N-1)	stream<<",";
       }
-      stream<<"}\n\t\t";
+      stream<<"}";
+      if(i!=N-1) stream<<"\n\t\t";
     }
     stream<<"}";
     return stream;

lib/tensors/Tensor_unary.h

@@ -27,11 +27,40 @@ template<class obj,int N> inline auto func(const iMatrix<obj,N> &z) -> iMatrix<o
     return ret;\
 }
 
+
+#define BINARY_RSCALAR(func,scal)					\
+template<class obj> inline iScalar<obj> func(const iScalar<obj> &z,scal y)	\
+{\
+    iScalar<obj> ret;\
+    ret._internal = func(z._internal,y);	\
+    return ret;\
+}\
+ template<class obj,int N> inline iVector<obj,N> func(const iVector<obj,N> &z,scal y) \
+{\
+    iVector<obj,N> ret;\
+    for(int c1=0;c1<N;c1++){\
+      ret._internal[c1] = func(z._internal[c1],y);	\
+    }\
+    return ret;\
+}\
+ template<class obj,int N> inline  iMatrix<obj,N> func(const iMatrix<obj,N> &z, scal y) \
+{\
+    iMatrix<obj,N> ret;\
+    for(int c1=0;c1<N;c1++){\
+    for(int c2=0;c2<N;c2++){\
+      ret._internal[c1][c2] = func(z._internal[c1][c2],y);	\
+    }}\
+    return ret;\
+}
+
 UNARY_REAL(sqrt);
 UNARY_REAL(rsqrt);
 UNARY_REAL(sin);
 UNARY_REAL(cos);
 
+BINARY_RSCALAR(mod,Integer);
+BINARY_RSCALAR(pow,RealD);
+
 
 }
 #endif

scripts/linecount

@@ -1,3 +1,12 @@
 #!/bin/sh 
 
-wc -l lib/*.h lib/*/*.h lib/*/*/*.h lib/*.cc lib/*/*.cc lib/*/*/*.cc tests/*.cc benchmarks/*.cc lib/*/*/*/*.cc lib/*/*/*/*.h
+wc -l lib/*.h lib/*/*.h lib/*/*/*.h lib/*.cc lib/*/*.cc lib/*/*/*.cc tests/*.cc benchmarks/*.cc lib/*/*/*/*.cc lib/*/*/*/*.h
+
+echo "lib"
+wc -l lib/*.h lib/*/*.h lib/*/*/*.h lib/*.cc lib/*/*.cc lib/*/*/*.cc lib/*/*/*/*.cc lib/*/*/*/*.h | grep total
+
+echo "tests"
+wc -l tests/*.cc | grep total
+
+echo "benchmarks"
+wc -l benchmarks/*.cc | grep total

tests/Test_cayley_coarsen_support.cc

@@ -23,7 +23,7 @@ int main (int argc, char ** argv)
 
 
 
-  const int Ls=8;
+  const int Ls=4;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
@@ -50,13 +50,18 @@ int main (int argc, char ** argv)
   LatticeFermion    tmp(FGrid);
   LatticeFermion    err(FGrid);
   LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
-
+#if 0
   std::vector<LatticeColourMatrix> U(4,UGrid);
-
+  Umu=zero;
-  for(int mu=0;mu<Nd;mu++){
+  Complex cone(1.0,0.0);
-    U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+  for(int nn=0;nn<Nd;nn++){
+    if(1) {
+      if (nn!=0) { U[nn]=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl; }
+      else       { U[nn] = cone;std::cout << "unit gauge field in dir "<<nn<<std::endl; }
+    }
+    pokeIndex<LorentzIndex>(Umu,U[nn],nn);
   }
-  
+#endif  
   RealD mass=0.5;
   RealD M5=1.8;
 
@@ -75,7 +80,7 @@ int main (int argc, char ** argv)
   err = result-ref;
   std::cout<<"Error "<<norm2(err)<<std::endl;
 
-  const int nbasis = 8;
+  const int nbasis = 2;
   std::vector<LatticeFermion> subspace(nbasis,FGrid);
   
   for(int b=0;b<nbasis;b++){

tests/Test_cshift.cc

@@ -34,7 +34,6 @@ int main (int argc, char ** argv)
   }
 
 
-
   TComplex cm;
   
   for(int dir=0;dir<4;dir++){