Further tuning and lanczos

Grid/algorithms/CoarsenedMatrix.h

@@ -323,7 +323,7 @@ public:
 
     // New normalised noise
     std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
-    std::cout << GridLogMessage<<" Chebyshev subspace pure noise  : nbasis"<<nn<<std::endl;
+    std::cout << GridLogMessage<<" Chebyshev subspace pure noise  : nbasis "<<nn<<std::endl;
 
 
     for(int b =0;b<nbasis;b++)
@@ -333,7 +333,8 @@ public:
       noise=noise*scale;
 
       // Initial matrix element
-      hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
+      hermop.Op(noise,Mn);
+      if(b==0) std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
       // Filter
       Chebyshev<FineField> Cheb(lo,hi,orderfilter);
       Cheb(hermop,noise,Mn);

Grid/algorithms/GeneralCoarsenedMatrix.h

@@ -220,6 +220,30 @@ public:
   GridBase      * FineGrid(void)       { return _FineGrid; };   // this is all the linalg routines need to know
   GridCartesian * CoarseGrid(void)     { return _CoarseGrid; };   // this is all the linalg routines need to know
 
+
+  void ProjectNearestNeighbour(RealD shift)
+  {
+    int Nd = geom.grid->Nd();
+    int point;
+    std::cout << "ProjectNearestNeighbour "<<std::endl;
+    for(int p=0;p<geom.npoint;p++){
+      int nhops = 0;
+      for(int s=0;s<Nd;s++){
+	nhops+=abs(geom.shifts[p][s]);
+      }
+      if(nhops>1) {
+	std::cout << "setting geom "<<p<<" shift "<<geom.shifts[p]<<" to zero "<<std::endl;
+	_A[p]=Zero();
+	_Adag[p]=Zero();
+      }
+      if(nhops==0) {
+	std::cout << " Adding IR shift "<<shift<<" to "<<geom.shifts[p]<<std::endl;
+	_A[p]=_A[p]+shift;
+	_Adag[p]=_Adag[p]+shift;
+      }
+    }
+  }
+  
   GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid)
     : geom(_geom),
       _FineGrid(FineGrid),
@@ -255,12 +279,19 @@ public:
   }
   void Mult (std::vector<CoarseMatrix> &A,const CoarseVector &in, CoarseVector &out)
   {
+    RealD ttot=0;
+    RealD tmult=0;
+    RealD texch=0;
+    RealD text=0;
+    ttot=-usecond();
     conformable(CoarseGrid(),in.Grid());
     conformable(in.Grid(),out.Grid());
     out.Checkerboard() = in.Checkerboard();
     CoarseVector tin=in;
 
+    texch-=usecond();
     CoarseVector pin  = Cell.Exchange(tin);
+    texch+=usecond();
 
     CoarseVector pout(pin.Grid());
     
@@ -281,11 +312,15 @@ public:
     typedef CComplex   calcComplex;
     
     int osites=pin.Grid()->oSites();
+    int gsites=pin.Grid()->gSites();
 
+    RealD flops = 1.0* npoint * nbasis * nbasis * 8 * gsites;
+      
     for(int point=0;point<npoint;point++){
       conformable(A[point],pin);
     }
     
+    tmult-=usecond();
     accelerator_for(sss, osites*nbasis, 1, {
       int ss = sss/nbasis;
       int b  = sss%nbasis;
@@ -313,10 +348,19 @@ public:
       }
       out_v[ss](b)=res;
     });
+    tmult+=usecond();
 
     for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
-
+    text-=usecond();
     out = Cell.Extract(pout);
+    text+=usecond();
+    ttot+=usecond();
+    std::cout << GridLogMessage<<"Coarse Mult exch "<<texch<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult mult "<<tmult<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult ext  "<<text<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Mult tot  "<<ttot<<" us"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse Kernel flops/s "<< flops/tmult<<" mflop/s"<<std::endl;
+    std::cout << GridLogMessage<<"Coarse flops/s "<< flops/ttot<<" mflop/s"<<std::endl;
   };
 
   void PopulateAdag(void)

Grid/algorithms/iterative/AdefGeneric.h

@@ -83,8 +83,8 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     coarsegrid = Aggregates.CoarseGrid;
     grid       = Aggregates.FineGrid;
   };
