diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index a83e7014..1ce4e5d9 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -265,34 +265,36 @@ class KrylovSchur { RealD approxLambdaMax; RealD beta_k; Field u; // Residual vector perpendicular to Krylov space (u_{k+1} in notes) - Eigen::VectorXcd b; // b vector in Schur decomposition (e_{k+1} in Arnoldi). + Eigen::VectorXcd b; // b vector in Schur decomposition (e_{k+1} in Arnoldi). std::vector basis; // orthonormal Krylov basis - Eigen::MatrixXcd Rayleigh; // Rayleigh quotient of size Nbasis (after construction) - Eigen::MatrixXcd Qt; // Transpose of basis rotation which projects out high modes. + Eigen::MatrixXcd Rayleigh; // Rayleigh quotient of size Nbasis (after construction) + Eigen::MatrixXcd Qt; // Transpose of basis rotation which projects out high modes. - Eigen::VectorXcd evals; // evals of Rayleigh quotient - Eigen::MatrixXcd littleEvecs; // Nm x Nm evecs matrix + Eigen::VectorXcd evals; // evals of Rayleigh quotient + std::vector ritzEstimates; // corresponding ritz estimates for evals + Eigen::MatrixXcd littleEvecs; // Nm x Nm evecs matrix std::vector evecs; // Vector of evec fields - RitzFilter ritzFilter; // how to sort evals + RitzFilter ritzFilter; // how to sort evals public: KrylovSchur(LinearOperatorBase &_Linop, GridBase *_Grid, RealD _Tolerance, RitzFilter filter = EvalReSmall) : Linop(_Linop), Grid(_Grid), Tolerance(_Tolerance), ritzFilter(filter), u(_Grid), MaxIter(-1), Nm(-1), Nk(-1), Nstop (-1), - evals (0), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0) + evals (0), ritzEstimates (), evecs (), ssq (0.0), rtol (0.0), beta_k (0.0), approxLambdaMax (0.0) { u = Zero(); }; /* Getters */ - int getNk() { return Nk; } - Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } - Field getU() { return u; } - std::vector getBasis() { return basis; } - Eigen::VectorXcd getEvals() { return evals; } - std::vector getEvecs() { return evecs; } + int getNk() { return Nk; } + Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } + Field getU() { return u; } + std::vector getBasis() { return basis; } + Eigen::VectorXcd getEvals() { return evals; } + std::vector getRitzEstimates() { return ritzEstimates; } + std::vector getEvecs() { return evecs; } /** * Runs the Krylov-Schur loop. @@ -365,19 +367,13 @@ class KrylovSchur { std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl; std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl; - - // Rayleigh = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); + Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); Rayleigh = RayTmp; - // basis = std::vector (basis.begin(), basis.begin() + Nk); std::vector basisTmp = std::vector (basis.begin(), basis.begin() + Nk); basis = basisTmp; - // evecs = std::vector (evecs.begin(), evecs.begin() + Nk); - // littleEvecs = littleEvecs(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); - - // b = b.head(Nk); Eigen::VectorXcd btmp = b.head(Nk); b = btmp; @@ -580,9 +576,11 @@ class KrylovSchur { int Nconv = 0; std::cout << GridLogDebug << "b: " << b << std::endl; Field tmp (Grid); Field fullEvec (Grid); + ritzEstimates.clear(); for (int k = 0; k < evecs.size(); k++) { Eigen::VectorXcd evec_k = littleEvecs.col(k); RealD ritzEstimate = std::abs(b.dot(evec_k)); // b^\dagger s + ritzEstimates.push_back(ritzEstimate); std::cout << GridLogMessage << "Ritz estimate for evec " << k << " = " << ritzEstimate << std::endl; if (ritzEstimate < rtol) { Nconv++; diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 8ae1ddfd..64f98a49 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -1,12 +1,35 @@ /************************************************************************************* + Runs the Krylov-Schur algorithm on a (pre-conditioned) domain-wall fermion operator + to determine part of its spectrum. + + Usage : + $ ./Example_spec_kryschur + + Nm = Maximum size of approximation subspace. + Nk = Size of truncation subspace + maxiter = Maximum number of iterations. + Nstop = Stop when Nstop eigenvalues have converged. + inFile = Gauge configuration to read in. + outDir = Directory to write output to. + rf = (Optional) RitzFilter to sort with. Takes in any string in + {EvalNormSmall, EvalNormLarge, EvalReSmall, EvalReLarge, EvalImSmall, EvalImLarge} + + Output: + ${outDir}/evals.txt = Contains all eigenvalues. Each line is formatted as `$idx $eval $ritz`, where: + - $idx is the index of the eigenvalue. + - $eval is the eigenvalue, formated as "(re,im)". + - $ritz is the Ritz estimate of the eigenvalue (deviation from being a true eigenvalue) + ${outDir}/evec${idx} = Eigenvector $idx written out in SCIDAC format (if LIME is enabled). + Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_padded_cell.cc Copyright (C) 2023 -Author: Peter Boyle + Author: Peter Boyle + Author: Patrick Oare This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -26,10 +49,6 @@ Author: Peter Boyle *************************************************************************************/ /* END LEGAL */ -// copied here from Test_general_coarse_pvdagm.cc - -// copied here from Test_general_coarse_pvdagm.cc - #include #include @@ -76,9 +95,10 @@ void writeEigensystem(KrylovSchur KS, std::string outDir) { std::ofstream fEval; fEval.open(evalPath); Eigen::VectorXcd evals = KS.getEvals(); + std::vector ritz = KS.getRitzEstimates(); for (int i = 0; i < Nk; i++) { - // write Eigenvalues - fEval << i << " " << evals(i); + // write eigenvalues and Ritz estimates + fEval << i << " " << evals(i) << " " << ritz[i]; if (i < Nk - 1) { fEval << "\n"; } } fEval.close(); @@ -236,113 +256,11 @@ ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shif } }; -template -class MGPreconditioner : public LinearFunction< Lattice > { -public: - using LinearFunction >::operator(); - - typedef Aggregation Aggregates; - typedef typename Aggregation::FineField FineField; - typedef typename Aggregation::CoarseVector CoarseVector; - typedef typename Aggregation::CoarseMatrix CoarseMatrix; - typedef LinearOperatorBase FineOperator; - typedef LinearFunction FineSmoother; - typedef LinearOperatorBase CoarseOperator; - typedef LinearFunction CoarseSolver; - Aggregates & _Aggregates; - FineOperator & _FineOperator; - FineSmoother & _PreSmoother; - FineSmoother & _PostSmoother; - CoarseOperator & _CoarseOperator; - CoarseSolver & _CoarseSolve; - - int level; void Level(int lv) {level = lv; }; - - MGPreconditioner(Aggregates &Agg, - FineOperator &Fine, - FineSmoother &PreSmoother, - FineSmoother &PostSmoother, - CoarseOperator &CoarseOperator_, - CoarseSolver &CoarseSolve_) - : _Aggregates(Agg), - _FineOperator(Fine), - _PreSmoother(PreSmoother), - _PostSmoother(PostSmoother), - _CoarseOperator(CoarseOperator_), - _CoarseSolve(CoarseSolve_), - level(1) { } - - virtual void operator()(const FineField &in, FineField & out) - { - GridBase *CoarseGrid = _Aggregates.CoarseGrid; - // auto CoarseGrid = _CoarseOperator.Grid(); - CoarseVector Csrc(CoarseGrid); - CoarseVector Csol(CoarseGrid); - FineField vec1(in.Grid()); - FineField vec2(in.Grid()); - - std::cout< - // assert (argc == 5); + // Usage : $ ./Example_spec_kryschur std::string NmStr = argv[1]; std::string NkStr = argv[2]; std::string maxIterStr = argv[3]; @@ -374,34 +292,15 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Construct a coarsened grid - // poare TODO: replace this with the following line? - Coordinate clatt = lat_size; -// Coordinate clatt = GridDefaultLatt(); // [PO] initial line before I edited it - for(int d=0;d seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); - std::vector cseeds({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); LatticeFermion src(FGrid); random(RNG5,src); - LatticeFermion result(FGrid); result=Zero(); - LatticeFermion ref(FGrid); ref=Zero(); - LatticeFermion tmp(FGrid); - LatticeFermion err(FGrid); LatticeGaugeField Umu(UGrid); FieldMetaData header; - // std::string file ("/sdcc/u/poare/PETSc-Grid/ckpoint_EODWF_lat.125"); - // std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_EODWF_lat.125"); NerscIO::readConfiguration(Umu,header,file); // RealD mass=0.01; @@ -411,16 +310,6 @@ int main (int argc, char ** argv) DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5); - // const int nbasis = 20; // size of approximate basis for low-mode space - const int nbasis = 3; // size of approximate basis for low-mode space - const int cb = 0 ; - LatticeFermion prom(FGrid); - - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef LittleDiracOperator::CoarseVector CoarseVector; - - NextToNearestStencilGeometry5D geom(Coarse5d); - std::cout< Dsq (Ddwf); NonHermitianLinearOperator DLinOp (Ddwf); - // int Nm = 200; - // int Nk = 110; - // int maxIter = 2000; - // int Nstop = 100; - int Nm = std::stoi(NmStr); int Nk = std::stoi(NkStr); int maxIter = std::stoi(maxIterStr); @@ -446,8 +330,8 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << "Runnning Krylov Schur. Nm = " << Nm << ", Nk = " << Nk << ", maxIter = " << maxIter << ", Nstop = " << Nstop << std::endl; - // KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); // use preconditioned PV^\dag D_{dwf} - KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use D_{dwf} + KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); // use preconditioned PV^\dag D_{dwf} + // KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use D_{dwf} KrySchur(src, maxIter, Nm, Nk, Nstop); std::cout< + + Nm = Maximum size of approximation subspace. + Nk = Size of truncation subspace + maxiter = Maximum number of iterations. + Nstop = Stop when Nstop eigenvalues have converged. + inFile = Gauge configuration to read in. + outDir = Directory to write output to. + rf = (Optional) RitzFilter to sort with. Takes in any string in + {EvalNormSmall, EvalNormLarge, EvalReSmall, EvalReLarge, EvalImSmall, EvalImLarge} + + Output: + ${outDir}/evals.txt = Contains all eigenvalues. Each line is formatted as `$idx $eval $ritz`, where: + - $idx is the index of the eigenvalue. + - $eval is the eigenvalue, formated as "(re,im)". + - $ritz is the Ritz estimate of the eigenvalue (deviation from being a true eigenvalue) + ${outDir}/evec${idx} = Eigenvector $idx written out in SCIDAC format (if LIME is enabled). + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_padded_cell.cc + + Copyright (C) 2023 + + Author: Peter Boyle + Author: Patrick Oare + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +using namespace std; +using namespace Grid; + +template void writeFile(T& in, std::string const fname){ + #ifdef HAVE_LIME + // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 + std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in,record,0); // Lexico + WR.close(); + #endif +} + +/** + * Writes the eigensystem of a Krylov Schur object to a directory. + * + * Parameters + * ---------- + * std::string path + * Directory to write to. + */ +template +void writeEigensystem(KrylovSchur KS, std::string outDir) { + int Nk = KS.getNk(); + std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + + // Write evals + std::string evalPath = outDir + "/evals.txt"; + std::ofstream fEval; + fEval.open(evalPath); + Eigen::VectorXcd evals = KS.getEvals(); + std::vector ritz = KS.getRitzEstimates(); + for (int i = 0; i < Nk; i++) { + // write eigenvalues and Ritz estimates + fEval << i << " " << evals(i) << " " << ritz[i]; + if (i < Nk - 1) { fEval << "\n"; } + } + fEval.close(); + + // Write evecs + std::vector evecs = KS.getEvecs(); + for (int i = 0; i < Nk; i++) { + std::string fName = outDir + "/evec" + std::to_string(i); + writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc + } +} + +// Hermitize a DWF operator by squaring it +template +class SquaredLinearOperator : public LinearOperatorBase { + + public: + Matrix &_Mat; + + public: + SquaredLinearOperator(Matrix &Mat): _Mat(Mat) {}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + // std::cout << "Op is overloaded as HermOp" << std::endl; + HermOp(in, out); + } + void AdjOp (const Field &in, Field &out){ + HermOp(in, out); + } + void _Op (const Field &in, Field &out){ + // std::cout << "Op: M "< +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + RealD shift; +public: + ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< +class ShiftedComplexPVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; + ComplexD shift; +public: +ShiftedComplexPVdagMLinearOperator(ComplexD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; + + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + std::cout << "Op: PVdag M "< + std::string NmStr = argv[1]; + std::string NkStr = argv[2]; + std::string maxIterStr = argv[3]; + std::string NstopStr = argv[4]; + std::string file = argv[5]; + std::string outDir = argv[6]; + + RitzFilter RF; + if (argc == 8) { + std::string rf = argv[7]; + RF = selectRitzFilter(rf); + } else { + RF = EvalReSmall; + } + std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; + + //const int Ls=16; + const int Ls = 8; + +// GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + //std::vector lat_size {16, 16, 16, 32}; + std::vector lat_size {8, 8, 8, 8}; + std::cout << "Lattice size: " << lat_size << std::endl; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(lat_size, + GridDefaultSimd(Nd,vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = UGrid; + GridRedBlackCartesian* FrbGrid = UrbGrid; + + std::vector seeds4({1,2,3,4}); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + LatticeFermion src(FGrid); random(RNG4, src); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,file); + + RealD mass = -0.1; + WilsonFermionD DWilson (Umu,*FGrid,*FrbGrid,mass); + + std::cout< Dsq (DWilson); + NonHermitianLinearOperator DLinOp (DWilson); + + int Nm = std::stoi(NmStr); + int Nk = std::stoi(NkStr); + int maxIter = std::stoi(maxIterStr); + int Nstop = std::stoi(NstopStr); + + std::cout << GridLogMessage << "Runnning Krylov Schur. Nm = " << Nm << ", Nk = " << Nk << ", maxIter = " << maxIter + << ", Nstop = " << Nstop << std::endl; + + KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use DWilson + KrySchur(src, maxIter, Nm, Nk, Nstop); + + std::cout<