Merge branch 'develop' of https://github.com/paboyle/Grid into feature/staggering

tests/Test_stencil.cc

@@ -48,7 +48,7 @@ int main(int argc, char ** argv) {
   double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
 
   GridCartesian Fine(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian rbFine(&Fine);
   GridParallelRNG       fRNG(&Fine);
 
   //  fRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

tests/core/Test_cshift_red_black.cc

@@ -47,7 +47,7 @@ int main (int argc, char ** argv)
   mask[0]=0;
 
   GridCartesian         Fine  (latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1);
+  GridRedBlackCartesian RBFine(&Fine,mask,1);
 
   GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
 

tests/core/Test_cshift_red_black_rotate.cc

@@ -47,7 +47,7 @@ int main (int argc, char ** argv)
   mask[0]=0;
 
   GridCartesian         Fine  (latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1);
+  GridRedBlackCartesian RBFine(&Fine,mask,1);
 
   GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
 

tests/core/Test_fft.cc

@@ -47,7 +47,7 @@ int main (int argc, char ** argv)
     vol = vol * latt_size[d];
   }
   GridCartesian         GRID(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian RBGRID(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian RBGRID(&GRID);
 
   LatticeComplexD     one(&GRID);
   LatticeComplexD      zz(&GRID);

tests/core/Test_gpwilson_even_odd.cc

@@ -40,7 +40,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/core/Test_main.cc

@@ -84,7 +84,7 @@ int main(int argc, char **argv) {
       double volume = latt_size[0] * latt_size[1] * latt_size[2] * latt_size[3];
 
       GridCartesian Fine(latt_size, simd_layout, mpi_layout);
-      GridRedBlackCartesian rbFine(latt_size, simd_layout, mpi_layout);
+      GridRedBlackCartesian rbFine(&Fine);
       GridParallelRNG FineRNG(&Fine);
       GridSerialRNG SerialRNG;
       GridSerialRNG SerialRNG1;

tests/core/Test_staggered.cc

@@ -40,7 +40,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/core/Test_wilson_even_odd.cc

@@ -51,7 +51,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/core/Test_wilson_twisted_mass_even_odd.cc

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/forces/Test_gp_rect_force.cc

@@ -42,7 +42,7 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/forces/Test_laplacian_force.cc

@@ -42,7 +42,7 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/forces/Test_rect_force.cc

@@ -42,7 +42,7 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/forces/Test_wilson_force.cc

@@ -42,7 +42,7 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
 
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   int threads = GridThread::GetThreads();
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;

tests/smearing/Test_WilsonFlow.cc

@@ -71,7 +71,7 @@ int main(int argc, char **argv) {
   std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1, 2, 3, 4, 5});
   GridSerialRNG sRNG;
@@ -149,4 +149,4 @@ JSON
 }
 
 
-*/
+*/

tests/solver/Test_dwf_mrhs_cg.cc

@@ -0,0 +1,195 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_mrhs_cg.cc
+
+    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 */
+#include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  typedef typename DomainWallFermionR::FermionField FermionField; 
+  typedef typename DomainWallFermionR::ComplexField ComplexField; 
+  typename DomainWallFermionR::ImplParams params; 
+
+  const int Ls=8;
+
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  std::vector<int> mpi_split (mpi_layout.size(),1);
+
+  std::cout << "UGrid (world root)"<<std::endl;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+
+  std::cout << "FGrid (child of UGrid)"<<std::endl;
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+
+  int nrhs = UGrid->RankCount() ;
+
+  /////////////////////////////////////////////
+  // Split into 1^4 mpi communicators
+  /////////////////////////////////////////////
+  std::cout << "SGrid (world root)"<<std::endl;
+  GridCartesian         * SGrid = new GridCartesian(GridDefaultLatt(),
+						    GridDefaultSimd(Nd,vComplex::Nsimd()),
+						    mpi_split,
+						    *UGrid); 
+
+  GridCartesian         * SFGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,SGrid);
+  std::cout << "SFGrid"<<std::endl;
+  GridRedBlackCartesian * SrbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(SGrid);
+  std::cout << "SrbGrid"<<std::endl;
+  GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,SGrid);
+  std::cout << "SFrbGrid"<<std::endl;
+
+  ///////////////////////////////////////////////
+  // Set up the problem as a 4d spreadout job
+  ///////////////////////////////////////////////
+  std::vector<int> seeds({1,2,3,4});
+
+  GridParallelRNG pRNG(UGrid );  pRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
+  std::vector<FermionField>    src(nrhs,FGrid);
+  std::vector<FermionField> result(nrhs,FGrid);
+
+  for(int s=0;s<nrhs;s++) random(pRNG5,src[s]);
+  for(int s=0;s<nrhs;s++) result[s] = zero;
+
+  LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu);
+
+  ///////////////////////////////////////////////////////////////
+  // Bounce these fields to disk
+  ///////////////////////////////////////////////////////////////
+
+  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+  std::cout << GridLogMessage << " Writing out in parallel view "<<std::endl;
+  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+  emptyUserRecord record;
+  std::string file("./scratch.scidac");
+  std::string filef("./scratch.scidac.ferm");
+  int me = UGrid->ThisRank();
+  LatticeGaugeField s_Umu(SGrid);
+  FermionField s_src(SFGrid);
+  FermionField s_res(SFGrid);
+
+  {
+    FGrid->Barrier();
+    ScidacWriter _ScidacWriter;
+    _ScidacWriter.open(file);
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    std::cout << GridLogMessage << " Writing out gauge field "<<std::endl;
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    _ScidacWriter.writeScidacFieldRecord(Umu,record);
+    _ScidacWriter.close();
+    FGrid->Barrier();
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    std::cout << GridLogMessage << " Reading in gauge field "<<std::endl;
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    ScidacReader  _ScidacReader;
+    _ScidacReader.open(file);
+    _ScidacReader.readScidacFieldRecord(s_Umu,record);
+    _ScidacReader.close();
+    FGrid->Barrier();
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    std::cout << GridLogMessage << " Read in gauge field "<<std::endl;
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+  }
+
+
+  {
+    for(int n=0;n<nrhs;n++){
+
+      std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+      std::cout << GridLogMessage << " Writing out record "<<n<<std::endl;
+      std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+
+      std::stringstream filefn;      filefn << filef << "."<< n;
+      ScidacWriter _ScidacWriter;
+      _ScidacWriter.open(filefn.str());
+      _ScidacWriter.writeScidacFieldRecord(src[n],record);
+      _ScidacWriter.close();
+    }
+      
+    FGrid->Barrier();
+
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+    std::cout << GridLogMessage << " Reading back in the single process view "<<std::endl;
+    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+      
+    for(int n=0;n<nrhs;n++){
+      if ( n==me ) { 
+	std::stringstream filefn;	filefn << filef << "."<< n;
+	ScidacReader  _ScidacReader;
+	_ScidacReader.open(filefn.str());
+	_ScidacReader.readScidacFieldRecord(s_src,record);
+	_ScidacReader.close();
+      }
+    }
+    FGrid->Barrier();
+  }
+
+
+  ///////////////////////////////////////////////////////////////
+  // Set up N-solvers as trivially parallel
+  ///////////////////////////////////////////////////////////////
+
+  RealD mass=0.01;
+  RealD M5=1.8;
+  DomainWallFermionR Ddwf(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5);
+
+  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+  std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
+  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+
+  MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
+  ConjugateGradient<FermionField> CG((1.0e-8/(me+1)),10000);
+  s_res = zero;
+  CG(HermOp,s_src,s_res);
+
+  ///////////////////////////////////////
+  // Share the information
+  ///////////////////////////////////////
+  std::vector<uint32_t> iterations(nrhs,0);
+  iterations[me] = CG.IterationsToComplete;
+
+  for(int n=0;n<nrhs;n++){
+    UGrid->GlobalSum(iterations[n]);
+  }
+
+  /////////////////////////////////////////////////////////////
+  // Report how long they all took
+  /////////////////////////////////////////////////////////////
+  for(int r=0;r<nrhs;r++){
+    std::cout << GridLogMessage<<" Rank "<<r<<" "<< iterations[r]<<" CG iterations"<<std::endl;
+  }
+  Grid_finalize();
+}

tests/solver/Test_laplacian.cc

@@ -40,7 +40,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4,5});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_staggered_block_cg_unprec.cc

@@ -27,7 +27,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
-#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 
 using namespace std;
 using namespace Grid;

tests/solver/Test_staggered_cg_prec.cc

@@ -57,7 +57,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_staggered_cg_unprec.cc

@@ -57,7 +57,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_wilson_cg_prec.cc

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_wilson_cg_schur.cc

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_wilson_cg_unprec.cc

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);

tests/solver/Test_wilson_cr_unprec.cc

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
   GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(&Grid);
 
   std::vector<int> seeds({1,2,3,4});
   GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);