-
+  
-  void Inflexible(Field &src,Field &psi)
+  void Inflexible(const Field &src,Field &psi)
   {
     Field resid(grid);
     RealD f;
@@ -99,11 +99,13 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     Field r  (grid);
     Field mu (grid);
     Field rp (grid);
-
+    
     //Initial residual computation & set up
     RealD guess = norm2(psi);
     double tn;
 
+    GridStopWatch HDCGTimer;
+    HDCGTimer.Start();
     //////////////////////////
     // x0 = Vstart -- possibly modify guess
     //////////////////////////
@@ -168,7 +170,8 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
       // Stopping condition
       if ( rn <= rsq ) { 
 
-	std::cout<<GridLogMessage<<"HDCG: Pcg converged in "<<k<<" iterations"<<std::endl;;
+	HDCGTimer.Stop();
+	std::cout<<GridLogMessage<<"HDCG: Pcg converged in "<<k<<" iterations and "<<HDCGTimer.Elapsed()<<std::endl;;
 
 	_FineLinop.HermOp(x,mmp);			  
 	axpy(tmp,-1.0,src,mmp);
@@ -189,126 +192,9 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     return ;
   }
   
-  // The Pcg routine is common to all, but the various matrices differ from derived 
+  virtual void operator() (const Field &in, Field &out)
-  // implementation to derived implmentation
+  {
-  void operator() (const Field &src, Field &psi){
+    this->Inflexible(in,out);
-
-    psi.Checkerboard() = src.Checkerboard();
-    grid             = src.Grid();
-
-    RealD f;
-    RealD rtzp,rtz,a,d,b;
-    RealD rptzp;
-    RealD tn;
-    RealD guess = norm2(psi);
-    RealD ssq   = norm2(src);
-    RealD rsq   = ssq*Tolerance*Tolerance;
-    
-    /////////////////////////////
-    // Set up history vectors
-    /////////////////////////////
-    std::vector<Field> p  (mmax,grid);
-    std::vector<Field> mmp(mmax,grid);
-    std::vector<RealD> pAp(mmax);
-
-    Field x  (grid);
-    Field z  (grid);
-    Field tmp(grid);
-    Field r  (grid);
-    Field mu (grid);
-  
-    //////////////////////////
-    // x0 = Vstart -- possibly modify guess
-    //////////////////////////
-    x=Zero();
-    Vstart(x,src);
-
-    // r0 = b -A x0
-    _FineLinop.HermOp(x,mmp[0]); // Fine operator
-    axpy (r, -1.0,mmp[0], src);    // Recomputes r=src-Ax0
-
-    //////////////////////////////////
-    // Compute z = M1 r
-    //////////////////////////////////
-    PcgM1(r,z);
-    rtzp =real(innerProduct(r,z));
-
-    ///////////////////////////////////////
-    // Solve for Mss mu = P A z and set p = z-mu
-    ///////////////////////////////////////
-    PcgM2(z,p[0]);
-
-    for (int k=0;k<=MaxIterations;k++){
-    
-      int peri_k  = k % mmax;
-      int peri_kp = (k+1) % mmax;
-
-      rtz=rtzp;
-      d= PcgM3(p[peri_k],mmp[peri_k]);
-      a = rtz/d;
-    
-      // Memorise this
-      pAp[peri_k] = d;
-      std::cout << GridLogMessage << " pCG d "<< d<<std::endl;
-
-      axpy(x,a,p[peri_k],x);
-      //      std::cout << GridLogMessage << " pCG x "<< norm2(x)<<std::endl;
-      RealD rn = axpy_norm(r,-a,mmp[peri_k],r);
-
-      std::cout << GridLogMessage << " pCG rn "<< rn<<std::endl;
-
-      // Compute z = M x
-      PcgM1(r,z);
-      //      std::cout << GridLogMessage << " pCG z "<< norm2(z)<<std::endl;
-
-      rtzp =real(innerProduct(r,z));
-      std::cout << GridLogMessage << " pCG rtzp "<<rtzp<<std::endl;
-      //      std::cout << GridLogMessage << " pCG r "<<norm2(r)<<std::endl;
-
-      PcgM2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
-
-      //      std::cout << GridLogMessage << " pCG mu "<<norm2(mu)<<std::endl;
-      
-      p[peri_kp]=mu;
-
-      //      std::cout << GridLogMessage << " pCG p[peri_kp] "<<norm2(p[peri_kp])<<std::endl;
-
-      // Standard search direction  p -> z + b p 
-      b = (rtzp)/rtz;
-      std::cout << GridLogMessage << " pCG b "<< b<<std::endl;
-
-      int northog;
-      //    northog     = (peri_kp==0)?1:peri_kp; // This is the fCG(mmax) algorithm
-      northog     = (k>mmax-1)?(mmax-1):k;        // This is the fCG-Tr(mmax-1) algorithm
-    
-      for(int back=0; back < northog; back++){
-	int peri_back = (k-back)%mmax;
-	RealD pbApk= real(innerProduct(mmp[peri_back],p[peri_kp]));
-	RealD beta = -pbApk/pAp[peri_back];
-	axpy(p[peri_kp],beta,p[peri_back],p[peri_kp]);
-      }
-      //      std::cout << GridLogMessage << " pCG p[peri_kp] orthog "<< norm2(p[peri_kp])<<std::endl;
-
-      RealD rrn=sqrt(rn/ssq);
-      std::cout<<GridLogMessage<<"TwoLevelfPcg: k= "<<k<<" residual = "<<rrn<<std::endl;
-
-      // Stopping condition
-      if ( rn <= rsq ) { 
-
-	_FineLinop.HermOp(x,mmp[0]); // Shouldn't this be something else?
-	axpy(tmp,-1.0,src,mmp[0]);
-	
-	RealD psinorm = sqrt(norm2(x));
-	RealD srcnorm = sqrt(norm2(src));
-	RealD tmpnorm = sqrt(norm2(tmp));
-	RealD true_residual = tmpnorm/srcnorm;
-	std::cout<<GridLogMessage<<"TwoLevelfPcg:   true residual is "<<true_residual<<std::endl;
-	std::cout<<GridLogMessage<<"TwoLevelfPcg: target residual was"<<Tolerance<<std::endl;
-	return;
-      }
-    }
-    // Non-convergence
-    assert(0);
   }
 
  public:
@@ -322,17 +208,37 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
     CoarseField PleftProj(coarsegrid);
     CoarseField PleftMss_proj(coarsegrid);
 
+    GridStopWatch SmootherTimer;
+    GridStopWatch MatrixTimer;
+    SmootherTimer.Start();
     _Smoother(in,Min);
+    SmootherTimer.Stop();
 
+    MatrixTimer.Start();
     _FineLinop.HermOp(Min,out);
+    MatrixTimer.Stop();
     axpy(tmp,-1.0,out,in);          // tmp  = in - A Min
 
+    GridStopWatch ProjTimer;
+    GridStopWatch CoarseTimer;
+    GridStopWatch PromTimer;
+    ProjTimer.Start();
     _Aggregates.ProjectToSubspace(PleftProj,tmp);     
+    ProjTimer.Stop();
+    CoarseTimer.Start();
     _CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
+    CoarseTimer.Stop();
+    PromTimer.Start();
     _Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]  
+    PromTimer.Stop();
+    std::cout << GridLogMessage << "PcgM1 breakdown "<<std::endl;
+    std::cout << GridLogMessage << "\tSmoother   " << SmootherTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tProj       " << ProjTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tCoarse     " << CoarseTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tProm       " << PromTimer.Elapsed() <<std::endl;
 
     axpy(out,1.0,Min,tmp); // Min+tmp
-
   }
 
   virtual void PcgM2(const Field & in, Field & out) {

tests/debug/Test_general_coarse_hdcg.cc

@@ -177,17 +177,21 @@ int main (int argc, char ** argv)
   LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d);
   LittleDiracOp.CoarsenOperatorColoured(FineHermOp,Aggregates);
 
