diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 1deb4093..7a6ad665 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -411,7 +411,7 @@ class KrylovSchur { if (Nconv >= Nstop || i == MaxIter - 1) { std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " << i << "." << std::endl; - basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); + // basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); // Think this might have been the issue std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; // writeEigensystem(path); @@ -629,7 +629,7 @@ class KrylovSchur { } // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) || - // checkRitzEstimate(); + checkRitzEstimate(); return Nconv; } diff --git a/examples/Example_wilson_evecs.cc b/examples/Example_wilson_evecs.cc new file mode 100644 index 00000000..cc4c4a11 --- /dev/null +++ b/examples/Example_wilson_evecs.cc @@ -0,0 +1,383 @@ +/************************************************************************************* + + Script for studying the Wilson eigenvectors resulting from the Krylov-Schur process. + + 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: 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 +} + +template void readFile(T& out, 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 << "Reads at: " << fname << std::endl; + Grid::emptyUserRecord record; + // Grid::ScidacReader SR(out.Grid()->IsBoss()); + Grid::ScidacReader SR; + SR.open(fname); + SR.readScidacFieldRecord(out, record); + SR.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 + int Nevecs = Nk; // don't write all of them + std::vector evecs = KS.getEvecs(); + for (int i = 0; i < Nevecs; 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 "< lat_size {16, 16, 16, 32}; + std::vector lat_size {32, 32, 32, 32}; + 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 = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + // GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,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); + + std::cout << GridLogMessage << "Loaded configuration" << std::endl; + + // RealD mass = 0.01; + RealD M5 = 1.8; + + // Wilson mass + RealD mass = -1.6; + + std::cout << GridLogMessage << "masses specified" << std::endl; + + std::vector boundary = {1,1,1,-1}; + WilsonFermionD::ImplParams Params(boundary); + + // DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + // NonHermitianLinearOperator DLinOp (Ddwf); + + // WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass); + WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params); + NonHermitianLinearOperator DLinOp (Dwilson); + + std::cout << GridLogMessage << "Dirac operator defined" << std::endl; + + std::string eigenPath = "/home/poare/lqcd/multigrid/spectra/32cube-rho0.124-tau4/U_smr_3.000000/Nm72_Nk24_8111835.aurora-pbs-0001.hostmgmt.cm.aurora.alcf.anl.gov/"; + + std::cout << GridLogMessage << "Loading eigenvalues" << std::endl; + std::ifstream evalFile(eigenPath + "evals.txt"); + std::string str; + std::vector evals; + while (std::getline(evalFile, str)) { + std::cout << GridLogMessage << "Reading line: " << str << std::endl; + int i1 = str.find("(") + 1; + int i2 = str.find(",") + 1; + int i3 = str.find(")"); + std::cout << "i1,i2,i3 = " << i1 << "," << i2 << "," << i3 << std::endl; + std::string reStr = str.substr(i1, i2 - i1); + std::string imStr = str.substr(i2, i3 - i2); + std::cout << GridLogMessage << "Parsed re = " << reStr << " and im = " << imStr << std::endl; + // ComplexD z (std::stof(reStr), std::stof(imStr)); + ComplexD z (std::stod(reStr), std::stod(imStr)); + evals.push_back(z); + } + std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; + + int Nevecs = 24; + std::vector evecs; + LatticeFermion evec (FGrid); + for (int i = 0; i < Nevecs; i++) { + std::string evecPath = eigenPath + "evec" + std::to_string(i); + readFile(evec, evecPath); + evecs.push_back(evec); + } + std::cout << GridLogMessage << "Evecs loaded" << std::endl; + + // Compute < evec | D - \lambda | evec > + std::cout << GridLogMessage << "Testing eigenvectors" << std::endl; + LatticeFermion Devec (FGrid); + ComplexD ritz; + for (int i = 0; i < Nevecs; i++) { + Devec = Zero(); + DLinOp.Op(evecs[i], Devec); + ritz = std::sqrt(norm2(Devec - evals[i] * evecs[i])); + std::cout << GridLogMessage << "i = " << i << ", || (D - lambda) |vi> || = " << ritz << std::endl; + } + // Eigen::MatrixXcd Dw_evecs; + // Dw_evecs = Eigen::MatrixXcd::Zero(Nevecs, Nevecs); + // for (int i = 0; i < Nevecs; i++) { + // Linop.Op(evecs[i], Devec); + // for (int j = 0; j < Nevecs; j++) { + + // } + // } + + std::cout< KS, std::string outDir) { fEval.close(); // Write evecs + int Nevecs = Nk; // don't write all of them std::vector evecs = KS.getEvecs(); - for (int i = 0; i < Nk; i++) { + for (int i = 0; i < Nevecs; i++) { std::string fName = outDir + "/evec" + std::to_string(i); writeFile(evecs[i], fName); // using method from Grid/HMC/ComputeWilsonFlow.cc } @@ -269,13 +270,16 @@ int main (int argc, char ** argv) 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; - } + // RitzFilter RF; + // if (argc == 8) { + // std::string rf = argv[7]; + // RF = selectRitzFilter(rf); + // } else { + // RF = EvalReSmall; + // } + // RitzFilter RF; + std::string rf = argv[7]; + RitzFilter RF = selectRitzFilter(rf); std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; const int Ls=16; @@ -288,8 +292,10 @@ int main (int argc, char ** argv) GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + // GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + // GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridCartesian * FGrid = UGrid; + GridRedBlackCartesian * FrbGrid = UrbGrid; std::vector seeds4({1,2,3,4}); GridParallelRNG RNG4(UGrid); @@ -303,21 +309,30 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << "Loaded configuration" << std::endl; - RealD mass = 0.01; + // RealD mass = 0.01; RealD M5 = 1.8; + // Wilson mass + RealD mass = -1.6; + std::cout << GridLogMessage << "masses specified" << std::endl; - // Define domain wall D. This is giving a floating point exception? - DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); - NonHermitianLinearOperator DLinOp (Ddwf); + std::vector boundary = {1,1,1,-1}; + WilsonFermionD::ImplParams Params(boundary); - std::cout << GridLogMessage << "DWF operator defined" << std::endl; + // DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + // NonHermitianLinearOperator DLinOp (Ddwf); + + // WilsonFermionD Dwilson(Umu, *FGrid, *FrbGrid, mass); + WilsonFermionD Dwilson(Umu, *UGrid, *UrbGrid, mass, Params); + NonHermitianLinearOperator DLinOp (Dwilson); + + std::cout << GridLogMessage << "Dirac operator defined" << std::endl; // Define PV^dag D (if we want) - DomainWallFermionD Dpv(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, 1.0, M5); - typedef PVdagMLinearOperator PVdagM_t; - PVdagM_t PVdagM(Ddwf, Dpv); + // DomainWallFermionD Dpv(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, 1.0, M5); + // typedef PVdagMLinearOperator PVdagM_t; + // PVdagM_t PVdagM(Ddwf, Dpv); std::cout<