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

lib/algorithms/LinearOperator.h

@@ -308,20 +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);
-	_Mat.Meooe(in,tmp);
+	Field tmp2(in._grid);
-	_Mat.MooeeInv(tmp,out);
+
-	_Mat.Meooe(out,tmp);
+	std::cout << GridLogIterative << " HermOp.Mpc "<<std::endl;
 	_Mat.Mooee(in,out);
-        return axpy_norm(out,-1.0,tmp,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;
+
+	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);

lib/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -168,8 +168,8 @@ void basisDeflate(const std::vector<Field> &_v,const std::vector<RealD>& eval,co
 template<class Field> class ImplicitlyRestartedLanczosTester 
 {
  public:
-  virtual int TestConvergence(int j,RealD resid,Field &evec, RealD &eval,RealD evalMaxApprox);
+  virtual int TestConvergence(int j,RealD resid,Field &evec, RealD &eval,RealD evalMaxApprox)=0;
-  virtual int ReconstructEval(int j,RealD resid,Field &evec, RealD &eval,RealD evalMaxApprox);
+  virtual int ReconstructEval(int j,RealD resid,Field &evec, RealD &eval,RealD evalMaxApprox)=0;
 };
 
 enum IRLdiagonalisation { 

lib/algorithms/iterative/SchurRedBlack.h

@@ -90,7 +90,7 @@ namespace Grid {
   // Take a matrix and form a Red Black solver calling a Herm solver
   // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
-
+  // Now make the norm reflect extra factor of Mee
   template<class Field> class SchurRedBlackStaggeredSolve {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
@@ -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)
@@ -136,14 +139,15 @@ namespace Grid {
       _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
       tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
 
-      src_o = tmp;     assert(src_o.checkerboard ==Odd);
+      //src_o = tmp;     assert(src_o.checkerboard ==Odd);
-      //  _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source
+      _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm.
 
       //////////////////////////////////////////////////////////////
       // Call the red-black solver
       //////////////////////////////////////////////////////////////
       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;
     }     
   };

lib/communicator/Communicator_mpi.cc

@@ -57,7 +57,7 @@ CartesianCommunicator::~CartesianCommunicator()
 {
   int MPI_is_finalised;
   MPI_Finalized(&MPI_is_finalised);
-  if (communicator && MPI_is_finalised)
+  if (communicator && !MPI_is_finalised)
     MPI_Comm_free(&communicator);
 }
 

lib/communicator/Communicator_mpi3.cc

@@ -604,7 +604,7 @@ CartesianCommunicator::~CartesianCommunicator()
 {
   int MPI_is_finalised;
   MPI_Finalized(&MPI_is_finalised);
-  if (communicator && MPI_is_finalised) {
+  if (communicator && !MPI_is_finalised) {
     MPI_Comm_free(&communicator);
     for(int i=0;i<communicator_halo.size();i++){
       MPI_Comm_free(&communicator_halo[i]);

lib/communicator/Communicator_mpit.cc

@@ -57,7 +57,7 @@ CartesianCommunicator::~CartesianCommunicator()
 {
   int MPI_is_finalised;
   MPI_Finalized(&MPI_is_finalised);
-  if (communicator && MPI_is_finalised) {
+  if (communicator && !MPI_is_finalised){
     MPI_Comm_free(&communicator);
     for(int i=0;i<  communicator_halo.size();i++){
       MPI_Comm_free(&communicator_halo[i]);
@@ -65,7 +65,6 @@ CartesianCommunicator::~CartesianCommunicator()
   }  
 }
 
-
 void CartesianCommunicator::GlobalSum(uint32_t &u){
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
   assert(ierr==0);
@@ -247,7 +246,7 @@ void CartesianCommunicator::StencilSendToRecvFromComplete(std::vector<CommsReque
 { 
   int nreq=waitall.size();
   MPI_Waitall(nreq, &waitall[0], MPI_STATUSES_IGNORE);
-};
+}
 double CartesianCommunicator::StencilSendToRecvFrom(void *xmit,
 						    int xmit_to_rank,
 						    void *recv,

lib/communicator/Communicator_shmem.cc

@@ -75,6 +75,8 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
   ShmInitGeneric();
 }
 
+CartesianCommunicator::~CartesianCommunicator(){}
+
 CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors,const CartesianCommunicator &parent) 
   : CartesianCommunicator(processors) 
 {

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
+    parallel_for(int ss=0;ss<full._grid->oSites();ss++){
-    for(int ss=0;ss<full._grid->oSites();ss++){
-      std::vector<int> coor;
       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(int ss=0;ss<full._grid->oSites();ss++){
-    //parallel_for
-    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++;
       }
     }
   }

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();
 }