From 7907e4ca03a2232d6db973659bd90caef3b86df7 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Sat, 11 Jul 2015 23:06:31 +0900
Subject: [PATCH 01/10] This file drives me crazy

---
 INSTALL | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/INSTALL b/INSTALL
index f812f5a2..80a61507 120000
--- a/INSTALL
+++ b/INSTALL
@@ -1 +1 @@
-/usr/share/automake-1.14/INSTALL
\ No newline at end of file
+/opt/local/share/automake-1.15/INSTALL
\ No newline at end of file

From 8a7b7f1e2b2a953f5ce81449b22f024158f45b05 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 12:12:29 +0900
Subject: [PATCH 02/10] verbose

---
 lib/algorithms/iterative/ConjugateGradient.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h
index 3829dad8..e281d765 100644
--- a/lib/algorithms/iterative/ConjugateGradient.h
+++ b/lib/algorithms/iterative/ConjugateGradient.h
@@ -15,7 +15,7 @@ public:
     Integer MaxIterations;
     int verbose;
     ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { 
-      verbose=1;
+      verbose=0;
     };
 
 

From c515d069cdef2584e33e590cd12acd6249313d95 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 12:13:03 +0900
Subject: [PATCH 03/10] Tweaks to subspace set up to put in g5 r5 hermiticity

---
 lib/algorithms/CoarsenedMatrix.h | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h
index e257294c..998292b2 100644
--- a/lib/algorithms/CoarsenedMatrix.h
+++ b/lib/algorithms/CoarsenedMatrix.h
@@ -117,15 +117,15 @@ namespace Grid {
       }
       Orthogonalise();
     }
