From 7ee66bf4533309bc1be3783fd278a85bfbb1237d Mon Sep 17 00:00:00 2001
From: Michael Marshall <43034299+mmphys@users.noreply.github.com>
Date: Sun, 19 Sep 2021 19:45:20 +0100
Subject: [PATCH 1/5] Make sure H5NS has empty definition if HDF5 built without
 C++ namespace. Add comment in Hdf5IO.cc indicating likely source of error
 using H5NS, i.e. lack of --enable-cxx in hdf5 configure.

---
 Grid/serialisation/Hdf5IO.cc  | 2 +-
 Grid/serialisation/Hdf5IO.h   | 4 ----
 Grid/serialisation/Hdf5Type.h | 4 +++-
 3 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/Grid/serialisation/Hdf5IO.cc b/Grid/serialisation/Hdf5IO.cc
index d3c7061e..db78df8e 100644
--- a/Grid/serialisation/Hdf5IO.cc
+++ b/Grid/serialisation/Hdf5IO.cc
@@ -33,7 +33,7 @@
 
 using namespace Grid;
 #ifndef H5_NO_NAMESPACE
-using namespace H5NS;
+using namespace H5NS; // Compile error here? Try adding --enable-cxx to hdf5 configure
 #endif
 
 // Writer implementation ///////////////////////////////////////////////////////
diff --git a/Grid/serialisation/Hdf5IO.h b/Grid/serialisation/Hdf5IO.h
index 46cb07e1..ae5e740b 100644
--- a/Grid/serialisation/Hdf5IO.h
+++ b/Grid/serialisation/Hdf5IO.h
@@ -40,10 +40,6 @@
 #include <Grid/tensors/Tensors.h>
 #include "Hdf5Type.h"
 
-#ifndef H5_NO_NAMESPACE
-#define H5NS H5
-#endif
-
 // default thresold above which datasets are used instead of attributes
 #ifndef HDF5_DEF_DATASET_THRES
 #define HDF5_DEF_DATASET_THRES 6u
diff --git a/Grid/serialisation/Hdf5Type.h b/Grid/serialisation/Hdf5Type.h
index 64dda349..d8a0dd22 100644
--- a/Grid/serialisation/Hdf5Type.h
+++ b/Grid/serialisation/Hdf5Type.h
@@ -5,7 +5,9 @@
 #include <complex>
 #include <memory>
 
-#ifndef H5_NO_NAMESPACE
+#ifdef H5_NO_NAMESPACE
+#define H5NS
+#else
 #define H5NS H5
 #endif
 

From 67e08aa952539c92cd03b9449924bad5c96e7f36 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Thu, 23 Sep 2021 23:39:55 +0200
Subject: [PATCH 2/5] New file not run yet

---
 examples/Example_wall_wall_spectrum.cc | 404 +++++++++++++++++++++++++
 1 file changed, 404 insertions(+)
 create mode 100644 examples/Example_wall_wall_spectrum.cc

