Eigenvectors created. Still need to correctly set parameters for test.

Hadrons/Modules/MDistil/Distil.hpp

@@ -90,6 +90,73 @@ inline void SliceShare( GridBase * gridLowDim, GridBase * gridHighDim, void * Bu
 //#endif
   }
 }
+
+/*************************************************************************************
+ 
+ -Grad^2 (Peardon, 2009, pg 2, equation 3)
+ Field      Type of field the operator will be applied to
+ GaugeField Gauge field the operator will smear using
+ 
+ TODO CANDIDATE for integration into laplacian operator
+ should just require adding number of dimensions to act on to constructor,
+ where the default=all dimensions, but we could specify 3 spatial dimensions
+ 
+ *************************************************************************************/
+
+template<typename Field, typename GaugeField=LatticeGaugeFieldD>
+class LinOpPeardonNabla : public LinearOperatorBase<Field>, public LinearFunction<Field> {
+  typedef typename GaugeField::vector_type vCoeff_t;
+protected: // I don't really mind if _gf is messed with ... so make this public?
+  //GaugeField & _gf;
+  int          nd; // number of spatial dimensions
+  std::vector<Lattice<iColourMatrix<vCoeff_t> > > U;
+public:
+  // Construct this operator given a gauge field and the number of dimensions it should act on
+  LinOpPeardonNabla( GaugeField& gf, int dimSpatial = Grid::QCD::Tdir ) : /*_gf(gf),*/ nd{dimSpatial} {
+    assert(dimSpatial>=1);
+    for( int mu = 0 ; mu < nd ; mu++ )
+      U.push_back(PeekIndex<LorentzIndex>(gf,mu));
+      }
+  
+  // Apply this operator to "in", return result in "out"
+  void operator()(const Field& in, Field& out) {
+    assert( nd <= in._grid->Nd() );
+    conformable( in, out );
+    out = ( ( Real ) ( 2 * nd ) ) * in;
+    Field _tmp(in._grid);
+    typedef typename GaugeField::vector_type vCoeff_t;
+    //Lattice<iColourMatrix<vCoeff_t> > U(in._grid);
+    for( int mu = 0 ; mu < nd ; mu++ ) {
+      //U = PeekIndex<LorentzIndex>(_gf,mu);
+      out -= U[mu] * Cshift( in, mu, 1);
+      _tmp = adj( U[mu] ) * in;
+      out -= Cshift(_tmp,mu,-1);
+    }
+  }
+  
+  void OpDiag (const Field &in, Field &out) { assert(0); };
+  void OpDir  (const Field &in, Field &out,int dir,int disp) { assert(0); };
+  void Op     (const Field &in, Field &out) { assert(0); };
+  void AdjOp  (const Field &in, Field &out) { assert(0); };
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) { assert(0); };
+  void HermOp(const Field &in, Field &out) { operator()(in,out); };
+};
+
+template<typename Field>
+class LinOpPeardonNablaHerm : public LinearFunction<Field> {
+public:
+  OperatorFunction<Field>   & _poly;
+  LinearOperatorBase<Field> &_Linop;
+  
+  LinOpPeardonNablaHerm(OperatorFunction<Field> & poly,LinearOperatorBase<Field>& linop) : _poly(poly), _Linop(linop) {
+  }
+  
+  void operator()(const Field& in, Field& out) {
+    _poly(_Linop,in,out);
+  }
+};
+
+
 END_MODULE_NAMESPACE // Grid
 
 /******************************************************************************
@@ -100,6 +167,9 @@ BEGIN_HADRONS_NAMESPACE
 
 BEGIN_MODULE_NAMESPACE(MDistil)
 
+typedef Grid::Hadrons::EigenPack<LatticeColourVector>       DistilEP;
+typedef std::vector<std::vector<std::vector<SpinVector> > > DistilNoises;
+
 /******************************************************************************
  Make a lower dimensional grid
  ******************************************************************************/
@@ -386,6 +456,41 @@ public:
   }
 };
 
