From 78e8704eacb41fae706e50c24ae0baa6b17b9481 Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Mon, 30 Oct 2017 00:25:31 +0000
Subject: [PATCH] Shaking out

---
 tests/solver/Test_dwf_mrhs_cg_mpi.cc | 99 +++++++++++++++++++++++++---
 1 file changed, 89 insertions(+), 10 deletions(-)

diff --git a/tests/solver/Test_dwf_mrhs_cg_mpi.cc b/tests/solver/Test_dwf_mrhs_cg_mpi.cc
index fbc6dd32..f640edff 100644
--- a/tests/solver/Test_dwf_mrhs_cg_mpi.cc
+++ b/tests/solver/Test_dwf_mrhs_cg_mpi.cc
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+   /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -47,20 +47,36 @@ int main (int argc, char ** argv)
   std::vector<int> mpi_layout  = GridDefaultMpi();
   std::vector<int> mpi_split (mpi_layout.size(),1);
 
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * rbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
-  int nrhs = UGrid->RankCount() ;
-
   /////////////////////////////////////////////
   // Split into 1^4 mpi communicators
   /////////////////////////////////////////////
+
+  for(int i=0;i<argc;i++){
+    if(std::string(argv[i]) == "--split"){
+      for(int k=0;k<mpi_layout.size();k++){
+	std::stringstream ss; 
+	ss << argv[i+1+k]; 
+	ss >> mpi_split[k];
+      }
+      break;
+    }
+  }
+
+  int nrhs = 1;
+  int me;
+  for(int i=0;i<mpi_layout.size();i++) nrhs *= (mpi_layout[i]/mpi_split[i]);
+
   GridCartesian         * SGrid = new GridCartesian(GridDefaultLatt(),
 						    GridDefaultSimd(Nd,vComplex::Nsimd()),
 						    mpi_split,
-						    *UGrid); 
+						    *UGrid,me); 
 
   GridCartesian         * SFGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,SGrid);
   GridRedBlackCartesian * SrbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(SGrid);
@@ -78,16 +94,46 @@ int main (int argc, char ** argv)
   std::vector<FermionField> result(nrhs,FGrid);
   FermionField tmp(FGrid);
 
-  for(int s=0;s<nrhs;s++) random(pRNG5,src[s]);
   for(int s=0;s<nrhs;s++) result[s]=zero;
+#undef LEXICO_TEST
+#ifdef LEXICO_TEST
+  {
+    LatticeFermion lex(FGrid);  lex = zero;
+    LatticeFermion ftmp(FGrid);
+    Integer stride =10000;
+    double nrm;
+    LatticeComplex coor(FGrid);
+    for(int d=0;d<5;d++){
+      LatticeCoordinate(coor,d);
+      ftmp = stride;
+      ftmp = ftmp * coor;
+      lex = lex + ftmp;
+      stride=stride/10;
+    }
+    for(int s=0;s<nrhs;s++) {
+      src[s]=lex;
+      ftmp = 1000*1000*s;
+      src[s] = src[s] + ftmp;
+    }    
+  }
+#else
+  for(int s=0;s<nrhs;s++) {
+    random(pRNG5,src[s]);
+    tmp = 100.0*s;
+    src[s] = (src[s] * 0.1) + tmp;
+    std::cout << " src ]"<<s<<"] "<<norm2(src[s])<<std::endl;
+  }
+#endif
+
+  for(int n =0 ; n< nrhs ; n++) { 
+    std::cout << " src"<<n<<"\n"<< src[n] <<std::endl;
+  }
 
   LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu);
 
   /////////////////
   // MPI only sends
   /////////////////
-  int me = UGrid->ThisRank();
-
   LatticeGaugeField s_Umu(SGrid);
   FermionField s_src(SFGrid);
   FermionField s_tmp(SFGrid);
@@ -98,6 +144,36 @@ int main (int argc, char ** argv)
   ///////////////////////////////////////////////////////////////
   Grid_split  (Umu,s_Umu);
   Grid_split  (src,s_src);
+  std::cout << " split rank  " <<me << " s_src "<<norm2(s_src)<<std::endl;
+  std::cout << " s_src\n "<< s_src <<std::endl;
+
+#ifdef LEXICO_TEST
+  FermionField s_src_tmp(SFGrid);
+  FermionField s_src_diff(SFGrid);
+  {
+    LatticeFermion lex(SFGrid);  lex = zero;
+    LatticeFermion ftmp(SFGrid);
+    Integer stride =10000;
+    double nrm;
+    LatticeComplex coor(SFGrid);
+    for(int d=0;d<5;d++){
+      LatticeCoordinate(coor,d);
+      ftmp = stride;
+      ftmp = ftmp * coor;
+      lex = lex + ftmp;
+      stride=stride/10;
+    }
+    s_src_tmp=lex;
+    ftmp = 1000*1000*me;
+    s_src_tmp = s_src_tmp + ftmp;
+  }
+  s_src_diff = s_src_tmp - s_src;
+  std::cout << " s_src_diff " << norm2(s_src_diff)<<std::endl;
+
+  std::cout << " s_src \n" << s_src << std::endl;
+  std::cout << " s_src_tmp \n" << s_src_tmp << std::endl;
+  std::cout << " s_src_diff \n" << s_src_diff << std::endl;
+#endif
 
   ///////////////////////////////////////////////////////////////
   // Set up N-solvers as trivially parallel
@@ -113,10 +189,11 @@ int main (int argc, char ** argv)
 
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOpCk(Dchk);
-  ConjugateGradient<FermionField> CG((1.0e-5/(me+1)),10000);
+  ConjugateGradient<FermionField> CG((1.0e-5),10000);
   s_res = zero;
   CG(HermOp,s_src,s_res);
 
+  std::cout << " s_res norm "<<norm2(s_res)<<std::endl;
   /////////////////////////////////////////////////////////////
   // Report how long they all took
   /////////////////////////////////////////////////////////////
@@ -134,10 +211,12 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage<< "Unsplitting the result"<<std::endl;
   Grid_unsplit(result,s_res);
 
+
   std::cout << GridLogMessage<< "Checking the residuals"<<std::endl;
   for(int n=0;n<nrhs;n++){
+    std::cout << " res["<<n<<"] norm "<<norm2(result[n])<<std::endl;
     HermOpCk.HermOp(result[n],tmp); tmp = tmp - src[n];
-    std::cout << GridLogMessage<<" resid["<<n<<"]  "<< norm2(tmp)<<std::endl;
+    std::cout << GridLogMessage<<" resid["<<n<<"]  "<< norm2(tmp)/norm2(src[n])<<std::endl;
   }
 
   Grid_finalize();