HDCR running on 16^3 with 2x-3x speed up.

lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h

@@ -47,6 +47,10 @@ namespace Grid {
     int mmax;
     int nstep;
     int steps;
+    GridStopWatch PrecTimer;
+    GridStopWatch MatTimer;
+    GridStopWatch LinalgTimer;
+
     LinearFunction<Field> &Preconditioner;
 
    PrecGeneralisedConjugateResidual(RealD tol,Integer maxit,LinearFunction<Field> &Prec,int _mmax,int _nstep) : 
@@ -68,6 +72,10 @@ namespace Grid {
       
       Field r(src._grid);
 
+        PrecTimer.Reset();
+         MatTimer.Reset();
+      LinalgTimer.Reset();
+
       GridStopWatch SolverTimer;
       SolverTimer.Start();
 
@@ -76,7 +84,7 @@ namespace Grid {
 
 	cp=GCRnStep(Linop,src,psi,rsq);
 
-	if ( verbose ) std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl;
+	std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl;
 
 	if(cp<rsq) {
 
@@ -89,7 +97,11 @@ namespace Grid {
 		   << " computed residual "<<sqrt(cp/ssq)
 	           << " true residual "    <<sqrt(tr/ssq)
 	           << " target "           <<Tolerance <<std::endl;
-	  std::cout<<GridLogMessage<<" Time elapsed: Total "<< SolverTimer.Elapsed() <<std::endl;
+
+	  std::cout<<GridLogMessage<<"VPGCR Time elapsed: Total  "<< SolverTimer.Elapsed() <<std::endl;
+	  std::cout<<GridLogMessage<<"VPGCR Time elapsed: Precon "<<   PrecTimer.Elapsed() <<std::endl;
+	  std::cout<<GridLogMessage<<"VPGCR Time elapsed: Matrix "<<    MatTimer.Elapsed() <<std::endl;
+	  std::cout<<GridLogMessage<<"VPGCR Time elapsed: Linalg "<< LinalgTimer.Elapsed() <<std::endl;
 	  return;
 	}
 
@@ -97,6 +109,7 @@ namespace Grid {
       std::cout<<GridLogMessage<<"Variable Preconditioned GCR did not converge"<<std::endl;
       assert(0);
     }
+
     RealD GCRnStep(LinearOperatorBase<Field> &Linop,const Field &src, Field &psi,RealD rsq){
 
       RealD cp;
@@ -123,24 +136,25 @@ namespace Grid {
       // initial guess x0 is taken as nonzero.
       // r0=src-A x0 = src
       //////////////////////////////////
+      MatTimer.Start();
       Linop.HermOpAndNorm(psi,Az,zAz,zAAz); 
+      MatTimer.Stop();
       r=src-Az;
       
       /////////////////////
       // p = Prec(r)
       /////////////////////
+      PrecTimer.Start();
       Preconditioner(r,z);
+      PrecTimer.Stop();
 
-      std::cout<<GridLogMessage<< " Preconditioner in " << norm2(r)<<std::endl; 
+      MatTimer.Start();
-      std::cout<<GridLogMessage<< " Preconditioner out " << norm2(z)<<std::endl; 
-      
       Linop.HermOp(z,tmp); 
+      MatTimer.Stop();
 
-      std::cout<<GridLogMessage<< " Preconditioner Aout " << norm2(tmp)<<std::endl; 
       ttmp=tmp;
       tmp=tmp-r;
 
-      std::cout<<GridLogMessage<< " Preconditioner resid " << std::sqrt(norm2(tmp)/norm2(r))<<std::endl; 
       /*
       std::cout<<GridLogMessage<<r<<std::endl;
       std::cout<<GridLogMessage<<z<<std::endl;
@@ -148,7 +162,9 @@ namespace Grid {
       std::cout<<GridLogMessage<<tmp<<std::endl;
       */
 
+      MatTimer.Start();
       Linop.HermOpAndNorm(z,Az,zAz,zAAz); 
+      MatTimer.Stop();
 
       //p[0],q[0],qq[0] 
       p[0]= z;
@@ -172,18 +188,22 @@ namespace Grid {
 
 	cp = axpy_norm(r,-a,q[peri_k],r);  
 
-	std::cout<<GridLogMessage<< " VPGCR_step resid" <<sqrt(cp/rsq)<<std::endl; 
 	if((k==nstep-1)||(cp<rsq)){
 	  return cp;
 	}
 
+	std::cout<<GridLogMessage<< " VPGCR_step     resid " <<sqrt(cp/rsq)<<std::endl; 
+
+	PrecTimer.Start();
 	Preconditioner(r,z);// solve Az = r
+	PrecTimer.Stop();
+
+	MatTimer.Start();
 	Linop.HermOpAndNorm(z,Az,zAz,zAAz);
-
-
 	Linop.HermOp(z,tmp);
+	MatTimer.Stop();
         tmp=tmp-r;
-	std::cout<<GridLogMessage<< " Preconditioner resid" <<sqrt(norm2(tmp)/norm2(r))<<std::endl; 
+	std::cout<<GridLogMessage<< " Preconditioner resid " <<sqrt(norm2(tmp)/norm2(r))<<std::endl; 
 
 	q[peri_kp]=Az;
 	p[peri_kp]=z;

tests/Test_dwf_hdcr.cc

@@ -263,11 +263,6 @@ public:
       resid = norm2(r) /norm2(src); 
       std::cout << "SAP "<<i<<" resid "<<resid<<std::endl;
 
-
-// Npoly*outer*2 1/2 vol matmuls.
-// 71 iters => 20*71 = 1400 matmuls.
-// 2*71 = 140 comms.
-
       // Even domain solve
       r= where(subset==(Integer)0,r,zz);
       _SmootherOperator.AdjOp(r,vec1);
@@ -332,7 +327,7 @@ public:
     CoarseVector Ctmp(_CoarseOperator.Grid());
     CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
 
-    ConjugateGradient<CoarseVector>  CG(1.0e-3,100000);
+    ConjugateGradient<CoarseVector>  CG(1.0e-2,100000);
     //    ConjugateGradient<FineField>    fCG(3.0e-2,1000);
 
     HermitianLinearOperator<CoarseOperator,CoarseVector>  HermOp(_CoarseOperator);
@@ -345,14 +340,14 @@ public:
 
     //    Chebyshev<FineField> Cheby    (0.5,70.0,30,InverseApproximation);
     //    Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
-    Chebyshev<FineField> Cheby    (2.0,70.0,15,InverseApproximation);
+    Chebyshev<FineField> Cheby    (params.lo,params.hi,params.order,InverseApproximation);
-    Chebyshev<FineField> ChebyAccu(2.0,70.0,15,InverseApproximation);
+    Chebyshev<FineField> ChebyAccu(params.lo,params.hi,params.order,InverseApproximation);
     //    Cheby.JacksonSmooth();
     //    ChebyAccu.JacksonSmooth();
 
-    _Aggregates.ProjectToSubspace  (Csrc,in);
+    //    _Aggregates.ProjectToSubspace  (Csrc,in);
-    _Aggregates.PromoteFromSubspace(Csrc,out);
+    //    _Aggregates.PromoteFromSubspace(Csrc,out);
-    std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
+    //    std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
     
     //    ofstream fout("smoother");
     //    Cheby.csv(fout);
@@ -479,7 +474,7 @@ int main (int argc, char ** argv)
   read(RD,"params",params);
   std::cout<<"Params: Order "<<params.order<<"["<<params.lo<<","<<params.hi<<"]"<< " steps "<<params.steps<<std::endl;
 
-  const int Ls=8;
+  const int Ls=16;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
@@ -490,10 +485,12 @@ int main (int argc, char ** argv)
   ///////////////////////////////////////////////////
   // Construct a coarsened grid; utility for this?
   ///////////////////////////////////////////////////
-  const int block=2;
+  std::vector<int> block ({2,2,2,2});
+  const int nbasis= 32;
+
   std::vector<int> clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/block;
+    clatt[d] = clatt[d]/block[d];
   }
   GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
   GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
@@ -539,7 +536,7 @@ int main (int argc, char ** argv)
   //  SU3::HotConfiguration(RNG4,Umu);
   //  Umu=zero;
 
-  RealD mass=0.01;
+  RealD mass=0.001;
   RealD M5=1.8;
 
   std::cout<<GridLogMessage << "**************************************************"<< std::endl;
@@ -548,9 +545,6 @@ int main (int argc, char ** argv)
   DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   DomainWallFermionR DdwfDD(UmuDD,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
 
-  const int nbasis = 16;
-  //  const int nbasis = 4;
-
   typedef Aggregation<vSpinColourVector,vTComplex,nbasis>              Subspace;
   typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>          CoarseOperator;
   typedef CoarseOperator::CoarseVector                                 CoarseVector;