Working aas good as possible on 48^3 in double

Grid/algorithms/iterative/ConjugateGradient.h

@@ -54,11 +54,14 @@ public:
   ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true)
     : Tolerance(tol),
       MaxIterations(maxit),
-      ErrorOnNoConverge(err_on_no_conv){};
+      ErrorOnNoConverge(err_on_no_conv)
+  {};
 
   void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
 
     GRID_TRACE("ConjugateGradient");
+    GridStopWatch PreambleTimer;
+    PreambleTimer.Start();
     psi.Checkerboard() = src.Checkerboard();
 
     conformable(psi, src);
@@ -66,22 +69,26 @@ public:
     RealD cp, c, a, d, b, ssq, qq;
     //RealD b_pred;
 
-    Field p(src);
-    Field mmp(src);
-    Field r(src);
+    // Was doing copies
+    Field p(src.Grid());
+    Field mmp(src.Grid());
+    Field r(src.Grid());
 
     // Initial residual computation & set up
+    ssq = norm2(src);
     RealD guess = norm2(psi);
     assert(std::isnan(guess) == 0);
-    
-    Linop.HermOpAndNorm(psi, mmp, d, b);
-    
-    r = src - mmp;
-    p = r;
-
-    a = norm2(p);
+    if ( guess == 0.0 ) {
+      r = src;
+      p = r;
+      a = ssq;
+    } else { 
+      Linop.HermOpAndNorm(psi, mmp, d, b);
+      r = src - mmp;
+      p = r;
+      a = norm2(p);
+    }
     cp = a;
-    ssq = norm2(src);
 
     // Handle trivial case of zero src
     if (ssq == 0.){
@@ -103,7 +110,7 @@ public:
     // Check if guess is really REALLY good :)
     if (cp <= rsq) {
       TrueResidual = std::sqrt(a/ssq);
-      std::cout << GridLogMessage << "ConjugateGradient guess is converged already : cp " << cp <<" rsq "<<rsq <<" ssq "<<ssq<< std::endl;
+      std::cout << GridLogMessage << "ConjugateGradient guess is converged already " << std::endl;
       IterationsToComplete = 0;	
       return;
     }
@@ -111,6 +118,7 @@ public:
     std::cout << GridLogIterative << std::setprecision(8)
               << "ConjugateGradient: k=0 residual " << cp << " target " << rsq << std::endl;
 
+    PreambleTimer.Stop();
     GridStopWatch LinalgTimer;
     GridStopWatch InnerTimer;
     GridStopWatch AxpyNormTimer;
@@ -183,7 +191,8 @@ public:
 		  << "\tTrue residual " << true_residual
 		  << "\tTarget " << Tolerance << std::endl;
 
-	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	//	std::cout << GridLogMessage << "\tPreamble   " << PreambleTimer.Elapsed() <<std::endl;
+	std::cout << GridLogMessage << "\tSolver Elapsed    " << SolverTimer.Elapsed() <<std::endl;
         std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
 	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
 	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
@@ -202,13 +211,22 @@ public:
       }
     }
     // Failed. Calculate true residual before giving up                                                         
-    Linop.HermOpAndNorm(psi, mmp, d, qq);
-    p = mmp - src;
-
-    TrueResidual = sqrt(norm2(p)/ssq);
+    // Linop.HermOpAndNorm(psi, mmp, d, qq);
+    //    p = mmp - src;
+    //TrueResidual = sqrt(norm2(p)/ssq);
+    //    TrueResidual = 1;
 
     std::cout << GridLogMessage << "ConjugateGradient did NOT converge "<<k<<" / "<< MaxIterations
-	      <<" residual "<< TrueResidual<< std::endl;
+    	      <<" residual "<< std::sqrt(cp / ssq)<< std::endl;
+    SolverTimer.Stop();
+    std::cout << GridLogMessage << "\tPreamble   " << PreambleTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tSolver     " << SolverTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "Solver breakdown "<<std::endl;
+    std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage<< "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+    std::cout << GridLogPerformance << "\t\tInner      " << InnerTimer.Elapsed() <<std::endl;
+    std::cout << GridLogPerformance << "\t\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+    std::cout << GridLogPerformance << "\t\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
     if (ErrorOnNoConverge) assert(0);
     IterationsToComplete = k;

Grid/util/FlightRecorder.cc

