From 00b0f75b0d9add932110ed2edaf5fc1770edc00b Mon Sep 17 00:00:00 2001 From: Michael Marshall Date: Thu, 24 Jan 2019 12:44:06 +0000 Subject: [PATCH] Eigenvectors created. Still need to correctly set parameters for test. --- Hadrons/Modules/MDistil/Distil.hpp | 105 +++++++++++++++++++ Hadrons/Modules/MDistil/LapEvec.hpp | 149 +++++++++++++++++++++++++-- tests/hadrons/Test_hadrons_distil.cc | 21 +++- 3 files changed, 262 insertions(+), 13 deletions(-) diff --git a/Hadrons/Modules/MDistil/Distil.hpp b/Hadrons/Modules/MDistil/Distil.hpp index 74e35085..ec502544 100644 --- a/Hadrons/Modules/MDistil/Distil.hpp +++ b/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 +class LinOpPeardonNabla : public LinearOperatorBase, public LinearFunction { + 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 > > 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(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 > U(in._grid); + for( int mu = 0 ; mu < nd ; mu++ ) { + //U = PeekIndex(_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 +class LinOpPeardonNablaHerm : public LinearFunction { +public: + OperatorFunction & _poly; + LinearOperatorBase &_Linop; + + LinOpPeardonNablaHerm(OperatorFunction & poly,LinearOperatorBase& 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 DistilEP; +typedef std::vector > > 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 & evec) +{ + ColourVector cv0; + auto grid = evec[0]._grid; + std::vector 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 diff --git a/Hadrons/Modules/MDistil/LapEvec.hpp b/Hadrons/Modules/MDistil/LapEvec.hpp index 910c4315..c06cc300 100644 --- a/Hadrons/Modules/MDistil/LapEvec.hpp +++ b/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 { public: GAUGE_TYPE_ALIASES(FImpl,); - typedef std::vector > DistilEP; public: // constructor @@ -155,6 +155,8 @@ MODULE_REGISTER_TMP(LapEvec, TLapEvec, MDistil); TLapEvec implementation ******************************************************************************/ +//constexpr char szEigenPackSuffix[] = "_eigenPack"; + // constructor ///////////////////////////////////////////////////////////////// template TLapEvec::TLapEvec(const std::string name) : gridLD{nullptr}, Module(name) @@ -182,7 +184,7 @@ std::vector TLapEvec::getInput(void) template std::vector TLapEvec::getOutput(void) { - std::vector out = {getName()}; + std::vector out = {getName()}; // This is the higher dimensional eigenpack return out; } @@ -195,16 +197,19 @@ void TLapEvec::setup(void) Environment & e{env()}; gridHD = e.getGrid(); gridLD = MakeLowerDimGrid( gridHD ); - Nx = gridHD->_gdimensions[Xdir]; - Ny = gridHD->_gdimensions[Ydir]; - Nz = gridHD->_gdimensions[Zdir]; - Nt = gridHD->_gdimensions[Tdir]; + Nx = gridHD->_fdimensions[Xdir]; + Ny = gridHD->_fdimensions[Ydir]; + Nz = gridHD->_fdimensions[Zdir]; + 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, "eig",1,std::vector(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::Cleanup(void) delete gridLD; gridLD = nullptr; } + gridHD = nullptr; } +/****************************************************************************** + Calculate low-mode eigenvalues of the Laplacian + ******************************************************************************/ + // execution /////////////////////////////////////////////////////////////////// template void TLapEvec::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::avgPlaquette(Umu) << std::endl; + { + envGetTmp(GaugeField, Umu_stout); + const int &Steps{par().Stout.steps}; + Smear_Stout 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::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, 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 PeardonNabla(UmuNoTime); + std::cout << "Chebyshev preconditioning to order " << ChebPar.PolyOrder + << " with parameters (alpha,beta) = (" << ChebPar.alpha << "," << ChebPar.beta << ")" << std::endl; + Chebyshev 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 PeardonNablaCheby(Cheb,PeardonNabla); + ImplicitlyRestartedLanczos 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="Michael"; + eig[t].record.solverXml="Felix"; + eig[t].write(sEigenPackName,false,t); + //std::cout << GridLogMessage << "Written eigenvectors" << std::endl; + } + for (int i=0;i(szGaugeName); // Now make an instance of the LapEvec object - MDistil::LapEvecPar par; - par.Stout.steps = 173; - par.Stout.parm = -9.87654321; - par.gauge = szGaugeName; - application.createModule("LapEvec",par); + MDistil::LapEvecPar p; + p.gauge = szGaugeName; + p.EigenPackName = "ePack"; + p.Distil.TI = 8; + 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("LapEvec",p); } /////////////////////////////////////////////////////////////