+/*************************************************************************************
+ 
+ Rotate eigenvectors into our phase convention
+ First component of first eigenvector is real and positive
+ 
+ *************************************************************************************/
+
+inline void RotateEigen(std::vector<LatticeColourVector> & evec)
+{
+  ColourVector cv0;
+  auto grid = evec[0]._grid;
+  std::vector<int> siteFirst(grid->Nd(),0);
+  peekSite(cv0, evec[0], siteFirst);
+  auto & cplx0 = cv0()()(0);
+  if( std::imag(cplx0) == 0 )
+    std::cout << GridLogMessage << "RotateEigen() : Site 0 : " << cplx0 << " => already meets phase convention" << std::endl;
+  else {
+    const auto cplx0_mag{std::abs(cplx0)};
+    const auto phase{std::conj(cplx0 / cplx0_mag)};
+    std::cout << GridLogMessage << "RotateEigen() : Site 0 : |" << cplx0 << "|=" << cplx0_mag << " => phase=" << (std::arg(phase) / 3.14159265) << " pi" << std::endl;
+    {
+      // TODO: Only really needed on the master slice
+      for( int k = 0 ; k < evec.size() ; k++ )
+        evec[k] *= phase;
+      if(grid->IsBoss()){
+        for( int c = 0 ; c < Nc ; c++ )
+          cv0()()(c) *= phase;
+        cplx0.imag(0); // This assumes phase convention is real, positive (so I get rid of rounding error)
+        //pokeSite(cv0, evec[0], siteFirst);
+        pokeLocalSite(cv0, evec[0], siteFirst);
+      }
+    }
+  }
+}
+
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE

Hadrons/Modules/MDistil/LapEvec.hpp

@@ -107,6 +107,7 @@ class LapEvecPar: Serializable
 public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(LapEvecPar,
                                   std::string,         gauge,
+                                  std::string,         EigenPackName,
                                   StoutParameters,     Stout,
                                   ChebyshevParameters, Cheby,
                                   LanczosParameters,   Lanczos,
@@ -125,7 +126,6 @@ class TLapEvec: public Module<LapEvecPar>
 {
 public:
   GAUGE_TYPE_ALIASES(FImpl,);
-  typedef std::vector<Grid::Hadrons::EigenPack<LatticeColourVector> > DistilEP;
 
 public:
   // constructor
@@ -155,6 +155,8 @@ MODULE_REGISTER_TMP(LapEvec, TLapEvec<FIMPL>, MDistil);
  TLapEvec implementation
  ******************************************************************************/
 
+//constexpr char szEigenPackSuffix[] = "_eigenPack";
+
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
 TLapEvec<FImpl>::TLapEvec(const std::string name) : gridLD{nullptr}, Module<LapEvecPar>(name)
@@ -182,7 +184,7 @@ std::vector<std::string> TLapEvec<FImpl>::getInput(void)
 template <typename FImpl>
 std::vector<std::string> TLapEvec<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {getName()}; // This is the higher dimensional eigenpack
     
     return out;
 }
@@ -195,16 +197,19 @@ void TLapEvec<FImpl>::setup(void)
   Environment & e{env()};
   gridHD = e.getGrid();
   gridLD = MakeLowerDimGrid( gridHD );
-  Nx = gridHD->_gdimensions[Xdir];
+  Nx = gridHD->_fdimensions[Xdir];
-  Ny = gridHD->_gdimensions[Ydir];
+  Ny = gridHD->_fdimensions[Ydir];
-  Nz = gridHD->_gdimensions[Zdir];
+  Nz = gridHD->_fdimensions[Zdir];
-  Nt = gridHD->_gdimensions[Tdir];
+  Nt = gridHD->_fdimensions[Tdir];
   // Temporaries
   envTmpLat(GaugeField, "Umu");
   envTmpLat(GaugeField, "Umu_stout");
   envTmpLat(GaugeField, "Umu_smear");
+  envTmp(LatticeGaugeField, "UmuNoTime",1,LatticeGaugeField(gridLD));
+  envTmp(LatticeColourVector, "src",1,LatticeColourVector(gridLD));
+  envTmp(std::vector<DistilEP>, "eig",1,std::vector<DistilEP>(Nt));
   // Output objects
-  envCreate(DistilEP, getName(), 1, Nt);
+  envCreate(DistilEP, getName(), 1, par().Lanczos.Nvec, gridHD );
 }
 
 // clean up any temporaries created by setup (that aren't stored in the environment)