-    virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop) {
+    virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
 
       RealD scale;
 
-      ConjugateGradient<FineField> CG(1.0e-4,10000);
+      ConjugateGradient<FineField> CG(2.0e-3,10000);
       FineField noise(FineGrid);
       FineField Mn(FineGrid);
 
-      for(int b=0;b<nbasis;b++){
+      for(int b=0;b<nn;b++){
 	
 	gaussian(RNG,noise);
 	scale = std::pow(norm2(noise),-0.5); 
@@ -133,7 +133,7 @@ namespace Grid {
 
 	hermop.Op(noise,Mn); std::cout << "noise   ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
 
-	for(int i=0;i<2;i++){
+	for(int i=0;i<1;i++){
 
 	  CG(hermop,noise,subspace[b]);
 
@@ -144,7 +144,8 @@ namespace Grid {
 	}
 
 	hermop.Op(noise,Mn); std::cout << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
-	subspace[b] = noise;
+	subspace[b]   = noise;
+
       }
 
       Orthogonalise();

From a70093361138b4fc123a88710c9e864844bde4c8 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:48:57 +0900
Subject: [PATCH 04/10] More info

---
 lib/algorithms/approx/Chebyshev.h | 46 ++++++++++++++++---------------
 1 file changed, 24 insertions(+), 22 deletions(-)

diff --git a/lib/algorithms/approx/Chebyshev.h b/lib/algorithms/approx/Chebyshev.h
index 0dad7179..f629ab08 100644
--- a/lib/algorithms/approx/Chebyshev.h
+++ b/lib/algorithms/approx/Chebyshev.h
@@ -50,6 +50,15 @@ namespace Grid {
       return;
     }
 
+    // Convenience for plotting the approximation
+    void   PlotApprox(std::ostream &out) {
+      out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
+      for(double x=lo;x<hi;x+=(hi-lo)/50.0){
+	out <<x<<"\t"<<approx(x)<<std::endl;
+      }
+    };
+
+
     Chebyshev(double _lo,double _hi,int _order, double (* func)(double) ){
       lo=_lo;
       hi=_hi;
@@ -95,46 +104,39 @@ namespace Grid {
       return sum;
     };
 
-    // Convenience for plotting the approximation
-    void   PlotApprox(std::ostream &out) {
-      out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
-      for(double x=lo;x<hi;x+=(hi-lo)/50.0){
-	out <<x<<"\t"<<approx(x)<<std::endl;
-      }
-    };
-
-    // Implement the required interface; could require Lattice base class
+    // Implement the required interface
     void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
 
-      Field T0 = in;
-      Field T1 = T0; // Field T1(T0._grid); more efficient but hardwires Lattice class
-      Field T2 = T1;
+      GridBase *grid=in._grid;
+
+      int vol=grid->gSites();
+
+      Field T0(grid); T0 = in;  
+      Field T1(grid); 
+      Field T2(grid);
+      Field y(grid);
       
-      // use a pointer trick to eliminate copies
       Field *Tnm = &T0;
       Field *Tn  = &T1;
       Field *Tnp = &T2;
-      Field y   = in;
-  
+
+      std::cout << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
+      // Tn=T1 = (xscale M + mscale)in
       double xscale = 2.0/(hi-lo);
       double mscale = -(hi+lo)/(hi-lo);
-
-      // Tn=T1 = (xscale M + mscale)in
-      Linop.Op(T0,y);
-
+      Linop.HermOp(T0,y);
       T1=y*xscale+in*mscale;
 
       // sum = .5 c[0] T0 + c[1] T1
       out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
-
       for(int n=2;n<order;n++){
 	
-	Linop.Op(*Tn,y);
+	Linop.HermOp(*Tn,y);
 
 	y=xscale*y+mscale*(*Tn);
 
 	*Tnp=2.0*y-(*Tnm);
-	
+
 	out=out+Coeffs[n]* (*Tnp);
 
 	// Cycle pointers to avoid copies

From 000766938150bc808b22cd06ce2b29d87c943729 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:51:56 +0900
Subject: [PATCH 05/10] Printing change

---
 lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h b/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
index 148359a8..fee41fa8 100644
--- a/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
+++ b/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
@@ -137,7 +137,7 @@ namespace Grid {
 
 	cp = axpy_norm(r,-a,q[peri_k],r);  
 
-	std::cout<< " VPCG_step resid" <<sqrt(cp/rsq)<<std::endl; 
+	std::cout<< " VPGCR_step resid" <<sqrt(cp/rsq)<<std::endl; 
 	if((k==nstep-1)||(cp<rsq)){
 	  return cp;
 	}

From 487fde849695be86ba20d34881fcf94e882d0262 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:52:23 +0900
Subject: [PATCH 06/10] This file is being developed and will remain hacky
 until the new algorithm is complete

---
 tests/Test_dwf_hdcr.cc | 176 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 169 insertions(+), 7 deletions(-)

diff --git a/tests/Test_dwf_hdcr.cc b/tests/Test_dwf_hdcr.cc
index bafdb639..488f7277 100644
--- a/tests/Test_dwf_hdcr.cc
+++ b/tests/Test_dwf_hdcr.cc
@@ -6,6 +6,10 @@ using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+RealD InverseApproximation(RealD x){
+  return 1.0/x;
+}
+
 template<class Fobj,class CComplex,int nbasis, class Matrix>
 class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
 public:
@@ -33,6 +37,27 @@ public:
       _Matrix(FineMatrix)
   {
   }
+
+  void PowerMethod(const FineField &in) {
+
+    FineField p1(in._grid);
+    FineField p2(in._grid);
+
+    MdagMLinearOperator<Matrix,FineField>   fMdagMOp(_Matrix);
+
+    p1=in;
+    RealD absp2;
+    for(int i=0;i<20;i++){
+      RealD absp1=std::sqrt(norm2(p1));
+      fMdagMOp.HermOp(p1,p2);// this is the G5 herm bit      
+      //      _FineOperator.Op(p1,p2);// this is the G5 herm bit      
+      RealD absp2=std::sqrt(norm2(p2));
+      if(i%10==9)
+	std::cout << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl;
+      p1=p2*(1.0/std::sqrt(absp2));
+    }
+  }
+
 #if 0
   void operator()(const FineField &in, FineField & out) {
 
@@ -49,7 +74,7 @@ public:
     std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
 
     // Build some solvers
-    ConjugateGradient<FineField>    fCG(1.0e-1,1000);
+    ConjugateGradient<FineField>    fCG(1.0e-3,1000);
     ConjugateGradient<CoarseVector>  CG(1.0e-8,100000);
 
     ////////////////////////////////////////////////////////////////////////
@@ -164,6 +189,7 @@ public:
   }
 #endif
   // ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in  
+#if 0
   void operator()(const FineField &in, FineField & out) {
 
     CoarseVector Csrc(_CoarseOperator.Grid());
@@ -171,11 +197,11 @@ public:
     CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
 
     ConjugateGradient<CoarseVector>  CG(1.0e-10,100000);
-    ConjugateGradient<FineField>    fCG(1.0e-3,1000);
+    ConjugateGradient<FineField>    fCG(3.0e-2,1000);
 
     HermitianLinearOperator<CoarseOperator,CoarseVector>  HermOp(_CoarseOperator);
     MdagMLinearOperator<CoarseOperator,CoarseVector>     MdagMOp(_CoarseOperator);
-    MdagMLinearOperator<Matrix,FineField>               fMdagMOp(_Matrix);
+    ShiftedMdagMLinearOperator<Matrix,FineField>        fMdagMOp(_Matrix,0.1);
 
     FineField tmp(in._grid);
     FineField res(in._grid);
@@ -191,7 +217,7 @@ public:
     CG(MdagMOp,Ctmp,Csol);
     _Aggregates.PromoteFromSubspace(Csol,Qin);
 
-
+    //    Qin=0;
     _FineOperator.Op(Qin,tmp);// A Q in
     tmp = in - tmp;            // in - A Q in
 
@@ -206,6 +232,116 @@ public:
     std::cout<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
 
   }
+#endif
+
+  void SmootherTest (const FineField & in){
+    
+    FineField vec1(in._grid);
+    FineField vec2(in._grid);
+    RealD lo[3] = { 0.5, 1.0, 2.0};
+
+    //    MdagMLinearOperator<Matrix,FineField>        fMdagMOp(_Matrix);
+    ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.5);
+
+    RealD Ni,r;
+
+    Ni = norm2(in);
+
+    for(int ilo=0;ilo<4;ilo++){
+      for(int ord=10;ord<60;ord+=10){
+
+	_FineOperator.AdjOp(in,vec1);
+
+	Chebyshev<FineField> Cheby  (lo[ilo],70.0,ord,InverseApproximation);
+	Cheby(fMdagMOp,vec1,vec2);    // solves  MdagM = g5 M g5M
+
+	_FineOperator.Op(vec2,vec1);// this is the G5 herm bit
+	vec1  = in - vec1;   // tmp  = in - A Min
+	r=norm2(vec1);
+	std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
+
+      }
+    }
+  }
+
+  void operator()(const FineField &in, FineField & out) {
+
+    CoarseVector Csrc(_CoarseOperator.Grid());
+    CoarseVector Ctmp(_CoarseOperator.Grid());
+    CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
+
+    ConjugateGradient<CoarseVector>  CG(1.0e-5,100000);
+    ConjugateGradient<FineField>    fCG(3.0e-2,1000);
+
+    HermitianLinearOperator<CoarseOperator,CoarseVector>  HermOp(_CoarseOperator);
+    MdagMLinearOperator<CoarseOperator,CoarseVector>     MdagMOp(_CoarseOperator);
+    //    MdagMLinearOperator<Matrix,FineField>        fMdagMOp(_Matrix);
+    ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.5);
+
+    FineField vec1(in._grid);
+    FineField vec2(in._grid);
+
+    //    Chebyshev<FineField> Cheby    (0.5,70.0,30,InverseApproximation);
+    //    Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
+    Chebyshev<FineField> Cheby    (1.0,70.0,20,InverseApproximation);
+    Chebyshev<FineField> ChebyAccu(1.0,70.0,20,InverseApproximation);
+
+    _Aggregates.ProjectToSubspace  (Csrc,in);
+    _Aggregates.PromoteFromSubspace(Csrc,out);
+    std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
+    
+    //    ofstream fout("smoother");
+    //    Cheby.csv(fout);
+
+    // V11 multigrid.
+    // Use a fixed chebyshev and hope hermiticity helps.
+
+    // To make a working smoother for indefinite operator
+    // must multiply by "Mdag" (ouch loses all low mode content)
+    // and apply to poly approx of (mdagm)^-1.
+    // so that we end up with an odd polynomial.
+
+    RealD Ni = norm2(in);
+
+    _FineOperator.AdjOp(in,vec1);// this is the G5 herm bit
+    ChebyAccu(fMdagMOp,vec1,out);    // solves  MdagM = g5 M g5M
+
+    std::cout << "Smoother norm "<<norm2(out)<<std::endl;
+
+    // Update with residual for out
+    _FineOperator.Op(out,vec1);// this is the G5 herm bit
+    vec1  = in - vec1;   // tmp  = in - A Min
+
+    RealD r = norm2(vec1);
+
+    std::cout << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
+    
+    _Aggregates.ProjectToSubspace  (Csrc,vec1);
+    HermOp.AdjOp(Csrc,Ctmp);// Normal equations
+    CG(MdagMOp,Ctmp,Csol);
+    _Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s
+                                             // Q = Q[in - A Min]  
+    out = out+vec1;
+
+    // Three preconditioner smoothing -- hermitian if C3 = C1
+    // Recompute error
+    _FineOperator.Op(out,vec1);// this is the G5 herm bit
+    vec1  = in - vec1;   // tmp  = in - A Min
+    r=norm2(vec1);
+
+    std::cout << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
+
+    // Reapply smoother
+    _FineOperator.Op(vec1,vec2);  // this is the G5 herm bit
+    ChebyAccu(fMdagMOp,vec2,vec1);    // solves  MdagM = g5 M g5M
+
+    out =out+vec1;
+    _FineOperator.Op(out,vec1);// this is the G5 herm bit
+    vec1  = in - vec1;   // tmp  = in - A Min
+    r=norm2(vec1);
+    std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
+
+  }
 
 };
 
@@ -224,7 +360,7 @@ int main (int argc, char ** argv)
   ///////////////////////////////////////////////////
   // Construct a coarsened grid; utility for this?
   ///////////////////////////////////////////////////
-  const int block=4;
+  const int block=2;
   std::vector<int> clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/block;
@@ -265,7 +401,8 @@ int main (int argc, char ** argv)
   std::cout << "**************************************************"<< std::endl;
   DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
 
-  const int nbasis = 6;
+  const int nbasis = 32;
+  //  const int nbasis = 4;
 
   typedef Aggregation<vSpinColourVector,vTComplex,nbasis>              Subspace;
   typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>          CoarseOperator;
@@ -276,7 +413,19 @@ int main (int argc, char ** argv)
   std::cout << "**************************************************"<< std::endl;
   MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
   Subspace Aggregates(Coarse5d,FGrid);
-  Aggregates.CreateSubspace(RNG5,HermDefOp);
+  //  Aggregates.CreateSubspace(RNG5,HermDefOp,nbasis);
+  assert ( (nbasis & 0x1)==0);
+  int nb=nbasis/2;
+  std::cout << " nbasis/2 = "<<nb<<std::endl;
+  Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
+  for(int n=0;n<nb;n++){
+    G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
+    std::cout<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
+  }
+  for(int n=0;n<nbasis;n++){
+    std::cout << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n])  <<std::endl;
+  }
+
 //  for(int i=0;i<nbasis;i++){
 //    result =     Aggregates.subspace[i];
 //    Aggregates.subspace[i]=result+g5*result;
@@ -319,12 +468,18 @@ int main (int argc, char ** argv)
   MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermion> Precon(Aggregates, LDOp,HermIndefOp,Ddwf);
   TrivialPrecon<LatticeFermion> simple;
 
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Testing smoother efficacy"<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  Precon.SmootherTest(src);
+
   std::cout << "**************************************************"<< std::endl;
   std::cout << "Unprec CG "<< std::endl;
   std::cout << "**************************************************"<< std::endl;
   //  TrivialPrecon<LatticeFermion> simple;
   ConjugateGradient<LatticeFermion> fCG(1.0e-8,100000);
   fCG(HermDefOp,src,result);
+  //  exit(0);
 
   std::cout << "**************************************************"<< std::endl;
   std::cout << "Testing GCR on indef matrix "<< std::endl;
@@ -332,6 +487,13 @@ int main (int argc, char ** argv)
   //  PrecGeneralisedConjugateResidual<LatticeFermion> UPGCR(1.0e-8,100000,simple,8,128);
   //  UPGCR(HermIndefOp,src,result);
 
+  
+  /// Get themax eval
+  std::cout << "**************************************************"<< std::endl;
+  std::cout <<" Applying power method to find spectral range      "<<std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  Precon.PowerMethod(src);
+
   std::cout << "**************************************************"<< std::endl;
   std::cout << "Building a two level PGCR "<< std::endl;
   std::cout << "**************************************************"<< std::endl;

From df2aac01f4e99cf4316cb3f800462b50449cbc33 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:52:59 +0900
Subject: [PATCH 07/10] This was needed to compile on gcc

---
 lib/qcd/hmc/integrators/Integrator.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 5eb319b2..90707d51 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -69,10 +69,10 @@ namespace Grid{
 	LatticeColourMatrix Umu(U._grid);
 	LatticeColourMatrix Pmu(U._grid);
 	for (int mu = 0; mu < Nd; mu++){
-	  Umu=peekLorentz(U, mu);
-	  Pmu=peekLorentz(*P, mu);
+	  Umu=PeekIndex<LorentzIndex>(U, mu);
+	  Pmu=PeekIndex<LorentzIndex>(*P, mu);
 	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
-	  pokeLorentz(U, Umu, mu);
+	  PokeIndex<LorentzIndex>(U, Umu, mu);
 	}
 
       }

From 8d654a86de43bc680f2a80ef3bff1da85ecc677f Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:53:23 +0900
Subject: [PATCH 08/10] Small pretty layout change

---
 lib/qcd/action/fermion/g5HermitianLinop.h | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/lib/qcd/action/fermion/g5HermitianLinop.h b/lib/qcd/action/fermion/g5HermitianLinop.h
index e5b691d8..864af86d 100644
--- a/lib/qcd/action/fermion/g5HermitianLinop.h
+++ b/lib/qcd/action/fermion/g5HermitianLinop.h
@@ -1,7 +1,9 @@
 #ifndef G5_HERMITIAN_LINOP
 #define G5_HERMITIAN_LINOP
+
 namespace Grid {
   namespace QCD {
+
 ////////////////////////////////////////////////////////////////////
 // Wrap an already herm matrix
 ////////////////////////////////////////////////////////////////////

From 5ac625f7167ec14ec0fff81c317f0a103f769612 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:54:09 +0900
Subject: [PATCH 09/10] No changes shown on git diff

---
 lib/algorithms/LinearOperator.h | 41 +++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h
index a804bceb..31e086ba 100644
--- a/lib/algorithms/LinearOperator.h
+++ b/lib/algorithms/LinearOperator.h
@@ -71,6 +71,47 @@ namespace Grid {
       }
     };
 
+    ////////////////////////////////////////////////////////////////////
+    // Construct herm op and shift it for mgrid smoother
+    ////////////////////////////////////////////////////////////////////
+    template<class Matrix,class Field>
+    class ShiftedMdagMLinearOperator : public LinearOperatorBase<Field> {
+      Matrix &_Mat;
+      RealD _shift;
+    public:
+    ShiftedMdagMLinearOperator(Matrix &Mat,RealD shift): _Mat(Mat), _shift(shift){};
+      // Support for coarsening to a multigrid
+      void OpDiag (const Field &in, Field &out) {
+	_Mat.Mdiag(in,out);
+	assert(0);
+      }
+      void OpDir  (const Field &in, Field &out,int dir,int disp) {
+	_Mat.Mdir(in,out,dir,disp);
+	assert(0);
+      }
+      void Op     (const Field &in, Field &out){
+	_Mat.M(in,out);
+	assert(0);
+      }
+      void AdjOp     (const Field &in, Field &out){
+	_Mat.Mdag(in,out);
+	assert(0);
+      }
+      void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+	_Mat.MdagM(in,out,n1,n2);
+	out = out + _shift*in;
+
+	ComplexD dot;	
+	dot= innerProduct(in,out);
+	n1=real(dot);
+	n2=norm2(out);
+      }
+      void HermOp(const Field &in, Field &out){
+	RealD n1,n2;
+	HermOpAndNorm(in,out,n1,n2);
+      }
+    };
+
     ////////////////////////////////////////////////////////////////////
     // Wrap an already herm matrix
     ////////////////////////////////////////////////////////////////////

From 4e94ddad46ba32509ba5f5821f733440cc2a1565 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Tue, 21 Jul 2015 13:56:22 +0900
Subject: [PATCH 10/10] Merge

---
 lib/algorithms/approx/Remez.h             |   2 +-
 lib/qcd/action/fermion/WilsonKernels.cc   |   2 -
 lib/qcd/hmc/integrators/Integrator.cc     |   2 +-
 lib/qcd/hmc/integrators/Integrator_base.h | 212 ++++++++++++++++++++++
 tests/Test_multishift_sqrt.cc             |  35 +++-
 5 files changed, 247 insertions(+), 6 deletions(-)
 create mode 100644 lib/qcd/hmc/integrators/Integrator_base.h

diff --git a/lib/algorithms/approx/Remez.h b/lib/algorithms/approx/Remez.h
index e6a973df..336382d2 100644
--- a/lib/algorithms/approx/Remez.h
+++ b/lib/algorithms/approx/Remez.h
@@ -15,7 +15,7 @@
 #ifndef INCLUDED_ALG_REMEZ_H
 #define INCLUDED_ALG_REMEZ_H
 
-#include <algorithms/approx/bigfloat_double.h>
+#include <algorithms/approx/bigfloat.h>
 
 #define JMAX 10000 //Maximum number of iterations of Newton's approximation
 #define SUM_MAX 10 // Maximum number of terms in exponential
diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc
index 4054280a..bdc18eee 100644
--- a/lib/qcd/action/fermion/WilsonKernels.cc
+++ b/lib/qcd/action/fermion/WilsonKernels.cc
@@ -1,5 +1,4 @@
 #include <Grid.h>
-
 namespace Grid {
 namespace QCD {
 
@@ -13,7 +12,6 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     vHalfSpinColourVector Uchi;
     int offset,local,perm, ptype;
 
-
     // Xp
     int ss = sF;
     offset = st._offsets [Xp][ss];
diff --git a/lib/qcd/hmc/integrators/Integrator.cc b/lib/qcd/hmc/integrators/Integrator.cc
index 8947405e..16a6513e 100644
--- a/lib/qcd/hmc/integrators/Integrator.cc
+++ b/lib/qcd/hmc/integrators/Integrator.cc
@@ -18,7 +18,7 @@ namespace Grid{
       Pmu = zero;
       for(int mu=0;mu<Nd;mu++){
 	SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
-	pokeLorentz(P, Pmu, mu);
+	PokeIndex<LorentzIndex>(P, Pmu, mu);
       }
       
     }
diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h
new file mode 100644
index 00000000..38798d67
--- /dev/null
+++ b/lib/qcd/hmc/integrators/Integrator_base.h
@@ -0,0 +1,212 @@
+//--------------------------------------------------------------------
+/*! @file Integrator_base.h
+ * @brief Declaration of classes for the abstract Molecular Dynamics integrator
+ *
+ * @author Guido Cossu
+ */
+//--------------------------------------------------------------------
+
+#ifndef INTEGRATOR_INCLUDED
+#define INTEGRATOR_INCLUDED
+
+#include <memory>
+
+class Observer;
+
+
+/*! @brief Abstract base class for Molecular Dynamics management */
+
+namespace Grid{
+  namespace QCD{
+
+    typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same size as the rest of the code
+    typedef std::vector<ActPtr> ActionLevel;
+    typedef std::vector<ActionLevel> ActionSet;
+    typedef std::vector<Observer*> ObserverList;
+    
+    class LeapFrog;
+
+
+    struct IntegratorParameters{
+      int Nexp;
+      int MDsteps;  // number of outer steps
+      RealD trajL;  // trajectory length 
+      RealD stepsize;
+
+      IntegratorParameters(int Nexp_,
+			   int MDsteps_, 
+			   RealD trajL_):
+	Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){};
+    };
+
+
+    namespace MDutils{
+      void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&);
+      void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
+    }
+
+    
+    
+
+
+    template< class IntegratorPolicy >
+    class Integrator{
+    private:
+      IntegratorParameters Params;
+      const ActionSet as;
+      const std::vector<int> Nrel; //relative step size per level
+      //ObserverList observers; // not yet
+      std::unique_ptr<LatticeLorentzColourMatrix> P;
+      
+      IntegratorPolicy TheIntegrator;// contains parameters too
+
+
+      void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
+	for(int a=0; a<as[level].size(); ++a){
+	  LatticeLorentzColourMatrix force(U._grid);
+	  as[level].at(a)->deriv(U,force);
+	  *P -= force*ep;
+	}
+      }
+
+
+      void update_U(LatticeLorentzColourMatrix&U, double ep){
+	//rewrite exponential to deal with the lorentz index?
+	LatticeColourMatrix Umu(U._grid);
+	LatticeColourMatrix Pmu(U._grid);
+	for (int mu = 0; mu < Nd; mu++){
+	  Umu=PeekIndex<LorentzIndex>(U, mu);
+	  Pmu=PeekIndex<LorentzIndex>(*P, mu);
+	  Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
+	}
+
+      }
+      
+      void register_observers();
+      void notify_observers();
+      
+      friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, 
+				     int level, std::vector<int>& clock,
+				     Integrator<LeapFrog>* Integ);
+    public:
+      Integrator(IntegratorParameters Par,
+		 ActionSet& Aset, std::vector<int> Nrel_):
+	Params(Par),as(Aset),Nrel(Nrel_){
+	assert(as.size() == Nrel.size());
+      };
+      
+      ~Integrator(){}
+
+
+      //Initialization of momenta and actions
+      void init(LatticeLorentzColourMatrix& U,
+		GridParallelRNG& pRNG){
+	std::cout<< "Integrator init\n";
+	if (!P)
+	  P = new LatticeLorentzColourMatrix(U._grid);
+	MDutils::generate_momenta(*P,pRNG);
+	
+	for(int level=0; level< as.size(); ++level){
+	  for(int actionID=0; actionID<as.at(level).size(); ++actionID){
+	    as[level].at(actionID)->init(U, pRNG);
+	  }
+	}
+      }
+
+      
+
+      RealD S(LatticeLorentzColourMatrix& U){
+	// Momenta
+	LatticeComplex Hloc = - trace((*P)*adj(*P));
+	Complex Hsum = sum(Hloc);
+	
+	RealD H = Hsum.real();
+
+	// Actions
+	for(int level=0; level<as.size(); ++level)
+	  for(int actionID=0; actionID<as.at(level).size(); ++actionID)
+	    H += as[level].at(actionID)->S(U);
+	
+	return H;
+      }
+      
+      
+
+      void integrate(LatticeLorentzColourMatrix& U, int level){
+	std::vector<int> clock;
+	clock.resize(as.size(),0);
+	for(int step=0; step< Params.MDsteps; ++step)   // MD step
+	  TheIntegrator.step(U,0,clock, *(this));
+      }
+    };
+    
+
+    class MinimumNorm2{
+      const double lambda = 0.1931833275037836;
+    public:
+      void step (LatticeLorentzColourMatrix& U, int level, std::vector<int>& clock);
+      
+    };
+    
+    class LeapFrog{
+    public:
+      void step (LatticeLorentzColourMatrix& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<LeapFrog>* Integ){
+	// cl  : current level
+	// fl  : final level
+	// eps : current step size
+	
+	int fl = Integ->as.size() -1;
+	double eps = Integ->Params.stepsize;
+	
+	// Get current level step size
+	for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l];
+	
+	int fin = 1;
+	for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l];
+	fin = 2*Integ->Params.MDsteps*fin - 1;
+	
+	for(int e=0; e<Integ->Nrel[level]; ++e){
+	  
+	  if(clock[level] == 0){    // initial half step
+	    Integ->update_P(U, level,eps/2);
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	  if(level == fl){          // lowest level
+	    Integ->update_U(U, eps);
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call
+	    step(U, level+1,clock, Integ);
+	  }
+	  if(clock[level] == fin){  // final half step
+	    Integ->update_P(U, level,eps/2);
+	    
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }else{                  // bulk step
+	    Integ->update_P(U, level,eps);
+	    
+	    clock[level]+=2;
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
+	  }
+	}
+
+
+
+
+      }
+    };
+
+
+
+
+    
+  }
+}
+#endif//INTEGRATOR_INCLUDED
diff --git a/tests/Test_multishift_sqrt.cc b/tests/Test_multishift_sqrt.cc
index 3d69fdf1..4fff263f 100644
--- a/tests/Test_multishift_sqrt.cc
+++ b/tests/Test_multishift_sqrt.cc
@@ -17,8 +17,12 @@ public:
     pRNG.SeedFixedIntegers(seeds);
 
     random(pRNG,sqrtscale);
