Reorganise a little to let the PV inverter be defined outside

the Reconstruct class.

This lets the multiple choices for PV inversion be composed without
changing the routine and no if/else case enumeration.

Implemented SchurDiagMooee PV inversion (red black) and Unprec PV inversion.
Red black cuts from 190 iterations to 90 iterations at 10^-12 on 8^4 test system

Will revisit multiple Schur options and add a Fourier based multishift PV inverse, similar
to the one Rudy Arthur did in BFM

.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/*
+Grid/Eigen
-
+Eigen/*
-# FFTW source #
-################
-lib/fftw/*
 
 # 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
+Grid/qcd/spin/gamma-gen/*.h
-lib/qcd/spin/gamma-gen/*.cc
+Grid/qcd/spin/gamma-gen/*.cc
-lib/version.h
-
-# vs code editor files #
-########################
-.vscode/
-.vscode/settings.json
-settings.json

.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,7 +21,7 @@ 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`
@@ -33,6 +38,7 @@ 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
@@ -49,12 +55,7 @@ script:
     - make -j4
     - make install
     - cd $CWD/build
-    - ../configure --enable-precision=single --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install
+    - ../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 --with-lime=$CWD/build/lime/install
-    - 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,10 +1,14 @@
 #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"
 #endif
-#include <Eigen/Dense>
+#include <Grid/Eigen/Dense>
 #if defined __GNUC__
 #pragma GCC diagnostic pop
 #endif

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)
+libGrid_adir                   = $(includedir)/Grid
-nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) $(eigen_unsupp_files) Config.h
+nobase_dist_pkginclude_HEADERS = $(HFILES) $(eigen_files) $(eigen_unsupp_files)

Grid/algorithms/Algorithms.h

@@ -52,6 +52,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/CoarsenedMatrix.h>
 #include <Grid/algorithms/FFT.h>
 
+
 // EigCg
 // Pcg
 // Hdcg

Grid/algorithms/iterative/Deflation.h

@@ -31,21 +31,13 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 namespace Grid { 
 
 template<class Field>
-class Guesser {
+class ZeroGuesser: public LinearFunction<Field> {
-public:
-  Guesser(void) = default;
-  virtual ~Guesser(void) = default;
-  virtual void operator()(const Field &src, Field &guess) = 0;
-};
-
-template<class Field>
-class ZeroGuesser: public Guesser<Field> {
 public:
   virtual void operator()(const Field &src, Field &guess) { guess = zero; };
 };
 
 template<class Field>
-class SourceGuesser: public Guesser<Field> {
+class SourceGuesser: public LinearFunction<Field> {
 public:
   virtual void operator()(const Field &src, Field &guess) { guess = src; };
 };
@@ -54,7 +46,7 @@ public:
 // Fine grid deflation
 ////////////////////////////////
 template<class Field>
-class DeflatedGuesser: public Guesser<Field> {
+class DeflatedGuesser: public LinearFunction<Field> {
 private:
   const std::vector<Field> &evec;
   const std::vector<RealD> &eval;
@@ -76,7 +68,7 @@ public:
 };
 
 template<class FineField, class CoarseField>
-class LocalCoherenceDeflatedGuesser: public Guesser<FineField> {
+class LocalCoherenceDeflatedGuesser: public LinearFunction<FineField> {
 private:
   const std::vector<FineField>   &subspace;
   const std::vector<CoarseField> &evec_coarse;

Grid/algorithms/iterative/Reconstruct5Dprop.h

@@ -0,0 +1,186 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/SchurRedBlack.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+
+template<class Field>
+class PauliVillarsSolverUnprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverUnprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    MdagMLinearOperator<Matrix,Field> HermOp(_Matrix);
+
+    _Matrix.SetMass(1.0);
+    _Matrix.Mdag(src,A);
+    CG(HermOp,A,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field>
+class PauliVillarsSolverRBprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverRBprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    _Matrix.SetMass(1.0);
+    SchurRedBlackDiagMooeeSolve<Field> SchurSolver(CG);
+    SchurSolver(_Matrix,src,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class PVinverter> class Reconstruct5DfromPhysical {
+ private:
+  PVinverter & PauliVillarsSolver;
+ public:
+
+ /////////////////////////////////////////////////////
+ // First cut works, 10 Oct 2018.
+ //
+ // Must form a plan to get this into production for Zmobius acceleration
+ // of the Mobius exact AMA corrections.
+ //
+ // TODO : understand absence of contact term in eqns in Hantao's thesis
+ //        sol4 is contact term subtracted.
+ //
+ // Options
+ // a) Defect correction approach:
+ //    1) Compute defect from current soln (initially guess).
+ //       This is ...... outerToInner check !!!!
+ //    2) Deflated Zmobius solve to get 4d soln
+ //       Ensure deflation is working
+ //    3) Refine 5d Outer using the inner 4d delta soln
+ // 
+ // Step 1: localise PV inverse in a routine. [DONE]
+ // Step 2: Schur based PV inverse            [DONE]
+ // Step 3: Fourier accelerated PV inverse
+ // Step 4: 
+ /////////////////////////////////////////////////////
+ 
+  Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) 
+    : PauliVillarsSolver(_PauliVillarsSolver) 
+  { 
+  };
+
+
+   template<class Matrix>
+   void PV(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.M(src,sol);
+     _Matrix.SetMass(m);
+   }
+   template<class Matrix>
+   void PVdag(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.Mdag(src,sol);
+     _Matrix.SetMass(m);
+   }
+  template<class Matrix>
+  void operator() (Matrix & _Matrix,const Field &sol4,const Field &src4, Field &sol5){
+
+    int Ls =  _Matrix.Ls;
+
+    Field psi4(_Matrix.GaugeGrid());
+    Field psi(_Matrix.FermionGrid());
+    Field A  (_Matrix.FermionGrid());
+    Field B  (_Matrix.FermionGrid());
+    Field c  (_Matrix.FermionGrid());
+
+    typedef typename Matrix::Coeff_t Coeff_t;
+
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< " Reconstruct5Dprop: c.f. MADWF algorithm         " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    _Matrix.ImportPhysicalFermionSource(src4,B); 
+
+    ///////////////////////////////////////
+    // Set up c from src4
+    ///////////////////////////////////////
+    PauliVillarsSolver(_Matrix,B,A);
+    _Matrix.Pdag(A,c);
+
+    //////////////////////////////////////
+    // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+    //////////////////////////////////////
+    psi4 = - sol4;
+    InsertSlice(psi4, psi, 0   , 0);
+    for (int s=1;s<Ls;s++) {
+      ExtractSlice(psi4,c,s,0);
+       InsertSlice(psi4,psi,s,0);
+    }
+
+    /////////////////////////////
+    // Pdag PV^-1 Dm P 
+    /////////////////////////////
+    _Matrix.P(psi,B);
+    _Matrix.M(B,A);
+    PauliVillarsSolver(_Matrix,A,B);
+    _Matrix.Pdag(B,A);
+
+    //////////////////////////////
+    // Reinsert surface prop
+    //////////////////////////////
+    InsertSlice(sol4,A,0,0);
+
+    //////////////////////////////
+    // Convert from y back to x 
+    //////////////////////////////
+    _Matrix.P(A,sol5);
+    
+  }
+};
+
+}
+}

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_rng.h

@@ -251,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));
@@ -283,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));
     }
@@ -291,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));
     }
@@ -299,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));
     }
@@ -317,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 {
@@ -377,6 +390,18 @@ namespace Grid {
       _time_counter += usecond()- inner_time_counter;
     };
 
+    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 parallel RNG with unique string '" 
+                << s << "'" << std::endl;
+      std::cout << GridLogMessage << "Seed SHA256: " << sha.str() << std::endl;
+      SeedFixedIntegers(seeds);
+    }
     void SeedFixedIntegers(const std::vector<int> &seeds){
 
       // Everyone generates the same seed_seq based on input seeds

Grid/parallelIO/IldgIO.h

@@ -250,8 +250,7 @@ class GridLimeReader : public BinaryIO {
   ////////////////////////////////////////////
   // Read a generic serialisable object
   ////////////////////////////////////////////
-  template<class serialisable_object>
+  void readLimeObject(std::string &xmlstring,std::string record_name)
-  void readLimeObject(serialisable_object &object,std::string object_name,std::string record_name)
   {
     // should this be a do while; can we miss a first record??
     while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
@@ -266,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]);
+   xmlstring = std::string(&xmlc[0]);
-	XmlReader RD(xmlstring, true, "");
-	read(RD,object_name,object);
 	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 

Grid/perfmon/Timer.h

@@ -102,6 +102,9 @@ public:
     assert(running == false);
     return (uint64_t) accumulator.count();
   }
+  bool isRunning(void){
+    return running;
+  }
 };
 
 }