From 9dcd7ca7613899c3e0b1340b60fb0d04189af911 Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Mon, 8 Sep 2025 12:59:48 -0400 Subject: [PATCH 01/12] added IO for evecs / evals --- Grid/algorithms/iterative/KrylovSchur.h | 65 ++++++++++++++++++++++-- examples/Example_spec_kryschur.cc | 67 ++++++++++++++++++++++--- 2 files changed, 123 insertions(+), 9 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 0c660585..9906741b 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -32,6 +32,12 @@ See the full license in the file "LICENSE" in the top level distribution directo #ifndef GRID_KRYLOVSCHUR_H #define GRID_KRYLOVSCHUR_H +// #include +// #include +// #include + +// #include + NAMESPACE_BEGIN(Grid); /** @@ -167,6 +173,15 @@ class ComplexSchurDecomposition { }; +// template +// inline void writeFile(const Field &field, const std::string &fname) { +// emptyUserRecord record; +// ScidacWriter WR(field.Grid()->IsBoss()); +// WR.open(fname); +// WR.writeScidacFieldRecord(field, record, 0); // 0 = Lexico +// WR.close(); +// } + /** * Implementation of the Krylov-Schur algorithm. */ @@ -212,6 +227,7 @@ class KrylovSchur { /* Getters */ + int getNk() { return Nk; } Eigen::MatrixXcd getRayleighQuotient() { return Rayleigh; } Field getU() { return u; } std::vector getBasis() { return basis; } @@ -321,6 +337,9 @@ class KrylovSchur { << i << "." << std::endl; basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; + + // writeEigensystem(path); + return; } } @@ -413,8 +432,6 @@ class KrylovSchur { * the Rayleigh quotient. Assumes that the Rayleigh quotient has already been constructed (by * calling the operator() function). * - * TODO implement in parent class eventually. - * * Parameters * ---------- * Eigen::MatrixXcd& S @@ -653,7 +670,49 @@ class KrylovSchur { return; } + // void writeEvec(Field& 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 << GridLogMessage << "Writing evec 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 + // // What is the appropriate way to throw error? + // } + + // /** + // * Writes the eigensystem of a Krylov Schur object to a directory. + // * + // * Parameters + // * ---------- + // * std::string path + // * Directory to write to. + // */ + // void writeEigensystem(std::string outDir) { + // std::cout << GridLogMessage << "Writing output to directory: " << outDir << std::endl; + // // TODO write a scidac density file so that we can easily integrate with visualization toolkit + // std::string evalPath = outDir + "/evals.txt"; + // std::ofstream fEval; + // fEval.open(evalPath); + // for (int i = 0; i < Nk; i++) { + // // write Eigenvalues + // fEval << i << " " << evals(i); + // if (i < Nk - 1) { fEval << "\n"; } + // } + // fEval.close(); + + // 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 + // // writeEvec(evecs[i], fName); + // } + + // } + }; - + NAMESPACE_END(Grid); #endif diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 6b3160ca..61c6ab30 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -40,9 +40,57 @@ Author: Peter Boyle #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(); + for (int i = 0; i < Nk; i++) { + // write Eigenvalues + fEval << i << " " << evals(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 { @@ -293,12 +341,14 @@ int main (int argc, char ** argv) { Grid_init(&argc,&argv); - // Usage : $ ./Example_spec_kryschur + // Usage : $ ./Example_spec_kryschur // assert (argc == 5); - std::string NmStr = argv[1]; - std::string NkStr = argv[2]; + std::string NmStr = argv[1]; + std::string NkStr = argv[2]; std::string maxIterStr = argv[3]; - std::string NstopStr = argv[4]; + std::string NstopStr = argv[4]; + std::string file = argv[5]; + std::string outDir = argv[6]; //const int Ls=16; const int Ls = 8; @@ -342,10 +392,11 @@ int main (int argc, char ** argv) 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"); + // std::string file("/Users/patrickoare/libraries/PETSc-Grid/ckpoint_EODWF_lat.125"); NerscIO::readConfiguration(Umu,header,file); - RealD mass=0.01; + // RealD mass=0.01; + RealD mass=0.001; RealD M5=1.8; DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); @@ -397,6 +448,10 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << "Krylov Schur eigenvalues: " << std::endl << KrySchur.getEvals() << std::endl; //std::cout << GridLogMessage << "Lanczos eigenvalues: " << std::endl << levals << std::endl; + // KrySchur.writeEigensystem(outDir); + + writeEigensystem(KrySchur, outDir); + std::cout< Date: Tue, 9 Sep 2025 13:02:11 -0400 Subject: [PATCH 02/12] updated RitzFilter enum and the input to run krylov schur --- Grid/algorithms/Algorithms.h | 2 +- Grid/algorithms/iterative/Arnoldi.h | 33 ----------- Grid/algorithms/iterative/KrylovSchur.h | 76 ++++++++++++++++++++++--- examples/Example_spec_kryschur.cc | 11 +++- 4 files changed, 79 insertions(+), 43 deletions(-) diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 85c4767f..4c1b8d10 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -81,7 +81,7 @@ NAMESPACE_CHECK(PowerMethod); NAMESPACE_CHECK(multigrid); #include -#include #include +#include #endif diff --git a/Grid/algorithms/iterative/Arnoldi.h b/Grid/algorithms/iterative/Arnoldi.h index 19b5f6c7..6232681d 100644 --- a/Grid/algorithms/iterative/Arnoldi.h +++ b/Grid/algorithms/iterative/Arnoldi.h @@ -8,8 +8,6 @@ Copyright (C) 2015 Author: Peter Boyle Author: paboyle -Author: Chulwoo Jung -Author: Christoph Lehner Author: Patrick Oare This program is free software; you can redistribute it and/or modify @@ -34,37 +32,6 @@ See the full license in the file "LICENSE" in the top level distribution directo NAMESPACE_BEGIN(Grid); -/** - * Options for which Ritz values to keep in implicit restart. - */ -enum RitzFilter { - EvalNormSmall, // Keep evals with smallest norm - EvalNormLarge, // Keep evals with largest norm - EvalReSmall, // Keep evals with smallest real part - EvalReLarge // Keep evals with largest real part -}; - -// Select comparison function from RitzFilter -struct ComplexComparator -{ - RitzFilter f; - ComplexComparator (RitzFilter _f) : f(_f) {} - bool operator()(std::complex z1, std::complex z2) { - switch (f) { - case EvalNormSmall: - return std::abs(z1) < std::abs(z2); - case EvalNormLarge: - return std::abs(z1) > std::abs(z2); - case EvalReSmall: - return std::abs(std::real(z1)) < std::abs(std::real(z2)); - case EvalReLarge: - return std::abs(std::real(z1)) > std::abs(std::real(z2)); - default: - assert(0); - } - } -}; - /** * Implementation of the Arnoldi algorithm. */ diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 9906741b..a83e7014 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -8,8 +8,6 @@ Copyright (C) 2015 Author: Peter Boyle Author: paboyle -Author: Chulwoo Jung -Author: Christoph Lehner Author: Patrick Oare This program is free software; you can redistribute it and/or modify @@ -32,14 +30,76 @@ See the full license in the file "LICENSE" in the top level distribution directo #ifndef GRID_KRYLOVSCHUR_H #define GRID_KRYLOVSCHUR_H -// #include -// #include -// #include - -// #include - NAMESPACE_BEGIN(Grid); +/** + * Options for which Ritz values to keep in implicit restart. TODO move this and utilities into a new file + */ +enum RitzFilter { + EvalNormSmall, // Keep evals with smallest norm + EvalNormLarge, // Keep evals with largest norm + EvalReSmall, // Keep evals with smallest real part + EvalReLarge, // Keep evals with largest real part + EvalImSmall, // Keep evals with smallest imaginary part + EvalImLarge // Keep evals with largest imaginary part +}; + +/** Selects the RitzFilter corresponding to the input string. */ +RitzFilter selectRitzFilter(std::string s) { + if (s == "EvalNormSmall") { return EvalNormSmall; } else + if (s == "EvalNormLarge") { return EvalNormLarge; } else + if (s == "EvalReSmall") { return EvalReSmall; } else + if (s == "EvalReLarge") { return EvalReLarge; } else + if (s == "EvalImSmall") { return EvalImSmall; } else + if (s == "EvalImLarge") { return EvalImLarge; } else + { assert(0); } +} + +/** Returns a string saying which RitzFilter it is. */ +std::string rfToString(RitzFilter RF) { + switch (RF) { + case EvalNormSmall: + return "EvalNormSmall"; + case EvalNormLarge: + return "EvalNormLarge"; + case EvalReSmall: + return "EvalReSmall"; + case EvalReLarge: + return "EvalReLarge"; + case EvalImSmall: + return "EvalImSmall"; + case EvalImLarge: + return "EvalImLarge"; + default: + assert(0); + } +} + +// Select comparison function from RitzFilter +struct ComplexComparator +{ + RitzFilter RF; + ComplexComparator (RitzFilter _rf) : RF(_rf) {} + bool operator()(std::complex z1, std::complex z2) { + switch (RF) { + case EvalNormSmall: + return std::abs(z1) < std::abs(z2); + case EvalNormLarge: + return std::abs(z1) > std::abs(z2); + case EvalReSmall: + return std::abs(std::real(z1)) < std::abs(std::real(z2)); + case EvalReLarge: + return std::abs(std::real(z1)) > std::abs(std::real(z2)); + case EvalImSmall: + return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + case EvalImLarge: + return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); + default: + assert(0); + } + } +}; + /** * Computes a complex Schur decomposition of a complex matrix A using Eigen's matrix library. The Schur decomposition, * A = Q^\dag S Q diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 61c6ab30..aad64dbd 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -350,6 +350,15 @@ 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; + } + std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; + //const int Ls=16; const int Ls = 8; @@ -438,7 +447,7 @@ int main (int argc, char ** argv) << ", Nstop = " << Nstop << std::endl; // Arnoldi Arn(PVdagM, FGrid, 1e-8); // Arn(src, maxIter, Nm, Nk, Nstop); - KrylovSchur KrySchur (PVdagM, FGrid, 1e-8); + KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); KrySchur(src, maxIter, Nm, Nk, Nstop); std::cout< Date: Tue, 9 Sep 2025 15:14:11 -0400 Subject: [PATCH 03/12] added commented out line to run un-preconditioned DWF --- examples/Example_spec_kryschur.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index aad64dbd..8ae1ddfd 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -445,9 +445,9 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << "Runnning Krylov Schur. Nm = " << Nm << ", Nk = " << Nk << ", maxIter = " << maxIter << ", Nstop = " << Nstop << std::endl; - // Arnoldi Arn(PVdagM, FGrid, 1e-8); - // Arn(src, maxIter, Nm, Nk, Nstop); - KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); + + // 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< Date: Wed, 10 Sep 2025 15:41:00 -0400 Subject: [PATCH 04/12] added wilson spectrum example --- Grid/algorithms/iterative/KrylovSchur.h | 38 ++- examples/Example_spec_kryschur.cc | 179 +++---------- examples/Example_wilson_spectrum.cc | 340 ++++++++++++++++++++++++ 3 files changed, 388 insertions(+), 169 deletions(-) create mode 100644 examples/Example_wilson_spectrum.cc 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< Date: Wed, 10 Sep 2025 17:16:48 -0400 Subject: [PATCH 05/12] added inline to rf functions --- Grid/algorithms/iterative/KrylovSchur.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 1ce4e5d9..800bc826 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -45,7 +45,7 @@ enum RitzFilter { }; /** Selects the RitzFilter corresponding to the input string. */ -RitzFilter selectRitzFilter(std::string s) { +inline RitzFilter selectRitzFilter(std::string s) { if (s == "EvalNormSmall") { return EvalNormSmall; } else if (s == "EvalNormLarge") { return EvalNormLarge; } else if (s == "EvalReSmall") { return EvalReSmall; } else @@ -56,7 +56,7 @@ RitzFilter selectRitzFilter(std::string s) { } /** Returns a string saying which RitzFilter it is. */ -std::string rfToString(RitzFilter RF) { +inline std::string rfToString(RitzFilter RF) { switch (RF) { case EvalNormSmall: return "EvalNormSmall"; @@ -367,7 +367,7 @@ class KrylovSchur { std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl; std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl; - + Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); Rayleigh = RayTmp; From 0b92ef990c8f85243f9e347db70d60766c1ab15e Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Fri, 12 Sep 2025 13:31:39 -0400 Subject: [PATCH 06/12] found bug in unprec DWF: was using |\cdot| in comparison for the eigenvalue sorting --- Grid/algorithms/iterative/KrylovSchur.h | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 800bc826..c1de912a 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -86,14 +86,22 @@ struct ComplexComparator return std::abs(z1) < std::abs(z2); case EvalNormLarge: return std::abs(z1) > std::abs(z2); + // case EvalReSmall: + // return std::abs(std::real(z1)) < std::abs(std::real(z2)); // DELETE THE ABS HERE!!! + // case EvalReLarge: + // return std::abs(std::real(z1)) > std::abs(std::real(z2)); + // case EvalImSmall: + // return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + // case EvalImLarge: + // return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); case EvalReSmall: - return std::abs(std::real(z1)) < std::abs(std::real(z2)); + return std::real(z1) < std::real(z2); // DELETE THE ABS HERE!!! case EvalReLarge: - return std::abs(std::real(z1)) > std::abs(std::real(z2)); + return std::real(z1) > std::real(z2); case EvalImSmall: - return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + return std::imag(z1) < std::imag(z2); case EvalImLarge: - return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); + return std::imag(z1) > std::imag(z2); default: assert(0); } @@ -311,6 +319,7 @@ class KrylovSchur { ssq = norm2(v0); RealD approxLambdaMax = approxMaxEval(v0); rtol = Tolerance * approxLambdaMax; + std::cout << GridLogMessage << "Approximate max eigenvalue: " << approxLambdaMax << std::endl; // rtol = Tolerance; b = Eigen::VectorXcd::Zero(Nm); // start as e_{k+1} From 612049f1ddd53607531580440c8278e0577c7755 Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Thu, 18 Sep 2025 15:09:31 -0400 Subject: [PATCH 07/12] commented out evec writer because it was taking up all the space on SDCC --- examples/Example_spec_kryschur.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/Example_spec_kryschur.cc b/examples/Example_spec_kryschur.cc index 64f98a49..f87d32be 100644 --- a/examples/Example_spec_kryschur.cc +++ b/examples/Example_spec_kryschur.cc @@ -103,12 +103,12 @@ void writeEigensystem(KrylovSchur KS, std::string outDir) { } 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 - } + // Write evecs (TODO: very heavy on storage costs! Don't write them all out) + // 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 From fa30c791aac8f5cb1a23945fdcc403c270534a3b Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Tue, 23 Sep 2025 15:24:50 +0000 Subject: [PATCH 08/12] updated wilson spec --- examples/Example_wilson_spectrum.cc | 30 +++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/examples/Example_wilson_spectrum.cc b/examples/Example_wilson_spectrum.cc index 7397691b..6709cce7 100644 --- a/examples/Example_wilson_spectrum.cc +++ b/examples/Example_wilson_spectrum.cc @@ -2,6 +2,8 @@ Runs the Krylov-Schur algorithm on a Wilson fermion operator to determine part of its spectrum. + TODO rename this file: really is running the topology change jobs on Aurora. + Usage : $ ./Example_spec_kryschur @@ -276,12 +278,11 @@ int main (int argc, char ** argv) } std::cout << "Sorting eigenvalues using " << rfToString(RF) << std::endl; - //const int Ls=16; - const int Ls = 8; + const int Ls=16; // 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::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()), @@ -298,18 +299,26 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(UGrid); FieldMetaData header; - NerscIO::readConfiguration(Umu,header,file); + NerscIO::readConfiguration(Umu, header, file); - RealD mass = -0.1; - WilsonFermionD DWilson (Umu,*FGrid,*FrbGrid,mass); + RealD mass = 0.01; + RealD M5 = 1.8; + + // Define domain wall D + DomainWallFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5); + NonHermitianLinearOperator DLinOp (Ddwf); + + // 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); std::cout< Dsq (DWilson); - NonHermitianLinearOperator DLinOp (DWilson); + // SquaredLinearOperator Dsq (DWilson); + // NonHermitianLinearOperator DLinOp (DWilson); int Nm = std::stoi(NmStr); int Nk = std::stoi(NkStr); @@ -319,7 +328,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 (DLinOp, FGrid, 1e-8, RF); // use DWilson + // KrylovSchur KrySchur (PVdagM, FGrid, 1e-8, RF); // use PV^\dag M + KrylovSchur KrySchur (DLinOp, FGrid, 1e-8, RF); // use Ddwf KrySchur(src, maxIter, Nm, Nk, Nstop); std::cout< Date: Thu, 25 Sep 2025 15:36:42 -0400 Subject: [PATCH 09/12] small bug fix for wilson spectrum since we're actually running DWF --- examples/Example_wilson_spectrum.cc | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/examples/Example_wilson_spectrum.cc b/examples/Example_wilson_spectrum.cc index 6709cce7..bd356bc5 100644 --- a/examples/Example_wilson_spectrum.cc +++ b/examples/Example_wilson_spectrum.cc @@ -288,8 +288,8 @@ int main (int argc, char ** argv) GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - GridCartesian* FGrid = UGrid; - GridRedBlackCartesian* FrbGrid = UrbGrid; + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); std::vector seeds4({1,2,3,4}); GridParallelRNG RNG4(UGrid); @@ -301,13 +301,19 @@ int main (int argc, char ** argv) FieldMetaData header; NerscIO::readConfiguration(Umu, header, file); + std::cout << GridLogMessage << "Loaded configuration" << std::endl; + RealD mass = 0.01; RealD M5 = 1.8; - // Define domain wall D + 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::cout << GridLogMessage << "DWF 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; From 4042ebf1bf007d023c17646cae21e809ca6243bd Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Mon, 20 Oct 2025 19:01:53 +0000 Subject: [PATCH 10/12] added ImNorm to sort --- Grid/algorithms/iterative/KrylovSchur.h | 75 ++++++++++++++++++------- 1 file changed, 56 insertions(+), 19 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index c1de912a..1deb4093 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -41,17 +41,21 @@ enum RitzFilter { EvalReSmall, // Keep evals with smallest real part EvalReLarge, // Keep evals with largest real part EvalImSmall, // Keep evals with smallest imaginary part - EvalImLarge // Keep evals with largest imaginary part + EvalImLarge, // Keep evals with largest imaginary part + EvalImNormSmall, // Keep evals with smallest |imaginary| part + EvalImNormLarge, // Keep evals with largest |imaginary| part }; /** Selects the RitzFilter corresponding to the input string. */ inline RitzFilter selectRitzFilter(std::string s) { - if (s == "EvalNormSmall") { return EvalNormSmall; } else - if (s == "EvalNormLarge") { return EvalNormLarge; } else - if (s == "EvalReSmall") { return EvalReSmall; } else - if (s == "EvalReLarge") { return EvalReLarge; } else - if (s == "EvalImSmall") { return EvalImSmall; } else - if (s == "EvalImLarge") { return EvalImLarge; } else + if (s == "EvalNormSmall") { return EvalNormSmall; } else + if (s == "EvalNormLarge") { return EvalNormLarge; } else + if (s == "EvalReSmall") { return EvalReSmall; } else + if (s == "EvalReLarge") { return EvalReLarge; } else + if (s == "EvalImSmall") { return EvalImSmall; } else + if (s == "EvalImLarge") { return EvalImLarge; } else + if (s == "EvalImNormSmall") { return EvalImNormSmall; } else + if (s == "EvalImNormLarge") { return EvalImNormLarge; } else { assert(0); } } @@ -70,6 +74,10 @@ inline std::string rfToString(RitzFilter RF) { return "EvalImSmall"; case EvalImLarge: return "EvalImLarge"; + case EvalImNormSmall: + return "EvalImNormSmall"; + case EvalImNormLarge: + return "EvalImNormLarge"; default: assert(0); } @@ -86,14 +94,6 @@ struct ComplexComparator return std::abs(z1) < std::abs(z2); case EvalNormLarge: return std::abs(z1) > std::abs(z2); - // case EvalReSmall: - // return std::abs(std::real(z1)) < std::abs(std::real(z2)); // DELETE THE ABS HERE!!! - // case EvalReLarge: - // return std::abs(std::real(z1)) > std::abs(std::real(z2)); - // case EvalImSmall: - // return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); - // case EvalImLarge: - // return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); case EvalReSmall: return std::real(z1) < std::real(z2); // DELETE THE ABS HERE!!! case EvalReLarge: @@ -102,6 +102,10 @@ struct ComplexComparator return std::imag(z1) < std::imag(z2); case EvalImLarge: return std::imag(z1) > std::imag(z2); + case EvalImNormSmall: + return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + case EvalImNormLarge: + return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); default: assert(0); } @@ -325,6 +329,9 @@ class KrylovSchur { b = Eigen::VectorXcd::Zero(Nm); // start as e_{k+1} b(Nm-1) = 1.0; + // basis = new std::vector (Nm, Grid); + // evecs.reserve(); + int start = 0; Field startVec = v0; littleEvecs = Eigen::MatrixXcd::Zero(Nm, Nm); @@ -360,6 +367,10 @@ class KrylovSchur { // basisRotate(evecs, Q, 0, Nm, 0, Nm, Nm); std::vector basis2; + // basis2.reserve(Nm); + // for (int i = start; i < Nm; i++) { + // basis2.emplace_back(Grid); + // } constructUR(basis2, basis, Qt, Nm); basis = basis2; @@ -431,14 +442,25 @@ class KrylovSchur { ComplexD coeff; Field w (Grid); // A acting on last Krylov vector. + // basis.reserve(Nm); + // for (int i = start; i < Nm; i++) { + // basis.emplace_back(Grid); + // } + // basis.assign(Nm, Field(Grid)); + // basis.resize(Nm); + // for (int i = start; i < Nm; i++) { + // basis[i] = Field(Grid); + // } + if (start == 0) { // initialize everything that we need. RealD v0Norm = 1 / std::sqrt(ssq); basis.push_back(v0Norm * v0); // normalized source + // basis[0] = v0Norm * v0; // normalized source Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); u = Zero(); } else { - assert( start == basis.size() ); // should be starting at the end of basis (start = Nk) + // assert( start == basis.size() ); // should be starting at the end of basis (start = Nk) std::cout << GridLogMessage << "Resetting Rayleigh and b" << std::endl; Eigen::MatrixXcd RayleighCp = Rayleigh; Rayleigh = Eigen::MatrixXcd::Zero(Nm, Nm); @@ -447,12 +469,14 @@ class KrylovSchur { // append b^\dag to Rayleigh, add u to basis Rayleigh(Nk, Eigen::seqN(0, Nk)) = b.adjoint(); basis.push_back(u); + // basis[start] = u; // TODO make sure this is correct b = Eigen::VectorXcd::Zero(Nm); } // Construct next Arnoldi vector by normalizing w_i = Dv_i - \sum_j v_j h_{ji} for (int i = start; i < Nm; i++) { Linop.Op(basis.back(), w); + // Linop.Op(basis[i], w); for (int j = 0; j < basis.size(); j++) { coeff = innerProduct(basis[j], w); // coeff = h_{ij}. Note that since {vi} is ONB it's OK to subtract it off after. Rayleigh(j, i) = coeff; @@ -475,6 +499,7 @@ class KrylovSchur { basis.push_back( (1.0/coeff) * w ); + // basis[i+1] = (1.0/coeff) * w; } // after iterations, update u and beta_k = ||u|| before norm @@ -506,8 +531,11 @@ class KrylovSchur { { std::cout << GridLogMessage << "Computing eigenvalues." << std::endl; - evecs.clear(); evals = S.diagonal(); + int n = evals.size(); // should be regular Nm + + evecs.clear(); + // evecs.assign(n, Field(Grid)); // TODO: is there a faster way to get the eigenvectors of a triangular matrix? // Rayleigh.triangularView tri; @@ -520,7 +548,7 @@ class KrylovSchur { std::cout << GridLogDebug << "Little evecs: " << littleEvecs << std::endl; // Convert evecs to lattice fields - for (int k = 0; k < evals.size(); k++) { + for (int k = 0; k < n; k++) { Eigen::VectorXcd vec = littleEvecs.col(k); Field tmp (basis[0].Grid()); tmp = Zero(); @@ -528,6 +556,7 @@ class KrylovSchur { tmp = tmp + vec[j] * basis[j]; } evecs.push_back(tmp); + // evecs[k] = tmp; } } @@ -583,13 +612,16 @@ class KrylovSchur { */ int converged() { int Nconv = 0; + int _Nm = evecs.size(); std::cout << GridLogDebug << "b: " << b << std::endl; Field tmp (Grid); Field fullEvec (Grid); ritzEstimates.clear(); - for (int k = 0; k < evecs.size(); k++) { + // ritzEstimates.resize(_Nm); + for (int k = 0; k < _Nm; k++) { Eigen::VectorXcd evec_k = littleEvecs.col(k); RealD ritzEstimate = std::abs(b.dot(evec_k)); // b^\dagger s ritzEstimates.push_back(ritzEstimate); + // ritzEstimates[k] = ritzEstimate; std::cout << GridLogMessage << "Ritz estimate for evec " << k << " = " << ritzEstimate << std::endl; if (ritzEstimate < rtol) { Nconv++; @@ -704,7 +736,9 @@ class KrylovSchur { */ void constructUR(std::vector& UR, std::vector &U, Eigen::MatrixXcd& R, int N) { Field tmp (Grid); + UR.clear(); + // UR.resize(N); std::cout << GridLogDebug << "R to rotate by (should be Rayleigh): " << R << std::endl; @@ -717,6 +751,7 @@ class KrylovSchur { } std::cout << GridLogDebug << "rotated norm at i = " << i << " is: " << norm2(tmp) << std::endl; UR.push_back(tmp); + // UR[i] = tmp; } return; } @@ -727,12 +762,14 @@ class KrylovSchur { void constructRU(std::vector& RU, std::vector &U, Eigen::MatrixXcd& R, int N) { Field tmp (Grid); RU.clear(); + // RU.resize(N); for (int i = 0; i < N; i++) { tmp = Zero(); for (int j = 0; j < N; j++) { tmp = tmp + R(i, j) * U[j]; } RU.push_back(tmp); + // RU[i] = tmp; } return; } From bf2a715ef74e5633c9fce7d451148daafc7d440d Mon Sep 17 00:00:00 2001 From: Patrick Oare Date: Fri, 31 Oct 2025 15:31:46 +0000 Subject: [PATCH 11/12] bug in wilson eigenvectors: ritz estimates not equalling deviation from being an evec --- Grid/algorithms/iterative/KrylovSchur.h | 4 +- examples/Example_wilson_evecs.cc | 383 ++++++++++++++++++++++++ examples/Example_wilson_spectrum.cc | 54 ++-- 3 files changed, 421 insertions(+), 20 deletions(-) create mode 100644 examples/Example_wilson_evecs.cc 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< Date: Fri, 31 Oct 2025 11:47:29 -0400 Subject: [PATCH 12/12] commented some slow code out --- Grid/algorithms/iterative/KrylovSchur.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 7a6ad665..337eba63 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -646,7 +646,7 @@ class KrylovSchur { // rotate basis by Rayleigh to construct UR // std::vector rotated; - std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl; + // std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl; std::vector rotated = basis; constructUR(rotated, basis, Rayleigh, k); // manually rotate @@ -665,10 +665,10 @@ class KrylovSchur { delta = delta / norm2(tmp); // relative tolerance deltaSum += delta; - std::cout << GridLogDebug << "Iteration " << i << std::endl; - std::cout << GridLogDebug << "Du = " << norm2(tmp) << std::endl; - std::cout << GridLogDebug << "rotated = " << norm2(rotated[i]) << std::endl; - std::cout << GridLogDebug << "b[i] = " << b(i) << std::endl; + // std::cout << GridLogDebug << "Iteration " << i << std::endl; + // std::cout << GridLogDebug << "Du = " << norm2(tmp) << std::endl; + // std::cout << GridLogDebug << "rotated = " << norm2(rotated[i]) << std::endl; + // std::cout << GridLogDebug << "b[i] = " << b(i) << std::endl; std::cout << GridLogMessage << "Deviation in decomp, column " << i << ": " << delta << std::endl; } std::cout << GridLogMessage << "Squared sum of relative deviations in decomposition: " << deltaSum << std::endl;