diff --git a/examples/Example_wall_wall_spectrum.cc b/examples/Example_wall_wall_spectrum.cc
new file mode 100644
index 00000000..0d70f351
--- /dev/null
+++ b/examples/Example_wall_wall_spectrum.cc
@@ -0,0 +1,404 @@
+/*
+ * Warning: This code illustrative only: not well tested, and not meant for production use
+ * without regression / tests being applied
+ */
+
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+typedef SpinColourMatrix Propagator;
+typedef SpinColourVector Fermion;
+
+template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GridBase *grid;
+  GaugeField U;
+  
+  CovariantLaplacianCshift(GaugeField &_U)    :
+    grid(_U.Grid()),
+    U(_U) {  };
+
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out=Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
+      out = out - Gimpl::CovShiftForward(Umu,mu,in);    
+      out = out - Gimpl::CovShiftBackward(Umu,mu,in);    
+      out = out + 2.0*in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+void MakePhase(Coordinate mom,LatticeComplex &phase)
+{
+  GridBase *grid = phase.Grid();
+  auto latt_size = grid->GlobalDimensions();
+  ComplexD ci(0.0,1.0);
+  phase=Zero();
+
+  LatticeComplex coor(phase.Grid());
+  for(int mu=0;mu<Nd;mu++){
+    RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
+    LatticeCoordinate(coor,mu);
+    phase = phase + (TwoPiL * mom[mu]) * coor;
+  }
+  phase = exp(phase*ci);
+}
+void PointSource(Coordinate &coor,LatticePropagator &source)
+{
+  //  Coordinate coor({0,0,0,0});
+  source=Zero();
+  SpinColourMatrix kronecker; kronecker=1.0;
+  pokeSite(kronecker,source,coor);
+}
+void GFWallSource(int tslice,LatticePropagator &source)
+{
+  GridBase *grid = source.Grid();
+  LatticeComplex one(grid); one = ComplexD(1.0,0.0);
+  LatticeComplex zz(grid); zz=Zero();
+  LatticeInteger t(grid);
+  LatticeCoordinate(t,Tdir);
+  one = where(t==Integer(tslice), one, zz);
+  source = 1.0;
+  source = source * one;
+}
+
+void Z2WallSource(GridParallelRNG &RNG,int tslice,LatticePropagator &source)
+{
+  GridBase *grid = source.Grid();
+  LatticeComplex noise(grid);
+  LatticeComplex zz(grid); zz=Zero();
+  LatticeInteger t(grid);
+
+  RealD nrm=1.0/sqrt(2);
+  bernoulli(RNG, noise); // 0,1 50:50
+
+  noise = (2.*noise - Complex(1,1))*nrm;
+
+  LatticeCoordinate(t,Tdir);
+  noise = where(t==Integer(tslice), noise, zz);
+
+  source = 1.0;
+  source = source*noise;
+  std::cout << " Z2 wall " << norm2(source) << std::endl;
+}
+void GaugeFix(LatticeGaugeField &U,LatticeGaugeField &Ufix)
+{
+  Real alpha=0.05;
+
+  Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(U);
+
+  std::cout << " Initial plaquette "<<plaq << std::endl;
+
+  LatticeColourMatrix   xform(U.Grid()); 
+  Ufix = U;
+  int orthog=Nd-1;
+  FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Ufix,xform,alpha,10000,1.0e-12, 1.0e-12,true,orthog);
+  
+  plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Ufix);
+
+  std::cout << " Final plaquette "<<plaq << std::endl;
+}
+template<class Field>
+void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared)
+{
+  typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
+  Laplacian_t Laplacian(U);
+
+  Integer Iterations = 40;
+  Real width = 2.0;
+  Real coeff = (width*width) / Real(4*Iterations);
+
+  Field tmp(U.Grid());
+  smeared=unsmeared;
+  //  chi = (1-p^2/2N)^N kronecker
+  for(int n = 0; n < Iterations; ++n) {
+    Laplacian.M(smeared,tmp);
+    smeared = smeared - coeff*tmp;
+    std::cout << " smear iter " << n<<" " <<norm2(smeared)<<std::endl;
+  }
+}
+void GaussianSource(Coordinate &site,LatticeGaugeField &U,LatticePropagator &source)
+{
+  LatticePropagator tmp(source.Grid());
+  PointSource(site,source);
+  std::cout << " GaussianSource Kronecker "<< norm2(source)<<std::endl;
+  tmp = source;
+  GaussianSmear(U,tmp,source);
+  std::cout << " GaussianSource Smeared "<< norm2(source)<<std::endl;
+}
+void GaussianWallSource(GridParallelRNG &RNG,int tslice,LatticeGaugeField &U,LatticePropagator &source)
+{
+  Z2WallSource(RNG,tslice,source);
+  auto tmp = source;
+  GaussianSmear(U,tmp,source);
+}
+void SequentialSource(int tslice,Coordinate &mom,LatticePropagator &spectator,LatticePropagator &source)
+{
+  assert(mom.size()==Nd);
+  assert(mom[Tdir] == 0);
+
+  GridBase * grid = spectator.Grid();
+
+  LatticeInteger ts(grid);
+  LatticeCoordinate(ts,Tdir);
+  source = Zero();
+  source = where(ts==Integer(tslice),spectator,source); // Stick in a slice of the spectator, zero everywhere else
+
+  LatticeComplex phase(grid);
+  MakePhase(mom,phase);
+
+  source = source *phase;
+}
+template<class Action>
+void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator)
+{
+  GridBase *UGrid = D.GaugeGrid();
+  GridBase *FGrid = D.FermionGrid();
+
+  LatticeFermion src4  (UGrid); 
+  LatticeFermion src5  (FGrid); 
+  LatticeFermion result5(FGrid);
+  LatticeFermion result4(UGrid);
+  
+  ConjugateGradient<LatticeFermion> CG(1.0e-8,100000);
+  SchurRedBlackDiagMooeeSolve<LatticeFermion> schur(CG);
+  ZeroGuesser<LatticeFermion> ZG; // Could be a DeflatedGuesser if have eigenvectors
+  for(int s=0;s<Nd;s++){
+    for(int c=0;c<Nc;c++){
+      PropToFerm<Action>(src4,source,s,c);
+
+      D.ImportPhysicalFermionSource(src4,src5);
+
+      result5=Zero();
+      schur(D,src5,result5,ZG);
+      std::cout<<GridLogMessage
+	       <<"spin "<<s<<" color "<<c
+	       <<" norm2(src5d) "   <<norm2(src5)
+               <<" norm2(result5d) "<<norm2(result5)<<std::endl;
+
+      D.ExportPhysicalFermionSolution(result5,result4);
+
+      FermToProp<Action>(propagator,result4,s,c);
+    }
+  }
+}
+
+class MesonFile: Serializable {
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFile, std::vector<std::vector<Complex> >, data);
+};
+
+void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase)
+{
+  const int nchannel=4;
+  Gamma::Algebra Gammas[nchannel][2] = {
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::GammaTGamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::GammaTGamma5}
+  };
+
+  Gamma G5(Gamma::Algebra::Gamma5);
+
+  LatticeComplex meson_CF(q1.Grid());
+  MesonFile MF;
+
+  for(int ch=0;ch<nchannel;ch++){
+
+    Gamma Gsrc(Gammas[ch][0]);
+    Gamma Gsnk(Gammas[ch][1]);
+
+    meson_CF = trace(G5*adj(q1)*G5*Gsnk*q2*adj(Gsrc));
+
+    std::vector<TComplex> meson_T;
+    sliceSum(meson_CF,meson_T, Tdir);
+
+    int nt=meson_T.size();
+
+    std::vector<Complex> corr(nt);
+    for(int t=0;t<nt;t++){
+      corr[t] = TensorRemove(meson_T[t]); // Yes this is ugly, not figured a work around
+      std::cout << " channel "<<ch<<" t "<<t<<" " <<corr[t]<<std::endl;
+    }
+    MF.data.push_back(corr);
+  }
+
+  {
+    XmlWriter WR(file);
+    write(WR,"MesonFile",MF);
+  }
+}
+
+
+void WallSinkMesonTrace(std::string file,std::vector<Propagator> &q1,std::vector<Propagator> &q2)
+{
+  const int nchannel=4;
+  Gamma::Algebra Gammas[nchannel][2] = {
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::GammaTGamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::GammaTGamma5}
+  };
+
+  Gamma G5(Gamma::Algebra::Gamma5);
+  int nt=q1.size();
+  std::vector<Complex> meson_CF(nt);
+  MesonFile MF;
+
+  for(int ch=0;ch<nchannel;ch++){
+
+    Gamma Gsrc(Gammas[ch][0]);
+    Gamma Gsnk(Gammas[ch][1]);
+
+    std::vector<Complex> corr(nt);
+    for(int t=0;t<nt;t++){
+      meson_CF[t] = trace(G5*adj(q1[t])*G5*Gsnk*q2[t]*adj(Gsrc));
+      corr[t] = TensorRemove(meson_CF[t]); // Yes this is ugly, not figured a work around
+      std::cout << " channel "<<ch<<" t "<<t<<" " <<corr[t]<<std::endl;
+    }
+    MF.data.push_back(corr);
+  }
+
+  {
+    XmlWriter WR(file);
+    write(WR,"MesonFile",MF);
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  // Double precision grids
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  //////////////////////////////////////////////////////////////////////
+  // You can manage seeds however you like.
+  // Recommend SeedUniqueString.
+  //////////////////////////////////////////////////////////////////////
+  std::vector<int> seeds4({1,2,3,4}); 
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField Umu(UGrid);
+  LatticeGaugeField Ufixed(UGrid);
+  std::string config;
+  if( argc > 1 && argv[1][0] != '-' )
+  {
+    std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, argv[1]);
+    config=argv[1];
+  }
+  else
+  {
+    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
+    SU<Nc>::ColdConfiguration(Umu);
+    //    SU<Nc>::HotConfiguration(RNG4,Umu);
+    config="HotConfig";
+  }
+  GaugeFix(Umu,Ufixed);
+  Umu=Ufixed;
+    
+
+  std::vector<RealD> masses({ 0.004,0.02477,0.447} ); // u/d, s, c ??
+  std::vector<RealD> M5s   ({ 1.8,1.8,1.0} ); 
+  std::vector<RealD> bs   ({ 1.0,1.0,1.5} );  // DDM
+  std::vector<RealD> cs   ({ 0.0,0.0,0.5} );  // DDM
+  std::vector<int>   Ls_s ({ 16,16,12} );
+  std::vector<GridCartesian *> FGrids;
+  std::vector<GridRedBlackCartesian *> FrbGrids;
+
+  int nmass = masses.size();
+
+  std::vector<MobiusFermionR *> FermActs;
+  
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+  std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+
+  for(int m=0;m<masses.size();m++) {
+
+    RealD mass = masses[m];
+    RealD M5   = M5s[m];
+    RealD b    = bs[m];
+    RealD c    = cs[m];
+    int   Ls   = Ls_s[m];
+
+    FGrids.push_back(SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid));
+    FrbGrids.push_back(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid));
+
+    FermActs.push_back(new MobiusFermionR(Umu,*FGrids[m],*FrbGrids[m],*UGrid,*UrbGrid,mass,M5,b,c));
+  }
+
+  LatticePropagator point_source(UGrid);
+  LatticePropagator z2wall_source(UGrid);
+  LatticePropagator gfwall_source(UGrid);
+
+  Coordinate Origin({0,0,0,0});
+  PointSource   (Origin,point_source);
+  Z2WallSource  (RNG4,0,z2wall_source);
+  GFWallSource  (0,gfwall_source);
+  
+  std::vector<LatticePropagator> PointProps(nmass,UGrid);
+  std::vector<LatticePropagator> GaussProps(nmass,UGrid);
+  std::vector<LatticePropagator> Z2Props   (nmass,UGrid);
+  std::vector<LatticePropagator> GFProps   (nmass,UGrid);
+
+  for(int m=0;m<nmass;m++) {
+
+    Solve(*FermActs[m],z2wall_source    ,Z2Props[m]);
+    Solve(*FermActs[m],gfwall_source    ,GFProps[m]);
+  
+  }
+
+  LatticeComplex phase(UGrid);
+  Coordinate mom({0,0,0,0});
+  MakePhase(mom,phase);
+
+  std::vector<std::vector<Propagator> > wsnk_z2Props(nmass);
+  std::vector<std::vector<Propagator> > wsnk_gfProps(nmass);
+  for(int m=0;m<nmass;m++){
+    sliceSum(Z2Props[m],wsnk_z2Props[m],Tdir);
+    sliceSum(GFProps[m],wsnk_gfProps[m],Tdir);
+  }
+  
+  for(int m1=0 ;m1<nmass;m1++) {
+  for(int m2=m1;m2<nmass;m2++) {
+    std::stringstream ssg,ssz;
+    std::stringstream wssg,wssz;
+
+    /// Point sinks
+    ssg<<config<< "_m" << m1 << "_m"<< m2 << "p_gf_meson.xml";
+    ssz<<config<< "_m" << m1 << "_m"<< m2 << "p_z2_meson.xml";
+
+    MesonTrace(ssz.str(),Z2Props[m1],Z2Props[m2],phase);
+
+    /// Wall sinks
+    wssg<<config<< "_m" << m1 << "_m"<< m2 << "w_gf_meson.xml";
+    wssz<<config<< "_m" << m1 << "_m"<< m2 << "w_z2_meson.xml";
+
+    WallSinkMesonTrace(wssg.str(),wsnk_gfProps[m1],wsnk_gfProps[m2]);
+    WallSinkMesonTrace(wssz.str(),wsnk_z2Props[m1],wsnk_z2Props[m2]);
+    
+  }}
+
+  Grid_finalize();
+}
+
+
+