@@ -285,12 +285,6 @@ void FlightRecorder::xmitLog(void *buf,uint64_t bytes)
     XmitLoggingCounter++;
   }
 #endif
-  } else {
-    uint64_t word = 1;
-    deviceVector<uint64_t> dev(1);
-    acceleratorCopyToDevice(&word,&dev[0],sizeof(uint64_t));
-    acceleratorCopySynchronise();
-    MPI_Barrier(MPI_COMM_WORLD);
   }
 }
 void FlightRecorder::recvLog(void *buf,uint64_t bytes,int rank)

HMC/site_autocorrelation.cc

@@ -0,0 +1,90 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: 
+
+Copyright (C) 2017
+
+Author: Peter Boyle
+
+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 <string>
+
+template <class T> void readFile(T& out, std::string const fname){
+  Grid::emptyUserRecord record;
+  Grid::ScidacReader RD;
+  RD.open(fname);
+  RD.readScidacFieldRecord(out,record);
+  RD.close();
+}
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
+
+  LatticeComplexD plaq1(&Grid), plaq2(&Grid);
+
+  FieldMetaData header;
+
+  double vol = plaq1.Grid()->gSites();
+  
+  std::string file1(argv[1]);
+  std::cout << "Reading "<<file1<<std::endl;
+  readFile(plaq1,file1);
+  std::string file2(argv[2]);
+  std::cout << "Reading "<<file2<<std::endl;
+  readFile(plaq2,file2);
+  
+  auto p1bar = TensorRemove(sum(plaq1));
+  auto p2bar = TensorRemove(sum(plaq2));
+
+  p1bar = p1bar / vol;
+  p2bar = p2bar / vol;
+
+  std::cout<< GridLogMessage << "p1bar = "<<p1bar<<std::endl;
+  std::cout<< GridLogMessage << "p2bar = "<<p2bar<<std::endl;
+
+  auto corr_site = plaq1 * plaq2 - p1bar * p2bar;
+  auto corr_bar  = TensorRemove(sum(corr_site))/vol;
+
+  auto cov1_site = plaq1 * plaq1 - p1bar * p1bar;
+  auto cov1_bar  = TensorRemove(sum(cov1_site))/vol;
+
+  auto cov2_site = plaq2 * plaq2 - p2bar * p2bar;
+  auto cov2_bar  = TensorRemove(sum(cov2_site))/vol;
+
+  std::cout<< GridLogMessage << "cov_bar = "<<corr_bar<<std::endl;
+
+  std::cout<< GridLogMessage << "corr_bar = "<<corr_bar/sqrt(cov1_bar*cov2_bar)<<std::endl;
+  
+  Grid_finalize();
+}  // main
+
+

HMC/site_plaquette.cc

@@ -0,0 +1,79 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: 
+
+Copyright (C) 2017
+
+Author: Peter Boyle
+
+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 <string>
+
+NAMESPACE_BEGIN(Grid);
+template <class T> void writeFile(T& out, std::string const fname){
+  emptyUserRecord record;
+  ScidacWriter WR(out.Grid()->IsBoss());
+  WR.open(fname);
+  WR.writeScidacFieldRecord(out,record,0,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+  WR.close();
+}
+NAMESPACE_END(Grid);
+int main(int argc, char **argv) {
+  using namespace Grid;
+  
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size, simd_layout, mpi_layout);
+
+  LatticeGaugeField Umu(&Grid);
+  std::vector<LatticeColourMatrix> U(4,&Grid);
+  LatticeComplexD plaq(&Grid);
+
+  FieldMetaData header;
+
+  double vol = Umu.Grid()->gSites();
+  double faces = (1.0 * Nd * (Nd - 1)) / 2.0;
+  double Ncdiv = 1.0/Nc;
+  
+  std::string file1(argv[1]);
+  std::string file2(argv[2]);
+  std::cout << "Reading "<<file1<<std::endl;
+  NerscIO::readConfiguration(Umu,header,file1);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+  SU3WilsonLoops::sitePlaquette(plaq,U);
+
+  plaq = plaq *(Ncdiv/faces);
+  
+  std::cout << "Writing "<<file2<<std::endl;
+  writeFile(plaq,file2);
+  
+  Grid_finalize();
+}  // main
+
+

tests/debug/Test_general_coarse_hdcg_phys48.cc

@@ -200,7 +200,7 @@ int main (int argc, char ** argv)
   Grid_init(&argc,&argv);
 
   const int Ls=24;
-  const int nbasis = 60;
+  const int nbasis = 62;
   const int cb = 0 ;
   RealD mass=0.00078;
   RealD M5=1.8;
@@ -274,7 +274,7 @@ int main (int argc, char ** argv)
   std::string ldop_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/LittleDiracOp.phys48.new.62");
   std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/evecs.scidac");
   std::string eval_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/eval.xml");
-  bool load_agg=true;
+  bool load_agg=false;
   bool load_refine=false;
   bool load_mat=false;
   bool load_evec=false;
@@ -287,8 +287,8 @@ int main (int argc, char ** argv)
   } else {
     //    Aggregates.CreateSubspaceMultishift(RNG5,HermOpEO,
     //					0.0003,1.0e-5,2000); // Lo, tol, maxit
-    //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500); <== last run
-    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); // 176 with refinement
+    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500); <== last run
+    //    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); // 176 with refinement
     //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.001,3000,1500,200,0.0); // Attempt to resurrect
     SaveBasis(Aggregates,subspace_file);
   }

tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc

@@ -41,8 +41,8 @@ void SaveBasis(aggregation &Agg,std::string file)
   ScidacWriter WR(Agg.FineGrid->IsBoss());
   WR.open(file);
   for(int b=0;b<Agg.subspace.size();b++){
-    //    WR.writeScidacFieldRecord(Agg.subspace[b],record,0,BINARYIO_LEXICOGRAPHIC);
-    WR.writeScidacFieldRecord(Agg.subspace[b],record);
+    WR.writeScidacFieldRecord(Agg.subspace[b],record,0,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+    //    WR.writeScidacFieldRecord(Agg.subspace[b],record);
   }
   WR.close();
 #endif
@@ -55,8 +55,8 @@ void LoadBasis(aggregation &Agg, std::string file)
   ScidacReader RD ;
   RD.open(file);
   for(int b=0;b<Agg.subspace.size();b++){
-    //RD.readScidacFieldRecord(Agg.subspace[b],record,BINARYIO_LEXICOGRAPHIC);
-    RD.readScidacFieldRecord(Agg.subspace[b],record,0);
+    RD.readScidacFieldRecord(Agg.subspace[b],record,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+    //    RD.readScidacFieldRecord(Agg.subspace[b],record,0);
   }    
   RD.close();
 #endif
@@ -160,7 +160,7 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
   // Construct a coarsened grid with 4^4 cell
-  Coordinate Block({4,4,6,4});
+  Coordinate Block({4,4,4,4});
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/Block[d];
@@ -210,14 +210,14 @@ int main (int argc, char ** argv)
   // Need to check about red-black grid coarsening
   ////////////////////////////////////////////////////////////
 
-  //  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.mixed.2500.60");
-  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.new.62");
-  std::string refine_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Refine.phys48.tmp.60");
+  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.mixed.2500.60");
+  //  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.new.62");
+  std::string refine_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Refine.phys48.mixed.2500.60");
   std::string ldop_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/LittleDiracOp.phys48.mixed.60");
   std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/evecs.scidac");
   std::string eval_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/eval.xml");
   bool load_agg=true;
-  bool load_refine=false;
+  bool load_refine=true;
   bool load_mat=false;
   bool load_evec=false;
 
@@ -227,7 +227,10 @@ int main (int argc, char ** argv)
       LoadBasis(Aggregates,subspace_file);
     }
   } else {
-    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); // 176 with refinement
+    //    Aggregates.CreateSubspaceMultishift(RNG5,HermOpEO,
+    //    					0.0003,1.0e-5,2000); // Lo, tol, maxit
+    //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500);// <== last run
+    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); 
     SaveBasis(Aggregates,subspace_file);
   }
 
@@ -236,7 +239,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
   ConjugateGradient<CoarseVector>  coarseCG(4.0e-2,20000,true);
     
-  const int nrhs=8;
+  const int nrhs=12;
     
   Coordinate mpi=GridDefaultMpi();
   Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
@@ -366,7 +369,7 @@ int main (int argc, char ** argv)
   HDCGmrhs(src_mrhs,res_mrhs);
 
   // Standard CG
-#if 0
+#if 1
   {
   std::cout << "**************************************"<<std::endl;
   std::cout << "Calling red black CG"<<std::endl;