+  // Try projecting to one hop only
+  LittleDiracOperator LittleDiracOpProj(LittleDiracOp);
+  LittleDiracOpProj.ProjectNearestNeighbour(0.5);
+
   typedef HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermMatrix;
   HermMatrix CoarseOp (LittleDiracOp);
 
   //////////////////////////////////////////
   // Build a coarse lanczos
   //////////////////////////////////////////
-  Chebyshev<CoarseVector>      IRLCheby(0.02,50.0,71);  // 1 iter
+  Chebyshev<CoarseVector>      IRLCheby(0.5,60.0,71);  // 1 iter
   FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,CoarseOp);
   PlainHermOp<CoarseVector>    IRLOp    (CoarseOp);
-  int Nk=64;
+  int Nk=48;
-  int Nm=128;
+  int Nm=64;
   int Nstop=Nk;
   ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-5,20);
 
@@ -195,6 +199,9 @@ int main (int argc, char ** argv)
   std::vector<RealD>            eval(Nm);
   std::vector<CoarseVector>     evec(Nm,Coarse5d);
   CoarseVector c_src(Coarse5d); c_src=1.0;
+
+  PowerMethod<CoarseVector>       cPM;   cPM(CoarseOp,c_src);
+
   IRL.calc(eval,evec,c_src,Nconv);
   DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
   
@@ -230,7 +237,9 @@ int main (int argc, char ** argv)
   //         use a limited stencil. Reread BFM code to check on evecs / deflation strategy with prec
   //
   std::vector<RealD> los({3.0}); // Nbasis 40 == 36,36 iters
-  std::vector<int> ords({7,8,10}); // Nbasis 40 == 40,38,36 iters (320,342,396 mults)
+
+  //  std::vector<int> ords({7,8,10}); // Nbasis 40 == 40,38,36 iters (320,342,396 mults)
+  std::vector<int> ords({7}); // Nbasis 40 == 40 iters (320 mults)  
 
   // Standard CG
   //      result=Zero();

tests/debug/Test_general_coarse_pvdagm.cc

@@ -189,7 +189,7 @@ int main (int argc, char ** argv)
   typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
   typedef LittleDiracOperator::CoarseVector CoarseVector;
 
-  NextToNearestStencilGeometry5D geom;
+  NextToNearestStencilGeometry5D geom(Coarse5d);
 
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;