-    sqrtscale = sqrtscale * adj(sqrtscale);// force real pos def
-    scale = sqrtscale * sqrtscale; 
+    sqrtscale = real(sqrtscale)*3.0+0.5;// force real pos def
+    scale = sqrtscale *sqrtscale; //scale should be bounded by 12.25
+
+    // 
+    //    sqrtscale = 2.0;
+    //    scale = 4.0;
   }
   // Support for coarsening to a multigrid
   void OpDiag (const Field &in, Field &out) {};
@@ -48,8 +52,18 @@ public:
   void ApplySqrt(const Field &in, Field &out){
     out = sqrtscale * in;
   }
+  void ApplyInverse(const Field &in, Field &out){
+    out = pow(scale,-1.0) * in;
+  }
 };
 
+RealD InverseApproximation(RealD x){
+  return 1.0/x;
+}
+RealD SqrtApproximation(RealD x){
+  return std::sqrt(x);
+}
+
 
 int main (int argc, char ** argv)
 {
@@ -114,5 +128,22 @@ int main (int argc, char ** argv)
   error = summed - combined;
   std::cout << " summed-combined "<<norm2(error)    <<std::endl;
 
+
+  src=1.0;
+  Chebyshev<LatticeFermion> Cheby(0.1,40.0,200,InverseApproximation);
+
+  std::cout<<"Chebuy approx vector "<<std::endl;
+  Cheby(Diagonal,src,combined);
+  std::ofstream of("cheby");
+  Cheby.csv(of);
+
+  Diagonal.ApplyInverse(src,reference);
+  error = reference - combined;
+
+  std::cout << "Chebyshev inverse test "<<std::endl;
+  std::cout << " Reference "<<norm2(reference)<<std::endl;
+  std::cout << " combined  "<<norm2(combined) <<std::endl;
+  std::cout << " error     "<<norm2(error)    <<std::endl;
+
   Grid_finalize();
 }