From e55397bc134ead26dbac8e2ef244406a8b9d6a3b Mon Sep 17 00:00:00 2001
From: azusayamaguchi <ayamaguc@staffmail.ed.ac.uk>
Date: Fri, 24 Nov 2017 14:18:30 +0000
Subject: [PATCH] Staggerd cg

---
 lib/algorithms/LinearOperator.h          | 20 +++++++++++---------
 lib/algorithms/iterative/SchurRedBlack.h |  9 +++++++--
 lib/lattice/Lattice_transfer.h           | 16 ++++++----------
 tests/solver/Test_staggered_cg_schur.cc  | 14 ++++++++++++++
 4 files changed, 38 insertions(+), 21 deletions(-)

diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h
index 0fa039c8..26746e6e 100644
--- a/lib/algorithms/LinearOperator.h
+++ b/lib/algorithms/LinearOperator.h
@@ -308,32 +308,34 @@ namespace Grid {
     public:
       SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){};
       virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+	GridLogIterative.TimingMode(1);
+	std::cout << GridLogIterative << " HermOpAndNorm "<<std::endl;
 	n2 = Mpc(in,out);
+	std::cout << GridLogIterative << " HermOpAndNorm.Mpc "<<std::endl;
 	ComplexD dot= innerProduct(in,out);
+	std::cout << GridLogIterative << " HermOpAndNorm.innerProduct "<<std::endl;
 	n1 = real(dot);
       }
       virtual void HermOp(const Field &in, Field &out){
+	std::cout << GridLogIterative << " HermOp "<<std::endl;
 	Mpc(in,out);
       }
       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;
 
 	_Mat.Meooe(in,out);
 	_Mat.Meooe(out,tmp2);
+	std::cout << GridLogIterative << " HermOp.MeooeMeooe "<<std::endl;
 
-	return axpy_norm(out,-1.0,tmp2,tmp);
-#if 0
-	//... much prefer conventional Schur norm
-	_Mat.Meooe(in,tmp);
-	_Mat.MooeeInv(tmp,out);
-	_Mat.Meooe(out,tmp);
-	_Mat.Mooee(in,out);
-        return axpy_norm(out,-1.0,tmp,out);
-#endif
+	RealD nn=axpy_norm(out,-1.0,tmp2,tmp);
+	std::cout << GridLogIterative << " HermOp.axpy_norm "<<std::endl;
+	return nn;
       }
       virtual  RealD MpcDag   (const Field &in, Field &out){
 	return Mpc(in,out);
diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h
index b9767aa8..5f5a8b66 100644
--- a/lib/algorithms/iterative/SchurRedBlack.h
+++ b/lib/algorithms/iterative/SchurRedBlack.h
@@ -123,11 +123,14 @@ namespace Grid {
       Field   tmp(grid);
       Field  Mtmp(grid);
       Field resid(fgrid);
-
+      
+      std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve " <<std::endl;
       pickCheckerboard(Even,src_e,in);
       pickCheckerboard(Odd ,src_o,in);
       pickCheckerboard(Even,sol_e,out);
       pickCheckerboard(Odd ,sol_o,out);
+
+      std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve checkerboards picked" <<std::endl;
     
       /////////////////////////////////////////////////////
       // src_o = (source_o - Moe MeeInv source_e)
@@ -144,6 +147,7 @@ namespace Grid {
       //////////////////////////////////////////////////////////////
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
       _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called  the Mpc solver" <<std::endl;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -152,15 +156,16 @@ namespace Grid {
       src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
       _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
      
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver reconstructed other CB" <<std::endl;
       setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
+      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;
     }     
   };
diff --git a/lib/lattice/Lattice_transfer.h b/lib/lattice/Lattice_transfer.h
index 1b09217b..43fdbe46 100644
--- a/lib/lattice/Lattice_transfer.h
+++ b/lib/lattice/Lattice_transfer.h
@@ -50,26 +50,22 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
   ////////////////////////////////////////////////////////////////////////////////////////////
   template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full){
     half.checkerboard = cb;
-    int ssh=0;
-    //parallel_for
-    for(int ss=0;ss<full._grid->oSites();ss++){
-      std::vector<int> coor;
+
+    parallel_for(int ss=0;ss<full._grid->oSites();ss++){
       int cbos;
-      
+      std::vector<int> coor;
       full._grid->oCoorFromOindex(coor,ss);
       cbos=half._grid->CheckerBoard(coor);
       
       if (cbos==cb) {
+	int ssh=half._grid->oIndex(coor);
 	half._odata[ssh] = full._odata[ss];
-	ssh++;
       }
     }
   }
   template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half){
     int cb = half.checkerboard;
-    int ssh=0;
-    //parallel_for
-    for(int ss=0;ss<full._grid->oSites();ss++){
+    parallel_for(int ss=0;ss<full._grid->oSites();ss++){
       std::vector<int> coor;
       int cbos;
 
@@ -77,8 +73,8 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
       cbos=half._grid->CheckerBoard(coor);
       
       if (cbos==cb) {
+	int ssh=half._grid->oIndex(coor);
 	full._odata[ss]=half._odata[ssh];
-	ssh++;
       }
     }
   }
diff --git a/tests/solver/Test_staggered_cg_schur.cc b/tests/solver/Test_staggered_cg_schur.cc
index 09044995..a5c25b85 100644
--- a/tests/solver/Test_staggered_cg_schur.cc
+++ b/tests/solver/Test_staggered_cg_schur.cc
@@ -70,7 +70,21 @@ int main (int argc, char ** argv)
   ConjugateGradient<FermionField> CG(1.0e-8,10000);
   SchurRedBlackStaggeredSolve<FermionField> SchurSolver(CG);
 
+  double volume=1.0;
+  for(int mu=0;mu<Nd;mu++){
+    volume=volume*latt_size[mu];
+  }  
+  double t1=usecond();
   SchurSolver(Ds,src,result);
+  double t2=usecond();
+
+  // Schur solver: uses DeoDoe => volume * 1146
+  double ncall=CG.IterationsToComplete;
+  double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 +  == 1146
+
+  std::cout<<GridLogMessage << "usec    =   "<< (t2-t1)<<std::endl;
+  std::cout<<GridLogMessage << "flop/s  =   "<< flops<<std::endl;
+  std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t2-t1)<<std::endl;
   
   Grid_finalize();
 }