@@ -215,14 +220,142 @@ void TLapEvec<FImpl>::Cleanup(void)
     delete gridLD;
     gridLD = nullptr;
   }
+  gridHD = nullptr;
 }
 
+/******************************************************************************
+ Calculate low-mode eigenvalues of the Laplacian
+ ******************************************************************************/
+
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TLapEvec<FImpl>::execute(void)
 {
   LOG(Message) << "execute() : start" << std::endl;
-  LOG(Message) << "Stout.steps=" << par().Stout.steps << ", Stout.parm=" << par().Stout.parm << std::endl;
+
+  // Alii for parameters
+  const int &TI{par().Distil.TI};
+  const int &LI{par().Distil.LI};
+  const int &nnoise{par().Distil.Nnoise};
+  const int &tsrc{par().Distil.tSrc};
+  const LanczosParameters &LPar{par().Lanczos};
+  const int &nvec{LPar.Nvec};
+  const bool exact_distillation{TI==Nt && LI==nvec};
+  const bool full_tdil{TI==Nt};
+  const int &Nt_inv{full_tdil ? 1 : TI};
+  const ChebyshevParameters &ChebPar{par().Cheby};
+
+  // Assertions on the parameters we read
+  assert(TI>1);
+  assert(LI>1);
+  if(exact_distillation)
+    assert(nnoise==1);
+  else
+    assert(nnoise>1);
+
+  // Stout smearing
+  envGetTmp(GaugeField, Umu);
+  envGetTmp(GaugeField, Umu_smear);
+  LOG(Message) << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
+  {
+    envGetTmp(GaugeField, Umu_stout);
+    const int &Steps{par().Stout.steps};
+    Smear_Stout<PeriodicGimplR> LS(par().Stout.parm);
+    for (int i = 0; i < Steps; i++) {
+      LS.smear(Umu_stout, Umu_smear);
+      Umu_smear = Umu_stout;
+    }
+  }
+  LOG(Message) << "Smeared plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu_smear) << std::endl;
+
+  // For debugging only, write logging output to a local file
+  std::ofstream * ll = nullptr;
+  const int rank{gridHD->ThisRank()};
+  if((0)) { // debug to a local log file
+    std::string filename{"Local_"};
+    filename.append(std::to_string(rank));
+    filename.append(".log");
+    ll = new std::ofstream(filename);
+  }
+
+  ////////////////////////////////////////////////////////////////////////
+  // Invert Peardon Nabla operator separately on each time-slice
+  ////////////////////////////////////////////////////////////////////////
+  
+  std::string sEigenPackName(par().EigenPackName);
+  bool bReturnValue = true;
+  auto & eig4d = envGet(DistilEP, getName() );
+  eig4d.resize(nvec,gridHD);
+  envGetTmp(std::vector<DistilEP>, eig);   // Eigenpack for each timeslice
+  envGetTmp(LatticeGaugeField, UmuNoTime); // Gauge field without time dimension
+  envGetTmp(LatticeColourVector, src);
+  const int Ntlocal{gridHD->LocalDimensions()[Tdir]};
+  const int Ntfirst{gridHD->LocalStarts()[Tdir]};
+  for(int t=Ntfirst;bReturnValue && t<Ntfirst+Ntlocal;t++){
+    std::cout << GridLogMessage << "------------------------------------------------------------" << std::endl;
+    std::cout << GridLogMessage << " Compute eigenpack, Timeslice  = " << t << std::endl;
+    std::cout << GridLogMessage << "------------------------------------------------------------" << std::endl;
+    
+    LOG(Message) << "eig.size()=" << eig.size() << std::endl;
+    eig[t].resize(LPar.Nk+LPar.Np,gridLD);
+    LOG(Message) << "After eig[t].resize" << std::endl;
+    
+    // Construct smearing operator
+    ExtractSliceLocal(UmuNoTime,Umu_smear,0,t-Ntfirst,Grid::QCD::Tdir); // switch to 3d/4d objects
+    LinOpPeardonNabla<LatticeColourVector> PeardonNabla(UmuNoTime);
+    std::cout << "Chebyshev preconditioning to order " << ChebPar.PolyOrder
+    << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl;
+    Chebyshev<LatticeColourVector> Cheb(ChebPar.alpha,ChebPar.beta,ChebPar.PolyOrder);
+    
+    //from Test_Cheby.cc
+    if ( ((0)) && Ntfirst == 0 && t==0) {
+      std::ofstream of("cheby_" + std::to_string(ChebPar.alpha) + "_" + std::to_string(ChebPar.beta) + "_" + std::to_string(ChebPar.PolyOrder));
+      Cheb.csv(of);
+    }
+
+    // Construct source vector according to Test_dwf_compressed_lanczos.cc
+    src=11.0;
+    RealD nn = norm2(src);
+    nn = Grid::sqrt(nn);
+    src = src * (1.0/nn);
+
+    GridLogIRL.Active(1);
+    LinOpPeardonNablaHerm<LatticeColourVector> PeardonNablaCheby(Cheb,PeardonNabla);
+    ImplicitlyRestartedLanczos<LatticeColourVector> IRL(PeardonNablaCheby,PeardonNabla,LPar.Nvec,LPar.Nk,LPar.Nk+LPar.Np,LPar.resid,LPar.MaxIt);
+    int Nconv = 0;
+    
+    if(ll) *ll << t << " : Before IRL.calc()" << std::endl;
+    IRL.calc(eig[t].eval,eig[t].evec,src,Nconv);
+    if(ll) *ll << t << " : After  IRL.calc()" << std::endl;
+    if( Nconv < LPar.Nvec ) {
+      bReturnValue = false;
+      if(ll) *ll << t << " : Convergence error : Only " << Nconv << " converged!" << std::endl;
+    } else {
+      if( Nconv > LPar.Nvec )
+        eig[t].resize( LPar.Nvec, gridLD );
+      std::cout << GridLogMessage << "Timeslice " << t << " has " << eig[t].eval.size() << " eigenvalues and " << eig[t].evec.size() << " eigenvectors." << std::endl;
+      
+      // Now rotate the eigenvectors into our phase convention
+      RotateEigen( eig[t].evec );
+      
+      // Write the eigenvectors and eigenvalues to disk
+      //std::cout << GridLogMessage << "Writing eigenvalues/vectors to " << pszEigenPack << std::endl;
+      eig[t].record.operatorXml="<OPERATOR>Michael</OPERATOR>";
+      eig[t].record.solverXml="<SOLVER>Felix</SOLVER>";
+      eig[t].write(sEigenPackName,false,t);
+      //std::cout << GridLogMessage << "Written eigenvectors" << std::endl;
+    }
+    for (int i=0;i<LPar.Nvec;i++){
+      std::cout << "Inserting Timeslice " << t << " into vector " << i << std::endl;
+      InsertSliceLocal(eig[t].evec[i],eig4d.evec[i],0,t,3);
+    }
+  }
+
+  // Close the local debugging log file
+  if( ll ) {
+    *ll << " Returning " << bReturnValue << std::endl;
+    delete ll;
+  }
   LOG(Message) << "execute() : end" << std::endl;
 }
 

tests/hadrons/Test_hadrons_distil.cc

@@ -235,11 +235,22 @@ void test_LapEvec(Application &application)
   // gauge field
   application.createModule<MGauge::Unit>(szGaugeName);
   // Now make an instance of the LapEvec object
-  MDistil::LapEvecPar par;
+  MDistil::LapEvecPar p;
-  par.Stout.steps = 173;
+  p.gauge = szGaugeName;
-  par.Stout.parm = -9.87654321;
+  p.EigenPackName = "ePack";
-  par.gauge = szGaugeName;
+  p.Distil.TI = 8;
-  application.createModule<MDistil::LapEvec>("LapEvec",par);
+  p.Distil.LI = 3;
+  p.Distil.Nnoise = 2;
+  p.Distil.tSrc = 0;
+  p.Stout.steps = 3;
+  p.Stout.parm = 0.2;
+  p.Cheby.PolyOrder = 11;
+  p.Cheby.alpha = 0.3;
+  p.Cheby.beta = 12.5;
+  p.Lanczos.Nvec = 5;
+  p.Lanczos.Nk = 6;
+  p.Lanczos.Np = 2;
+  application.createModule<MDistil::LapEvec>("LapEvec",p);
 }
 
 /////////////////////////////////////////////////////////////