From 73a95fa96f09b9dd63d713d952630e3cba4de541 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 28 Sep 2021 12:44:26 +0100
Subject: [PATCH 3/5] LinearFunction loops over vectors by default, can be
 overloaded

---
 Grid/algorithms/LinearOperator.h | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h
index 1add212c..815226d2 100644
--- a/Grid/algorithms/LinearOperator.h
+++ b/Grid/algorithms/LinearOperator.h
@@ -530,6 +530,16 @@ public:
 template<class Field> class LinearFunction {
 public:
   virtual void operator() (const Field &in, Field &out) = 0;
+
+  virtual void operator() (const std::vector<Field> &in, std::vector<Field> &out)
+  {
+    assert(in.size() == out.size());
+
+    for (unsigned int i = 0; i < in.size(); ++i)
+    {
+      (*this)(in[i], out[i]);
+    }
+  }
 };
 
 template<class Field> class IdentityLinearFunction : public LinearFunction<Field> {

From 9523ad3d7345f443c272e6c148f7c54311b64e68 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Tue, 28 Sep 2021 12:45:47 +0100
Subject: [PATCH 4/5] vector version of Schur solver use vector guesser

---
 Grid/algorithms/iterative/SchurRedBlack.h | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h
index 15ef95c7..d97e4993 100644
--- a/Grid/algorithms/iterative/SchurRedBlack.h
+++ b/Grid/algorithms/iterative/SchurRedBlack.h
@@ -185,16 +185,19 @@ namespace Grid {
       ////////////////////////////////////////////////
       if ( subGuess ) guess_save.resize(nblock,grid);
 
-      for(int b=0;b<nblock;b++){
-        if(useSolnAsInitGuess) {
+      
+      if(useSolnAsInitGuess) {
+        for(int b=0;b<nblock;b++){
           pickCheckerboard(Odd, sol_o[b], out[b]);
-        } else {
-          guess(src_o[b],sol_o[b]); 
         }
+      } else {
+        guess(src_o, sol_o); 
+      }
 
-	if ( subGuess ) { 
-	  guess_save[b] = sol_o[b];
-	}
+	    if ( subGuess ) { 
+        for(int b=0;b<nblock;b++){
+          guess_save[b] = sol_o[b];
+        }
       }
       //////////////////////////////////////////////////////////////
       // Call the block solver

From 6c66b8d99780dbe2b75aa226dcfa8f5ac34e9e75 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 30 Sep 2021 19:25:12 +0100
Subject: [PATCH 5/5] deflated guesser can optionally be used with less vectors
 than provided

---
 Grid/algorithms/iterative/Deflation.h | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/Grid/algorithms/iterative/Deflation.h b/Grid/algorithms/iterative/Deflation.h
index 626b60e3..43fe3e35 100644
--- a/Grid/algorithms/iterative/Deflation.h
+++ b/Grid/algorithms/iterative/Deflation.h
@@ -54,15 +54,23 @@ class DeflatedGuesser: public LinearFunction<Field> {
 private:
   const std::vector<Field> &evec;
   const std::vector<RealD> &eval;
+  const unsigned int       N;
 
 public:
 
-  DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval) : evec(_evec), eval(_eval) {};
+  DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval)
+  : DeflatedGuesser(_evec, _eval, _evec.size())
+  {}
+
+  DeflatedGuesser(const std::vector<Field> & _evec, const std::vector<RealD> & _eval, const unsigned int _N)
+  : evec(_evec), eval(_eval), N(_N)
+  {
+    assert(evec.size()==eval.size());
+    assert(N <= evec.size());
+  } 
 
   virtual void operator()(const Field &src,Field &guess) {
     guess = Zero();
-    assert(evec.size()==eval.size());
-    auto N = evec.size();
     for (int i=0;i<N;i++) {
       const Field& tmp = evec[i];
       axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess);