Const correctness patch

by Christoph. Differs from the one by Rudy in BFM since it vectorises the twisted
4D solves in pairs.

.gitignore

@@ -83,6 +83,7 @@ ltmain.sh
 .Trashes
 ehthumbs.db
 Thumbs.db
+.dirstamp
 
 # build directory #
 ###################
@@ -97,11 +98,8 @@ build.sh
 
 # Eigen source #
 ################
-lib/Eigen/*
-
-# FFTW source #
-################
-lib/fftw/*
+Grid/Eigen
+Eigen/*
 
 # libtool macros #
 ##################
@@ -112,21 +110,7 @@ m4/libtool.m4
 ################
 gh-pages/
 
-# Buck files #
-##############
-.buck*
-buck-out
-BUCK
-make-bin-BUCK.sh
-
 # generated sources #
 #####################
-lib/qcd/spin/gamma-gen/*.h
-lib/qcd/spin/gamma-gen/*.cc
-lib/version.h
-
-# vs code editor files #
-########################
-.vscode/
-.vscode/settings.json
-settings.json
+Grid/qcd/spin/gamma-gen/*.h
+Grid/qcd/spin/gamma-gen/*.cc

.travis.yml

@@ -9,6 +9,11 @@ matrix:
     - os:        osx
       osx_image: xcode8.3
       compiler: clang
+      env: PREC=single
+    - os:        osx
+      osx_image: xcode8.3
+      compiler: clang
+      env: PREC=double
       
 before_install:
     - export GRIDDIR=`pwd`
@@ -16,9 +21,11 @@ before_install:
     - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export PATH="${GRIDDIR}/clang/bin:${PATH}"; fi
     - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export LD_LIBRARY_PATH="${GRIDDIR}/clang/lib:${LD_LIBRARY_PATH}"; fi
     - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi
-    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install libmpc; fi
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install libmpc openssl; fi
     
 install:
+    - export CWD=`pwd`
+    - echo $CWD
     - export CC=$CC$VERSION
     - export CXX=$CXX$VERSION
     - echo $PATH
@@ -31,17 +38,24 @@ install:
     - which $CXX
     - $CXX --version
     - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export LDFLAGS='-L/usr/local/lib'; fi
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export EXTRACONF='--with-openssl=/usr/local/opt/openssl'; fi
     
 script:
     - ./bootstrap.sh
     - mkdir build
     - cd build
-    - ../configure --enable-precision=single --enable-simd=SSE4 --enable-comms=none
+    - mkdir lime
+    - cd lime
+    - mkdir build
+    - cd build
+    - wget http://usqcd-software.github.io/downloads/c-lime/lime-1.3.2.tar.gz
+    - tar xf lime-1.3.2.tar.gz
+    - cd lime-1.3.2
+    - ./configure --prefix=$CWD/build/lime/install
+    - make -j4
+    - make install
+    - cd $CWD/build
+    - ../configure --enable-precision=$PREC --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install ${EXTRACONF}
     - make -j4 
     - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
-    - echo make clean
-    - ../configure --enable-precision=double --enable-simd=SSE4 --enable-comms=none
-    - make -j4
-    - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
     - make check
-

Grid/GridCore.h

@@ -48,6 +48,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/serialisation/Serialisation.h>
 #include <Grid/threads/Threads.h>
 #include <Grid/util/Util.h>
+#include <Grid/util/Sha.h>
 #include <Grid/communicator/Communicator.h> 
 #include <Grid/cartesian/Cartesian.h>    
 #include <Grid/tensors/Tensors.h>      

Grid/Grid_Eigen_Dense.h

@@ -1,4 +1,9 @@
 #pragma once
+// Force Eigen to use MKL if Grid has been configured with --enable-mkl
+#ifdef USE_MKL
+#define EIGEN_USE_MKL_ALL
+#endif
+
 #if defined __GNUC__
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"

Grid/Makefile.am

@@ -21,6 +21,32 @@ if BUILD_HDF5
   extra_headers+=serialisation/Hdf5Type.h
 endif
 
+all: version-cache
+
+version-cache:
+	@if [ `git status --porcelain | grep -v '??' | wc -l` -gt 0 ]; then\
+		a="uncommited changes";\
+	else\
+		a="clean";\
+	fi;\
+	echo "`git log -n 1 --format=format:"#define GITHASH \\"%H:%d $$a\\"%n" HEAD`" > vertmp;\
+	if [ -e version-cache ]; then\
+		d=`diff vertmp version-cache`;\
+		if [ "$${d}" != "" ]; then\
+			mv vertmp version-cache;\
+			rm -f Version.h;\
+		fi;\
+	else\
+		mv vertmp version-cache;\
+		rm -f Version.h;\
+	fi;\
+	rm -f vertmp
+
+Version.h:
+	cp version-cache Version.h
+
+.PHONY: version-cache
+
 #
 # Libraries
 #
@@ -30,8 +56,8 @@ include Eigen.inc
 lib_LIBRARIES = libGrid.a
 
 CCFILES += $(extra_sources)
-HFILES  += $(extra_headers)
+HFILES  += $(extra_headers) Config.h Version.h
 
 libGrid_a_SOURCES              = $(CCFILES)
-libGrid_adir                   = $(pkgincludedir)
-nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) Config.h
+libGrid_adir                   = $(includedir)/Grid
+nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) $(eigen_unsupp_files)

Grid/algorithms/LinearOperator.h

@@ -51,7 +51,7 @@ namespace Grid {
 
       virtual void Op     (const Field &in, Field &out) = 0; // Abstract base
       virtual void AdjOp  (const Field &in, Field &out) = 0; // Abstract base
-      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0;
+      virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) = 0;
       virtual void HermOp(const Field &in, Field &out)=0;
     };
 
@@ -309,36 +309,59 @@ namespace Grid {
       class SchurStaggeredOperator :  public SchurOperatorBase<Field> {
     protected:
       Matrix &_Mat;
+      Field tmp;
+      RealD mass;
+      double tMpc;
+      double tIP;
+      double tMeo;
+      double taxpby_norm;
+      uint64_t ncall;
     public:
-      SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){};
+      void Report(void)
+      {
+	std::cout << GridLogMessage << " HermOpAndNorm.Mpc "<< tMpc/ncall<<" usec "<<std::endl;
+	std::cout << GridLogMessage << " HermOpAndNorm.IP  "<< tIP /ncall<<" usec "<<std::endl;
+	std::cout << GridLogMessage << " Mpc.MeoMoe        "<< tMeo/ncall<<" usec "<<std::endl;
+	std::cout << GridLogMessage << " Mpc.axpby_norm    "<< taxpby_norm/ncall<<" usec "<<std::endl;
+      }
+      SchurStaggeredOperator (Matrix &Mat): _Mat(Mat), tmp(_Mat.RedBlackGrid()) 
+      { 
+	assert( _Mat.isTrivialEE() );
+	mass = _Mat.Mass();
+	tMpc=0;
+	tIP =0;
+        tMeo=0;
+        taxpby_norm=0;
+	ncall=0;
+      }
       virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-	GridLogIterative.TimingMode(1);
-	std::cout << GridLogIterative << " HermOpAndNorm "<<std::endl;
+	ncall++;
+	tMpc-=usecond();
 	n2 = Mpc(in,out);
-	std::cout << GridLogIterative << " HermOpAndNorm.Mpc "<<std::endl;
+	tMpc+=usecond();
+	tIP-=usecond();
 	ComplexD dot= innerProduct(in,out);
-	std::cout << GridLogIterative << " HermOpAndNorm.innerProduct "<<std::endl;
+	tIP+=usecond();
 	n1 = real(dot);
       }
       virtual void HermOp(const Field &in, Field &out){
-	std::cout << GridLogIterative << " HermOp "<<std::endl;
-	Mpc(in,out);
+	ncall++;
+	tMpc-=usecond();
+	_Mat.Meooe(in,out);
+	_Mat.Meooe(out,tmp);
+	tMpc+=usecond();
+	taxpby_norm-=usecond();
+	axpby(out,-1.0,mass*mass,tmp,in);
+	taxpby_norm+=usecond();
       }
       virtual  RealD Mpc      (const Field &in, Field &out) {
-	Field tmp(in._grid);
-	Field tmp2(in._grid);
-
-	std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
-	_Mat.Mooee(in,out);
-	_Mat.Mooee(out,tmp);
-	std::cout << GridLogIterative << " HermOp.MooeeMooee "<<std::endl;
-
+	tMeo-=usecond();
 	_Mat.Meooe(in,out);
-	_Mat.Meooe(out,tmp2);
-	std::cout << GridLogIterative << " HermOp.MeooeMeooe "<<std::endl;
-
-	RealD nn=axpy_norm(out,-1.0,tmp2,tmp);
-	std::cout << GridLogIterative << " HermOp.axpy_norm "<<std::endl;
+	_Mat.Meooe(out,tmp);
+	tMeo+=usecond();
+	taxpby_norm-=usecond();
+	RealD nn=axpby_norm(out,-1.0,mass*mass,tmp,in);
+	taxpby_norm+=usecond();
 	return nn;
       }
       virtual  RealD MpcDag   (const Field &in, Field &out){

Grid/algorithms/iterative/ConjugateGradient.h

@@ -54,6 +54,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
 
   void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
 
+
     psi.checkerboard = src.checkerboard;
     conformable(psi, src);
 
@@ -69,7 +70,6 @@ class ConjugateGradient : public OperatorFunction<Field> {
 
     
     Linop.HermOpAndNorm(psi, mmp, d, b);
-    
 
     r = src - mmp;
     p = r;
@@ -96,38 +96,44 @@ class ConjugateGradient : public OperatorFunction<Field> {
               << "ConjugateGradient: k=0 residual " << cp << " target " << rsq << std::endl;
 
     GridStopWatch LinalgTimer;
+    GridStopWatch InnerTimer;
+    GridStopWatch AxpyNormTimer;
+    GridStopWatch LinearCombTimer;
     GridStopWatch MatrixTimer;
     GridStopWatch SolverTimer;
 
     SolverTimer.Start();
     int k;
-    for (k = 1; k <= MaxIterations; k++) {
+    for (k = 1; k <= MaxIterations*1000; k++) {
       c = cp;
 
       MatrixTimer.Start();
-      Linop.HermOpAndNorm(p, mmp, d, qq);
+      Linop.HermOp(p, mmp);
       MatrixTimer.Stop();
 
       LinalgTimer.Start();
-      //  RealD    qqck = norm2(mmp);
-      //  ComplexD dck  = innerProduct(p,mmp);
 
+      InnerTimer.Start();
+      ComplexD dc  = innerProduct(p,mmp);
+      InnerTimer.Stop();
+      d = dc.real();
       a = c / d;
-      b_pred = a * (a * qq - d) / c;
 
+      AxpyNormTimer.Start();
       cp = axpy_norm(r, -a, mmp, r);
+      AxpyNormTimer.Stop();
       b = cp / c;
 
-      // Fuse these loops ; should be really easy
-      psi = a * p + psi;
-      p = p * b + r;
-
+      LinearCombTimer.Start();
+      parallel_for(int ss=0;ss<src._grid->oSites();ss++){
+	vstream(psi[ss], a      *  p[ss] + psi[ss]);
+	vstream(p  [ss], b      *  p[ss] + r[ss]);
+      }
+      LinearCombTimer.Stop();
       LinalgTimer.Stop();
 
       std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k
                 << " residual " << cp << " target " << rsq << std::endl;
-      std::cout << GridLogDebug << "a = "<< a << " b_pred = "<< b_pred << "  b = "<< b << std::endl;
-      std::cout << GridLogDebug << "qq = "<< qq << " d = "<< d << "  c = "<< c << std::endl;
 
       // Stopping condition
       if (cp <= rsq) {
@@ -144,10 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
 	std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
 	std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
 
-        std::cout << GridLogMessage << "Time breakdown "<<std::endl;
-	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
+	std::cout << GridLogPerformance << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 

Grid/algorithms/iterative/ConjugateGradientMultiShift.h

@@ -43,6 +43,7 @@ namespace Grid {
 public:                                                
     RealD   Tolerance;
     Integer MaxIterations;
+    Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
     int verbose;
     MultiShiftFunction shifts;
 
@@ -163,7 +164,16 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
   for(int s=0;s<nshift;s++) {
     axpby(psi[s],0.,-bs[s]*alpha[s],src,src);
   }
-  
+ 
+  ///////////////////////////////////////
+  // Timers
+  ///////////////////////////////////////
+  GridStopWatch AXPYTimer;
+  GridStopWatch ShiftTimer;
+  GridStopWatch QRTimer;
+  GridStopWatch MatrixTimer;
+  GridStopWatch SolverTimer;
+  SolverTimer.Start();
   
   // Iteration loop
   int k;
@@ -171,7 +181,9 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
   for (k=1;k<=MaxIterations;k++){
     
     a = c /cp;
+    AXPYTimer.Start();
     axpy(p,a,p,r);
+    AXPYTimer.Stop();
     
     // Note to self - direction ps is iterated seperately
     // for each shift. Does not appear to have any scope
@@ -180,6 +192,7 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
     // However SAME r is used. Could load "r" and update
     // ALL ps[s]. 2/3 Bandwidth saving
     // New Kernel: Load r, vector of coeffs, vector of pointers ps
+    AXPYTimer.Start();
     for(int s=0;s<nshift;s++){
       if ( ! converged[s] ) { 
 	if (s==0){
@@ -190,22 +203,34 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
 	}
       }
     }
+    AXPYTimer.Stop();
     
     cp=c;
+    MatrixTimer.Start();  
+    //Linop.HermOpAndNorm(p,mmp,d,qq); // d is used
+    // The below is faster on KNL
+    Linop.HermOp(p,mmp); 
+    d=real(innerProduct(p,mmp));
     
-    Linop.HermOpAndNorm(p,mmp,d,qq);
+    MatrixTimer.Stop();  
+
+    AXPYTimer.Start();
     axpy(mmp,mass[0],p,mmp);
+    AXPYTimer.Stop();
     RealD rn = norm2(p);
     d += rn*mass[0];
     
     bp=b;
     b=-cp/d;
     
+    AXPYTimer.Start();
     c=axpy_norm(r,b,mmp,r);
+    AXPYTimer.Stop();
 
     // 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];
@@ -215,6 +240,7 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
 	bs[s] = b*z[s][iz]/z0; // NB sign  rel to Mike
       }
     }
+    ShiftTimer.Stop();
     
     for(int s=0;s<nshift;s++){
       int ss = s;
@@ -257,6 +283,9 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
     
     if ( all_converged ){
 
+    SolverTimer.Stop();
+
+
       std::cout<<GridLogMessage<< "CGMultiShift: All shifts have converged iteration "<<k<<std::endl;
       std::cout<<GridLogMessage<< "CGMultiShift: Checking solutions"<<std::endl;
       
@@ -269,8 +298,19 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
 	RealD cn = norm2(src);
 	std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<<std::sqrt(rn/cn)<<std::endl;
       }
+
+      std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
+      std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
+      std::cout << GridLogMessage << "\tAXPY    " << AXPYTimer.Elapsed()     <<std::endl;
+      std::cout << GridLogMessage << "\tMarix    " << MatrixTimer.Elapsed()     <<std::endl;
+      std::cout << GridLogMessage << "\tShift    " << ShiftTimer.Elapsed()     <<std::endl;
+
+      IterationsToComplete = k;	
+
       return;
     }
+
+   
   }
   // ugly hack
   std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;

Grid/algorithms/iterative/Deflation.h

@@ -30,22 +30,23 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 namespace Grid { 
 
-struct ZeroGuesser {
+template<class Field>
+class ZeroGuesser: public LinearFunction<Field> {
 public:
-  template<class Field> 
-  void operator()(const Field &src,Field &guess) { guess = Zero(); };
+  virtual void operator()(const Field &src, Field &guess) { guess = zero; };
 };
-struct SourceGuesser {
+
+template<class Field>
+class SourceGuesser: public LinearFunction<Field> {
 public:
-  template<class Field> 
-  void operator()(const Field &src,Field &guess) { guess = src; };
+  virtual void operator()(const Field &src, Field &guess) { guess = src; };
 };
 
 ////////////////////////////////
 // Fine grid deflation
 ////////////////////////////////
 template<class Field>
-struct DeflatedGuesser {
+class DeflatedGuesser: public LinearFunction<Field> {
 private:
   const std::vector<Field> &evec;
   const std::vector<RealD> &eval;
@@ -54,7 +55,7 @@ public:
 
   DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval) : evec(_evec), eval(_eval) {};
 
-  void operator()(const Field &src,Field &guess) { 
+  virtual void operator()(const Field &src,Field &guess) {
     guess = zero;
     assert(evec.size()==eval.size());
     auto N = evec.size();
@@ -62,11 +63,12 @@ public:
       const Field& tmp = evec[i];
       axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess);
     }
+    guess.checkerboard = src.checkerboard;
   }
 };
 
 template<class FineField, class CoarseField>
-class LocalCoherenceDeflatedGuesser {
+class LocalCoherenceDeflatedGuesser: public LinearFunction<FineField> {
 private:
   const std::vector<FineField>   &subspace;
   const std::vector<CoarseField> &evec_coarse;
@@ -92,6 +94,7 @@ public:
       axpy(guess_coarse,TensorRemove(innerProduct(tmp,src_coarse)) / eval_coarse[i],tmp,guess_coarse);
     }
     blockPromote(guess_coarse,guess,subspace);
+    guess.checkerboard = src.checkerboard;
   };
 };
 

Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -57,8 +57,9 @@ void basisRotate(std::vector<Field> &basis,Eigen::MatrixXd& Qt,int j0, int j1, i
       
   parallel_region
   {
-    std::vector < vobj > B(Nm); // Thread private
-        
+
+    std::vector < vobj , commAllocator<vobj> > B(Nm); // Thread private
+       
     parallel_for_internal(int ss=0;ss < grid->oSites();ss++){
       for(int j=j0; j<j1; ++j) B[j]=0.;
       

Grid/algorithms/iterative/LocalCoherenceLanczos.h

@@ -285,9 +285,11 @@ public:
   };
 
   void Orthogonalise(void ) {
-    CoarseScalar InnerProd(_CoarseGrid); 
-    blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
-    blockOrthogonalise(InnerProd,subspace);std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
+    CoarseScalar InnerProd(_CoarseGrid);
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
+    std::cout << GridLogMessage <<" Gramm-Schmidt pass 2"<<std::endl;
+    blockOrthogonalise(InnerProd,subspace);
   };
 
   template<typename T>  static RealD normalise(T& v) 
@@ -333,7 +335,7 @@ public:
     // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
     //////////////////////////////////////////////////////////////////////////////////////////////////
     Chebyshev<FineField>                          ChebySmooth(cheby_smooth);
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (ChebySmooth,_FineOp,_subspace);
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (ChebySmooth,_FineOp,subspace);
     ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);
 
     for(int k=0;k<evec_coarse.size();k++){
@@ -374,14 +376,14 @@ public:
 		  RealD MaxIt, RealD betastp, int MinRes)
   {
     Chebyshev<FineField>                          Cheby(cheby_op);
-    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,_subspace);
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,_subspace);
+    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace);
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace);
     //////////////////////////////////////////////////////////////////////////////////////////////////
     // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
     //////////////////////////////////////////////////////////////////////////////////////////////////
 
     Chebyshev<FineField>                                           ChebySmooth(cheby_smooth);
-    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,_subspace,relax);
+    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);
 
     evals_coarse.resize(Nm);
     evec_coarse.resize(Nm,_CoarseGrid);

Grid/algorithms/iterative/SchurRedBlack.h

@@ -95,20 +95,30 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
       CBfactorise=0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>
@@ -150,9 +160,12 @@ namespace Grid {
       // Call the red-black solver
       //////////////////////////////////////////////////////////////
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
-      guess(src_o,sol_o);
+      guess(src_o, sol_o);
+      Mtmp = sol_o;
       _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called  the Mpc solver" <<std::endl;
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -167,11 +180,15 @@ namespace Grid {
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl;
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
-      std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
+        std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
   template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
@@ -184,18 +201,28 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0)  :  _HermitianRBSolver(HermitianRBSolver) 
+  SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0, const bool initSubGuess = false)  :  _HermitianRBSolver(HermitianRBSolver) 
   { 
     CBfactorise=cb;
+    subtractGuess(initSubGuess);
   };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>
@@ -236,7 +263,10 @@ namespace Grid {
       //////////////////////////////////////////////////////////////
       std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
       guess(src_o,sol_o);
+      Mtmp = sol_o;
       _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -249,12 +279,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
 
@@ -267,20 +301,30 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
-      CBfactorise=0;
+      CBfactorise = 0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix,class Guesser>
@@ -322,8 +366,11 @@ namespace Grid {
       std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
 //      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
       guess(src_o,tmp);
+      Mtmp = tmp;
       _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
-      _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)      tmp = tmp-Mtmp;
+      _Matrix.MooeeInv(tmp,sol_o);       assert(  sol_o.checkerboard   ==Odd);
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -336,12 +383,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -352,20 +403,30 @@ namespace Grid {
   private:
     LinearFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
       CBfactorise=0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
-      ZeroGuesser guess;
+      ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
     template<class Matrix, class Guesser>
@@ -408,7 +469,10 @@ namespace Grid {
 //      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
 //      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
       guess(src_o,tmp);
-      _HermitianRBSolver(src_o,tmp);  assert(tmp.checkerboard==Odd);
+      Mtmp = tmp;
+      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)      tmp = tmp-Mtmp;
       _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
 
       ///////////////////////////////////////////////////
@@ -422,12 +486,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout << GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid " << std::sqrt(nr / ns) << " nr " << nr << " ns " << ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
 

Grid/communicator/SharedMemoryMPI.cc

@@ -132,7 +132,6 @@ int Log2Size(int TwoToPower,int MAXLOG2)
 }
 void GlobalSharedMemory::OptimalCommunicator(const std::vector<int> &processors,Grid_MPI_Comm & optimal_comm)
 {
-#undef HYPERCUBE 
 #ifdef HYPERCUBE
   ////////////////////////////////////////////////////////////////
   // Assert power of two shm_size.
@@ -175,7 +174,7 @@ void GlobalSharedMemory::OptimalCommunicator(const std::vector<int> &processors,
 
   std::string hname(name);
   std::cout << "hostname "<<hname<<std::endl;
-  std::cout << "R " << R << " I " << I << " N "<< N<<
+  std::cout << "R " << R << " I " << I << " N "<< N
             << " hypercoor 0x"<<std::hex<<hypercoor<<std::dec<<std::endl;
 
   //////////////////////////////////////////////////////////////////

Grid/lattice/Lattice_arith.h

@@ -244,19 +244,11 @@ namespace Grid {
 
   template<class sobj,class vobj> strong_inline
   RealD axpy_norm(Lattice<vobj> &ret,sobj a,const Lattice<vobj> &x,const Lattice<vobj> &y){
-    ret.checkerboard = x.checkerboard;
-    conformable(ret,x);
-    conformable(x,y);
-    axpy(ret,a,x,y);
-    return norm2(ret);
+    return axpy_norm_fast(ret,a,x,y);
   }
   template<class sobj,class vobj> strong_inline
   RealD axpby_norm(Lattice<vobj> &ret,sobj a,sobj b,const Lattice<vobj> &x,const Lattice<vobj> &y){
-    ret.checkerboard = x.checkerboard;
-    conformable(ret,x);
-    conformable(x,y);
-    axpby(ret,a,b,x,y);
-    return norm2(ret); // FIXME implement parallel norm in ss loop
+    return axpby_norm_fast(ret,a,b,x,y);
   }
 
 }

Grid/lattice/Lattice_reduction.h

@@ -33,41 +33,94 @@ namespace Grid {
   // Deterministic Reduction operations
   ////////////////////////////////////////////////////////////////////////////////////////////////////
 template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
-  ComplexD nrm = innerProduct(arg,arg);
+  auto nrm = innerProduct(arg,arg);
   return std::real(nrm); 
 }
 
 // Double inner product
 template<class vobj>
-inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
+inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right)
 {
   typedef typename vobj::scalar_type scalar_type;
   typedef typename vobj::vector_typeD vector_type;
-  scalar_type  nrm;
-  
   GridBase *grid = left._grid;
-  
-  std::vector<vector_type,alignedAllocator<vector_type> > sumarray(grid->SumArraySize());
-  
+  const int pad = 8;
+
+  ComplexD  inner;
+  Vector<ComplexD> sumarray(grid->SumArraySize()*pad);
+
   parallel_for(int thr=0;thr<grid->SumArraySize();thr++){
     int nwork, mywork, myoff;
     GridThread::GetWork(left._grid->oSites(),thr,mywork,myoff);
     
-    decltype(innerProductD(left._odata[0],right._odata[0])) vnrm=zero; // private to thread; sub summation
+    decltype(innerProductD(left._odata[0],right._odata[0])) vinner=zero; // private to thread; sub summation
     for(int ss=myoff;ss<mywork+myoff; ss++){
-      vnrm = vnrm + innerProductD(left._odata[ss],right._odata[ss]);
+      vinner = vinner + innerProductD(left._odata[ss],right._odata[ss]);
     }
-    sumarray[thr]=TensorRemove(vnrm) ;
+    // All threads sum across SIMD; reduce serial work at end
+    // one write per cacheline with streaming store
+    ComplexD tmp = Reduce(TensorRemove(vinner)) ;
+    vstream(sumarray[thr*pad],tmp);
   }
   
-  vector_type vvnrm; vvnrm=zero;  // sum across threads
+  inner=0.0;
   for(int i=0;i<grid->SumArraySize();i++){
-    vvnrm = vvnrm+sumarray[i];
+    inner = inner+sumarray[i*pad];
   } 
-  nrm = Reduce(vvnrm);// sum across simd
-  right._grid->GlobalSum(nrm);
-  return nrm;
+  right._grid->GlobalSum(inner);
+  return inner;
 }
+
+/////////////////////////
+// Fast axpby_norm
+// z = a x + b y
+// return norm z
+/////////////////////////
+template<class sobj,class vobj> strong_inline RealD 
+axpy_norm_fast(Lattice<vobj> &z,sobj a,const Lattice<vobj> &x,const Lattice<vobj> &y) 
+{
+  sobj one(1.0);
+  return axpby_norm_fast(z,a,one,x,y);
+}
+
+template<class sobj,class vobj> strong_inline RealD 
+axpby_norm_fast(Lattice<vobj> &z,sobj a,sobj b,const Lattice<vobj> &x,const Lattice<vobj> &y) 
+{
+  const int pad = 8;
+  z.checkerboard = x.checkerboard;
+  conformable(z,x);
+  conformable(x,y);
+
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_typeD vector_type;
+  RealD  nrm;
+  
+  GridBase *grid = x._grid;
+  
+  Vector<RealD> sumarray(grid->SumArraySize()*pad);
+  
+  parallel_for(int thr=0;thr<grid->SumArraySize();thr++){
+    int nwork, mywork, myoff;
+    GridThread::GetWork(x._grid->oSites(),thr,mywork,myoff);
+    
+    // private to thread; sub summation
+    decltype(innerProductD(z._odata[0],z._odata[0])) vnrm=zero; 
+    for(int ss=myoff;ss<mywork+myoff; ss++){
+      vobj tmp = a*x._odata[ss]+b*y._odata[ss];
+      vnrm = vnrm + innerProductD(tmp,tmp);
+      vstream(z._odata[ss],tmp);
+    }
+    vstream(sumarray[thr*pad],real(Reduce(TensorRemove(vnrm)))) ;
+  }
+  
+  nrm = 0.0; // sum across threads; linear in thread count but fast
+  for(int i=0;i<grid->SumArraySize();i++){
+    nrm = nrm+sumarray[i*pad];
+  } 
+  z._grid->GlobalSum(nrm);
+  return nrm; 
+}
+
  
 template<class Op,class T1>
 inline auto sum(const LatticeUnaryExpression<Op,T1> & expr)
@@ -221,6 +274,115 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   }
 }
 
+template<class vobj>
+static void mySliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
+{
+  // std::cout << GridLogMessage << "Start mySliceInnerProductVector" << std::endl;
+
+  typedef typename vobj::scalar_type scalar_type;
+  std::vector<scalar_type> lsSum;
+  localSliceInnerProductVector(result, lhs, rhs, lsSum, orthogdim);
+  globalSliceInnerProductVector(result, lhs, lsSum, orthogdim);
+  // std::cout << GridLogMessage << "End mySliceInnerProductVector" << std::endl;
+}
+
+template <class vobj>
+static void localSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, const Lattice<vobj> &rhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  // std::cout << GridLogMessage << "Start prep" << std::endl;
+  typedef typename vobj::vector_type   vector_type;
+  typedef typename vobj::scalar_type   scalar_type;
+  GridBase  *grid = lhs._grid;
+  assert(grid!=NULL);
+  conformable(grid,rhs._grid);
+
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  assert(orthogdim >= 0);
+  assert(orthogdim < Nd);
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+  // std::cout << GridLogMessage << "Start alloc" << std::endl;
+
+  std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd); // will locally sum vectors first
+  lsSum.resize(ld,scalar_type(0.0));                    // sum across these down to scalars
+  std::vector<iScalar<scalar_type>> extracted(Nsimd);   // splitting the SIMD  
+  // std::cout << GridLogMessage << "End alloc" << std::endl;
+
+  result.resize(fd); // And then global sum to return the same vector to every node for IO to file
+  for(int r=0;r<rd;r++){
+    lvSum[r]=zero;
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+  // std::cout << GridLogMessage << "End prep" << std::endl;
+  // std::cout << GridLogMessage << "Start parallel inner product, _rd = " << rd << std::endl;
+  vector_type vv;
+  parallel_for(int r=0;r<rd;r++)
+  {
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+        int ss = so + n * stride + b;
+        vv = TensorRemove(innerProduct(lhs._odata[ss], rhs._odata[ss]));
+        lvSum[r] = lvSum[r] + vv;
+      }
+    }
+  }
+  // std::cout << GridLogMessage << "End parallel inner product" << std::endl;
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  std::vector<int> icoor(Nd);
+  for(int rt=0;rt<rd;rt++){
+
+    iScalar<vector_type> temp; 
+    temp._internal = lvSum[rt];
+    extract(temp,extracted);
+
+    for(int idx=0;idx<Nsimd;idx++){
+
+      grid->iCoorFromIindex(icoor,idx);
+
+      int ldx =rt+icoor[orthogdim]*rd;
+
+      lsSum[ldx]=lsSum[ldx]+extracted[idx]._internal;
+
+    }
+  }
+  // std::cout << GridLogMessage << "End sum over simd lanes" << std::endl;
+}
+template <class vobj>
+static void globalSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  typedef typename vobj::scalar_type scalar_type;
+  GridBase *grid = lhs._grid;
+  int fd = result.size();
+  int ld = lsSum.size();
+  // sum over nodes.
+  std::vector<scalar_type> gsum;
+  gsum.resize(fd, scalar_type(0.0));
+  // std::cout << GridLogMessage << "Start of gsum[t] creation:" << std::endl;
+  for(int t=0;t<fd;t++){
+    int pt = t/ld; // processor plane
+    int lt = t%ld;
+    if ( pt == grid->_processor_coor[orthogdim] ) {
+      gsum[t]=lsSum[lt];
+    }
+  }
+  // std::cout << GridLogMessage << "End of gsum[t] creation:" << std::endl;
+  // std::cout << GridLogMessage << "Start of GlobalSumVector:" << std::endl;
+  grid->GlobalSumVector(&gsum[0], fd);
+  // std::cout << GridLogMessage << "End of GlobalSumVector:" << std::endl;
+
+  result = gsum;
+}
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {

Grid/lattice/Lattice_rng.h

@@ -158,10 +158,19 @@ namespace Grid {
       // tens of seconds per trajectory so this is clean in all reasonable cases,
       // and margin of safety is orders of magnitude.
       // We could hack Sitmo to skip in the higher order words of state if necessary
+      //
+      // Replace with 2^30 ; avoid problem on large volumes
+      //
       /////////////////////////////////////////////////////////////////////////////////////
       //      uint64_t skip = site+1;  //   Old init Skipped then drew.  Checked compat with faster init
+      const int shift = 30;
+
       uint64_t skip = site;
-      skip = skip<<40;
+
+      skip = skip<<shift;
+
+      assert((skip >> shift)==site); // check for overflow
+
       eng.discard(skip);
       //      std::cout << " Engine  " <<site << " state " <<eng<<std::endl;
     } 
@@ -242,7 +251,7 @@ namespace Grid {
 
       dist[0].reset();
       for(int idx=0;idx<words;idx++){
-  fillScalar(buf[idx],dist[0],_generators[0]);
+	fillScalar(buf[idx],dist[0],_generators[0]);
       }
 
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
@@ -274,7 +283,7 @@ namespace Grid {
       RealF *pointer=(RealF *)&l;
       dist[0].reset();
       for(int i=0;i<2*vComplexF::Nsimd();i++){
-  fillScalar(pointer[i],dist[0],_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
@@ -282,7 +291,7 @@ namespace Grid {
       RealD *pointer=(RealD *)&l;
       dist[0].reset();
       for(int i=0;i<2*vComplexD::Nsimd();i++){
-  fillScalar(pointer[i],dist[0],_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
@@ -290,7 +299,7 @@ namespace Grid {
       RealF *pointer=(RealF *)&l;
       dist[0].reset();
       for(int i=0;i<vRealF::Nsimd();i++){
-  fillScalar(pointer[i],dist[0],_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
@@ -308,6 +317,19 @@ namespace Grid {
       std::seed_seq src(seeds.begin(),seeds.end());
       Seed(src,0);
     }
+
+    void SeedUniqueString(const std::string &s){
+      std::vector<int> seeds;
+      std::stringstream sha;
+      seeds = GridChecksum::sha256_seeds(s);
+      for(int i=0;i<seeds.size();i++) { 
+        sha << std::hex << seeds[i];
+      }
+      std::cout << GridLogMessage << "Intialising serial RNG with unique string '" 
+                << s << "'" << std::endl;
+      std::cout << GridLogMessage << "Seed SHA256: " << sha.str() << std::endl;
+      SeedFixedIntegers(seeds);
+    }
   };
 
   class GridParallelRNG : public GridRNGbase {
@@ -368,6 +390,14 @@ namespace Grid {
       _time_counter += usecond()- inner_time_counter;
     };
 
+    void SeedUniqueString(const std::string &s){
+      std::vector<int> seeds;
+      seeds = GridChecksum::sha256_seeds(s);
+      std::cout << GridLogMessage << "Intialising parallel RNG with unique string '" 
+                << s << "'" << std::endl;
+      std::cout << GridLogMessage << "Seed SHA256: " << GridChecksum::sha256_string(seeds) << std::endl;
+      SeedFixedIntegers(seeds);
+    }
     void SeedFixedIntegers(const std::vector<int> &seeds){
 
       // Everyone generates the same seed_seq based on input seeds

Grid/lattice/Lattice_transfer.h

@@ -485,7 +485,7 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
 
 
 template<class vobj>
-void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
+void ExtractSliceLocal(Lattice<vobj> &lowDim, const Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
 {
   typedef typename vobj::scalar_object sobj;
 
@@ -520,7 +520,7 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slic
 
 
 template<class vobj>
-void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
+void Replicate(const Lattice<vobj> &coarse,Lattice<vobj> & fine)
 {
   typedef typename vobj::scalar_object sobj;
 

Grid/log/Log.h

@@ -86,7 +86,7 @@ protected:
   Colours &Painter;
   int active;
   int timing_mode;
-  int topWidth{-1};
+  int topWidth{-1}, chanWidth{-1};
   static int timestamp;
   std::string name, topName;
   std::string COLOUR;
@@ -126,6 +126,7 @@ public:
     }
   }
   void setTopWidth(const int w) {topWidth = w;}
+  void setChanWidth(const int w) {chanWidth = w;}
 
   friend std::ostream& operator<< (std::ostream& stream, Logger& log){
 
@@ -136,7 +137,12 @@ public:
         stream << std::setw(log.topWidth);
       }
       stream << log.topName << log.background()<< " : ";
-      stream << log.colour() <<  std::left << log.name << log.background() << " : ";
+      stream << log.colour() <<  std::left;
+      if (log.chanWidth > 0)
+      {
+        stream << std::setw(log.chanWidth);
+      }
+      stream << log.name << log.background() << " : ";
       if ( log.timestamp ) {
 	log.StopWatch->Stop();
 	GridTime now = log.StopWatch->Elapsed();

Grid/parallelIO/BinaryIO.h

@@ -358,17 +358,10 @@ PARALLEL_CRITICAL
 
       if ( (control & BINARYIO_LEXICOGRAPHIC) && (nrank > 1) ) {
 #ifdef USE_MPI_IO
-	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< " and offset " << offset << std::endl;
+	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< std::endl;
 	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);    assert(ierr==0);
 	ierr=MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);    assert(ierr==0);
 	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);    assert(ierr==0);
-
-	MPI_Offset os;
-	MPI_File_get_position(fh, &os);
-	MPI_File_get_byte_offset(fh, os, &disp);
-	offset = disp;
-
-
 	MPI_File_close(&fh);
 	MPI_Type_free(&fileArray);
 	MPI_Type_free(&localArray);
@@ -377,22 +370,20 @@ PARALLEL_CRITICAL
 #endif
       } else {
 	std::cout << GridLogMessage <<"IOobject: C++ read I/O " << file << " : "
-                  << iodata.size() * sizeof(fobj) << " bytes and offset " << offset << std::endl;
+                  << iodata.size() * sizeof(fobj) << " bytes" << std::endl;
         std::ifstream fin;
 	fin.open(file, std::ios::binary | std::ios::in);
         if (control & BINARYIO_MASTER_APPEND)
-	  {
-	    fin.seekg(-sizeof(fobj), fin.end);
-	  }
+        {
+          fin.seekg(-sizeof(fobj), fin.end);
+        }
         else
-	  {
-	    fin.seekg(offset + myrank * lsites * sizeof(fobj));
-	  }
+        {
+          fin.seekg(offset + myrank * lsites * sizeof(fobj));
+        }
         fin.read((char *)&iodata[0], iodata.size() * sizeof(fobj));
-        
         assert(fin.fail() == 0);
-        offset = fin.tellg();
-	fin.close();
+        fin.close();
       }
       timer.Stop();
 
@@ -647,11 +638,6 @@ PARALLEL_CRITICAL
     IOobject(w,grid,iodata,file,offset,format,BINARYIO_READ|BINARYIO_LEXICOGRAPHIC,
 	     nersc_csum,scidac_csuma,scidac_csumb);
 
-    std::cout << GridLogMessage << "RNG file nersc_checksum   " << std::hex << nersc_csum << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file scidac_checksuma " << std::hex << scidac_csuma << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file scidac_checksumb " << std::hex << scidac_csumb << std::dec << std::endl;
-
-
     timer.Start();
     parallel_for(uint64_t lidx=0;lidx<lsites;lidx++){
       std::vector<RngStateType> tmp(RngStateCount);
@@ -670,11 +656,6 @@ PARALLEL_CRITICAL
       serial.SetState(tmp,0);
     }
 
-    std::cout << GridLogMessage << "RNG file checksum t " << std::hex << nersc_csum_tmp    << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksuma t " << std::hex << scidac_csuma_tmp << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksumb t " << std::hex << scidac_csumb_tmp << std::dec << std::endl;
-
-
     nersc_csum   = nersc_csum   + nersc_csum_tmp;
     scidac_csuma = scidac_csuma ^ scidac_csuma_tmp;
     scidac_csumb = scidac_csumb ^ scidac_csumb_tmp;
@@ -725,11 +706,6 @@ PARALLEL_CRITICAL
 
     IOobject(w,grid,iodata,file,offset,format,BINARYIO_WRITE|BINARYIO_LEXICOGRAPHIC,
 	     nersc_csum,scidac_csuma,scidac_csumb);
-
-    std::cout << GridLogMessage << "RNG file checksum " << std::hex << nersc_csum    << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksuma " << std::hex << scidac_csuma << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksumb " << std::hex << scidac_csumb << std::dec << std::endl;
-   
     iodata.resize(1);
     {
       std::vector<RngStateType> tmp(RngStateCount);
@@ -739,11 +715,6 @@ PARALLEL_CRITICAL
     IOobject(w,grid,iodata,file,offset,format,BINARYIO_WRITE|BINARYIO_MASTER_APPEND,
 	     nersc_csum_tmp,scidac_csuma_tmp,scidac_csumb_tmp);
 
-    std::cout << GridLogMessage << "RNG file checksum t " << std::hex << nersc_csum_tmp    << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksuma t " << std::hex << scidac_csuma_tmp << std::dec << std::endl;
-    std::cout << GridLogMessage << "RNG file checksumb t " << std::hex << scidac_csumb_tmp << std::dec << std::endl;
- 
-
     nersc_csum   = nersc_csum   + nersc_csum_tmp;
     scidac_csuma = scidac_csuma ^ scidac_csuma_tmp;
     scidac_csumb = scidac_csumb ^ scidac_csumb_tmp;

Grid/parallelIO/IldgIO.h

@@ -182,6 +182,11 @@ class GridLimeReader : public BinaryIO {
    {
      filename= _filename;
      File = fopen(filename.c_str(), "r");
+     if (File == nullptr)
+     {
+       std::cerr << "cannot open file '" << filename << "'" << std::endl;
+       abort();
+     }
      LimeR = limeCreateReader(File);
    }
    /////////////////////////////////////////////
@@ -238,57 +243,14 @@ class GridLimeReader : public BinaryIO {
 	// Verify checksums
 	/////////////////////////////////////////////
 	assert(scidacChecksumVerify(scidacChecksum_,scidac_csuma,scidac_csumb)==1);
-  std::cout << GridLogMessage<< " readLimeLatticeBinaryObject checksums match ! " <<std::endl;
 	return;
       }
     }
   }
-
-  ////////////////////////////////////////////
-  // Read an RNG object and verify checksum
-  ////////////////////////////////////////////
-  void readLimeRNGObject(GridSerialRNG &sRNG, GridParallelRNG &pRNG,std::string record_name)
-  {
-    scidacChecksum scidacChecksum_;
-    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
-
-    while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
-      uint64_t file_bytes =limeReaderBytes(LimeR);
-
-      //      std::cout << GridLogMessage << limeReaderType(LimeR) << " "<< file_bytes <<" bytes "<<std::endl;
-      //      std::cout << GridLogMessage<< " readLimeObject seeking "<<  record_name <<" found record :" <<limeReaderType(LimeR) <<std::endl;
-
-      if ( !strncmp(limeReaderType(LimeR), record_name.c_str(),strlen(record_name.c_str()) )  ) {
-
-        const int RngStateCount = GridSerialRNG::RngStateCount;
-        typedef std::array<typename GridSerialRNG::RngStateType,RngStateCount> RNGstate;
-
-	      uint64_t PayloadSize = sizeof(RNGstate) * (pRNG._grid->_gsites+1);
-
-    	  assert(PayloadSize == file_bytes);// Must match or user error
-        uint64_t offset= ftello(File);
-	    	std::cout << GridLogDebug << " ReadLatticeObject from offset "<<offset << std::endl;
-	      BinaryIO::readRNG(sRNG, pRNG, filename, offset, nersc_csum,scidac_csuma,scidac_csumb);
-
-	      /////////////////////////////////////////////
-	      // Insist checksum is next record
-	      /////////////////////////////////////////////
-	      readLimeObject(scidacChecksum_,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM));
-
-	      /////////////////////////////////////////////
-	      // Verify checksums
-	      /////////////////////////////////////////////
-	      assert(scidacChecksumVerify(scidacChecksum_,scidac_csuma,scidac_csumb)==1);
-        std::cout << GridLogMessage<< " readLimeRNGObject checksums match ! " <<std::endl;
-	      return;
-      }
-    }
-  }
   ////////////////////////////////////////////
   // Read a generic serialisable object
   ////////////////////////////////////////////
-  template<class serialisable_object>
-  void readLimeObject(serialisable_object &object,std::string object_name,std::string record_name)
+  void readLimeObject(std::string &xmlstring,std::string record_name)
   {
     // should this be a do while; can we miss a first record??
     while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
@@ -303,15 +265,23 @@ class GridLimeReader : public BinaryIO {
 	limeReaderReadData((void *)&xmlc[0], &nbytes, LimeR);    
 	//	std::cout << GridLogMessage<< " readLimeObject matches XML " << &xmlc[0] <<std::endl;
 
-  std::string xmlstring(&xmlc[0]);
-	XmlReader RD(xmlstring, true, "");
-	read(RD,object_name,object);
+   xmlstring = std::string(&xmlc[0]);
 	return;
       }
 
     }  
     assert(0);
   }
+
+  template<class serialisable_object>
+  void readLimeObject(serialisable_object &object,std::string object_name,std::string record_name)
+  {
+    std::string xmlstring;
+
+    readLimeObject(xmlstring, record_name);
+	  XmlReader RD(xmlstring, true, "");
+	  read(RD,object_name,object);
+  }
 };
 
 class GridLimeWriter : public BinaryIO 
@@ -362,16 +332,11 @@ class GridLimeWriter : public BinaryIO
   ////////////////////////////////////////////
   // Write a generic serialisable object
   ////////////////////////////////////////////
-  template<class serialisable_object>
-  void writeLimeObject(int MB,int ME,serialisable_object &object,std::string object_name,std::string record_name)
+  void writeLimeObject(int MB,int ME,XmlWriter &writer,std::string object_name,std::string record_name)
   {
     if ( boss_node ) {
-      std::string xmlstring;
-      {
-	XmlWriter WR("","");
-	write(WR,object_name,object);
-	xmlstring = WR.XmlString();
-      }
+      std::string xmlstring = writer.docString();
+
       //    std::cout << "WriteLimeObject" << record_name <<std::endl;
       uint64_t nbytes = xmlstring.size();
       //    std::cout << " xmlstring "<< nbytes<< " " << xmlstring <<std::endl;
@@ -385,6 +350,20 @@ class GridLimeWriter : public BinaryIO
       limeDestroyHeader(h);
     }
   }
+
+  template<class serialisable_object>
+  void writeLimeObject(int MB,int ME,serialisable_object &object,std::string object_name,std::string record_name, const unsigned int scientificPrec = 0)
+  {
+    XmlWriter WR("","");
+
+    if (scientificPrec)
+    {
+      WR.scientificFormat(true);
+      WR.setPrecision(scientificPrec);
+    }
+    write(WR,object_name,object);
+    writeLimeObject(MB, ME, WR, object_name, record_name);
+  }
   ////////////////////////////////////////////////////
   // Write a generic lattice field and csum
   // This routine is Collectively called by all nodes
@@ -471,78 +450,6 @@ class GridLimeWriter : public BinaryIO
       writeLimeObject(0,1,checksum,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM));
     }
   }
-
-
-  ////////////////////////////////////////////////////
-  // Write an rng object and csum
-  // This routine is Collectively called by all nodes
-  // in communicator used by the field._grid
-  ////////////////////////////////////////////////////
-  void writeLimeRNGObject(GridSerialRNG &sRNG, GridParallelRNG &pRNG, std::string record_name)
-  {
-    GridBase *grid = pRNG._grid;
-    assert(boss_node == pRNG._grid->IsBoss() );
-
-    const int RngStateCount = GridSerialRNG::RngStateCount;
-    typedef std::array<typename GridSerialRNG::RngStateType,RngStateCount> RNGstate;
-
-    ////////////////////////////////////////////
-    // Create record header
-    ////////////////////////////////////////////
-    int err;
-    uint32_t nersc_csum,scidac_csuma,scidac_csumb;
-    uint64_t PayloadSize = sizeof(RNGstate) * (grid->_gsites+1);
-    std::cout << GridLogDebug << "Computed payload size " << PayloadSize << std::endl;
-    if ( boss_node ) {
-      createLimeRecordHeader(record_name, 0, 0, PayloadSize);
-      fflush(File);
-    }
-    
-    ////////////////////////////////////////////////
-    // Check all nodes agree on file position
-    ////////////////////////////////////////////////
-    uint64_t offset1;
-    if ( boss_node ) {
-      offset1 = ftello(File);    
-    }
-    grid->Broadcast(0,(void *)&offset1,sizeof(offset1));
-
-    ///////////////////////////////////////////
-    // The above is collective. Write by other means into the binary record
-    ///////////////////////////////////////////
-    uint64_t offset = offset1;
-    BinaryIO::writeRNG(sRNG, pRNG,filename, offset, nersc_csum,scidac_csuma,scidac_csumb);
-
-    ///////////////////////////////////////////
-    // Wind forward and close the record
-    ///////////////////////////////////////////
-    if ( boss_node ) {
-      fseek(File,0,SEEK_END);             
-      uint64_t offset2 = ftello(File);    
-      std::cout << GridLogDebug << " now at offset "<<offset2 << std::endl;
-      assert( (offset2-offset1) == PayloadSize);
-    }
-
-    /////////////////////////////////////////////////////////////
-    // Check MPI-2 I/O did what we expect to file
-    /////////////////////////////////////////////////////////////
-
-    if ( boss_node ) { 
-      err=limeWriterCloseRecord(LimeW);  assert(err>=0);
-    }
-    ////////////////////////////////////////
-    // Write checksum element, propagaing forward from the BinaryIO
-    // Always pair a checksum with a binary object, and close message
-    ////////////////////////////////////////
-    scidacChecksum checksum;
-    std::stringstream streama; streama << std::hex << scidac_csuma;
-    std::stringstream streamb; streamb << std::hex << scidac_csumb;
-    checksum.suma= streama.str();
-    checksum.sumb= streamb.str();
-    if ( boss_node ) { 
-      writeLimeObject(0,1,checksum,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM));
-    }
-  }
 };
 
 class ScidacWriter : public GridLimeWriter {
@@ -559,32 +466,12 @@ class ScidacWriter : public GridLimeWriter {
       writeLimeObject(0,1,_userFile,_userFile.SerialisableClassName(),std::string(SCIDAC_FILE_XML));
     }
   }
-
-  void writeScidacRNGRecord(GridSerialRNG &sRNG, GridParallelRNG &pRNG) 
-  {
-    GridBase *grid = pRNG._grid;
-    FieldMetaData header;
-
-    header.floating_point = "IEEE64BIG";
-    header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
-    GridMetaData(grid,header); 
-    MachineCharacteristics(header);
-    
-    //////////////////////////////////////////////
-    // Fill the Lime file record by record
-    //////////////////////////////////////////////
-    if ( this->boss_node ) {
-      writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
-    }
-    // Collective call
-    writeLimeRNGObject(sRNG, pRNG,std::string(ILDG_BINARY_DATA));      // Closes message with checksum
-  }
-  
   ////////////////////////////////////////////////
   // Write generic lattice field in scidac format
   ////////////////////////////////////////////////
   template <class vobj, class userRecord>
-  void writeScidacFieldRecord(Lattice<vobj> &field,userRecord _userRecord) 
+  void writeScidacFieldRecord(Lattice<vobj> &field,userRecord _userRecord,
+                              const unsigned int recordScientificPrec = 0) 
   {
     GridBase * grid = field._grid;
 
@@ -602,13 +489,12 @@ class ScidacWriter : public GridLimeWriter {
     //////////////////////////////////////////////
     if ( this->boss_node ) {
       writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
-      writeLimeObject(0,0,_userRecord,_userRecord.SerialisableClassName(),std::string(SCIDAC_RECORD_XML));
+      writeLimeObject(0,0,_userRecord,_userRecord.SerialisableClassName(),std::string(SCIDAC_RECORD_XML), recordScientificPrec);
       writeLimeObject(0,0,_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML));
     }
     // Collective call
     writeLimeLatticeBinaryObject(field,std::string(ILDG_BINARY_DATA));      // Closes message with checksum
   }
-  
 };
 
 
@@ -622,27 +508,8 @@ class ScidacReader : public GridLimeReader {
      readLimeObject(_scidacFile,_scidacFile.SerialisableClassName(),std::string(SCIDAC_PRIVATE_FILE_XML));
      readLimeObject(_userFile,_userFile.SerialisableClassName(),std::string(SCIDAC_FILE_XML));
    }
-
   ////////////////////////////////////////////////
-  // Read RNGobject in scidac format
-  ////////////////////////////////////////////////
-  void readScidacRNGRecord(GridSerialRNG &sRNG, GridParallelRNG &pRNG) 
-  {
-    GridBase * grid = pRNG._grid;
-
-    ////////////////////////////////////////
-    // fill the Grid header
-    ////////////////////////////////////////
-    FieldMetaData header;
- 
-    //////////////////////////////////////////////
-    // Fill the Lime file record by record
-    //////////////////////////////////////////////
-    readLimeObject(header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
-    readLimeRNGObject(sRNG, pRNG, std::string(ILDG_BINARY_DATA));
-  }
-  ////////////////////////////////////////////////
-  // Read generic lattice field in scidac format
+  // Write generic lattice field in scidac format
   ////////////////////////////////////////////////
   template <class vobj, class userRecord>
   void readScidacFieldRecord(Lattice<vobj> &field,userRecord &_userRecord) 

Grid/perfmon/Timer.h

@@ -49,7 +49,8 @@ inline double usecond(void) {
 
 typedef  std::chrono::system_clock          GridClock;
 typedef  std::chrono::time_point<GridClock> GridTimePoint;
-typedef  std::chrono::milliseconds          GridTime;
+typedef  std::chrono::milliseconds          GridMillisecs;
+typedef  std::chrono::microseconds          GridTime;
 typedef  std::chrono::microseconds          GridUsecs;
 
 inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milliseconds & time)
@@ -57,6 +58,11 @@ inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milli
   stream << time.count()<<" ms";
   return stream;
 }
+inline std::ostream& operator<< (std::ostream & stream, const std::chrono::microseconds & time)
+{
+  stream << time.count()<<" usec";
+  return stream;
+}
  
 class GridStopWatch {
 private:
@@ -96,6 +102,9 @@ public:
     assert(running == false);
     return (uint64_t) accumulator.count();
   }
+  bool isRunning(void){
+    return running;
+  }
 };
 
 }

Grid/qcd/QCD.h

@@ -90,17 +90,20 @@ namespace QCD {
     // That probably makes for GridRedBlack4dCartesian grid.
 
     // s,sp,c,spc,lc
-    template<typename vtype> using iSinglet                   = iScalar<iScalar<iScalar<vtype> > >;
-    template<typename vtype> using iSpinMatrix                = iScalar<iMatrix<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourMatrix              = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
-    template<typename vtype> using iSpinColourMatrix          = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
-    template<typename vtype> using iLorentzColourMatrix       = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
-    template<typename vtype> using iDoubleStoredColourMatrix  = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
-    template<typename vtype> using iSpinVector                = iScalar<iVector<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourVector              = iScalar<iScalar<iVector<vtype, Nc> > >;
-    template<typename vtype> using iSpinColourVector          = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
-    template<typename vtype> using iHalfSpinVector            = iScalar<iVector<iScalar<vtype>, Nhs> >;
-    template<typename vtype> using iHalfSpinColourVector      = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+
+    template<typename vtype> using iSinglet                     = iScalar<iScalar<iScalar<vtype> > >;
+    template<typename vtype> using iSpinMatrix                  = iScalar<iMatrix<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourMatrix                = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
+    template<typename vtype> using iSpinColourMatrix            = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
+    template<typename vtype> using iLorentzColourMatrix         = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
+    template<typename vtype> using iDoubleStoredColourMatrix    = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
+    template<typename vtype> using iSpinVector                  = iScalar<iVector<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourVector                = iScalar<iScalar<iVector<vtype, Nc> > >;
+    template<typename vtype> using iSpinColourVector            = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
+    template<typename vtype> using iHalfSpinVector              = iScalar<iVector<iScalar<vtype>, Nhs> >;
+    template<typename vtype> using iHalfSpinColourVector        = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+    template<typename vtype> using iSpinColourSpinColourMatrix  = iScalar<iMatrix<iMatrix<iMatrix<iMatrix<vtype, Nc>, Ns>, Nc>, Ns> >;
+
 
     template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Ns>, Ngp >;
     template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
@@ -127,10 +130,28 @@ namespace QCD {
     typedef iSpinColourMatrix<Complex  >    SpinColourMatrix;
     typedef iSpinColourMatrix<ComplexF >    SpinColourMatrixF;
     typedef iSpinColourMatrix<ComplexD >    SpinColourMatrixD;
-
+    
     typedef iSpinColourMatrix<vComplex >    vSpinColourMatrix;
     typedef iSpinColourMatrix<vComplexF>    vSpinColourMatrixF;
     typedef iSpinColourMatrix<vComplexD>    vSpinColourMatrixD;
+    
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
+
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
 
     // LorentzColour
     typedef iLorentzColourMatrix<Complex  > LorentzColourMatrix;
@@ -229,6 +250,9 @@ namespace QCD {
     typedef Lattice<vSpinColourMatrixF>     LatticeSpinColourMatrixF;
     typedef Lattice<vSpinColourMatrixD>     LatticeSpinColourMatrixD;
 
+    typedef Lattice<vSpinColourSpinColourMatrix>      LatticeSpinColourSpinColourMatrix;
+    typedef Lattice<vSpinColourSpinColourMatrixF>     LatticeSpinColourSpinColourMatrixF;
+    typedef Lattice<vSpinColourSpinColourMatrixD>     LatticeSpinColourSpinColourMatrixD;
 
     typedef Lattice<vLorentzColourMatrix>  LatticeLorentzColourMatrix;
     typedef Lattice<vLorentzColourMatrixF> LatticeLorentzColourMatrixF;