From d2189d59e0204cc5f02b4fac8b2750dd1f3bb2ef Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Tue, 16 Jun 2026 16:05:06 -0400 Subject: [PATCH] Adding parity flip and gamma 5 option for Block KS starting vector --- Grid/algorithms/Algorithms.h | 2 +- Grid/algorithms/iterative/BlockKrylovSchur.h | 67 +++++++++++++++++-- .../iterative/HarmonicBlockKrylovSchur.h | 50 +++++++++++++- examples/Example_krylov_schur.cc | 16 +++-- examples/LanParams.xml | 6 +- 5 files changed, 126 insertions(+), 15 deletions(-) diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index a6a0b794d..435d60bef 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -88,7 +88,7 @@ NAMESPACE_CHECK(multigrid); #include #include #include -//#include +#include #include //#include #include diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index 968942126..5a730cb2f 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -124,7 +124,17 @@ protected: public: std::vector evecs; - bool doEvalCheck = false; + bool doEvalCheck = false; + // When true (and Nblock even), only Nblock/2 starting vectors are required. + // The remaining Nblock/2 slots are filled by pairing each supplied vector + // with its parity-flipped partner (sign negated on odd checkerboard sites). + bool useParityFlip = false; + // Same pairing behaviour but the partner is produced by gamma5Func(v). + // gamma5Func must be set by the caller before operator() is called. + // It is a std::function so that BlockKrylovSchur.h does not need to know + // about the Grid QCD Gamma class (which is defined at a later include layer). + bool useGamma5 = false; + std::function gamma5Func; //-------------------------------------------------------------------- // Constructor @@ -161,10 +171,17 @@ public: Nstop = _Nstop; Nblock = _Nblock; - assert((int)v0.size() >= Nblock); + { + int divisor = 1; + if (useParityFlip) divisor *= 2; + if (useGamma5) divisor *= 2; + assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor); + std::cout << GridLogMessage << "BlockKrylovSchur: divisor = " <= Nstop || iter == MaxIter - 1) { std::cout << GridLogMessage << "BlockKrylovSchur: done after " << iter << " restarts, " << Nconv << " converged." << std::endl; - std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl; + std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; if (doEvalCheck) { Field w(Grid_); @@ -279,6 +323,21 @@ public: CVec getEvals() { return evals; } std::vector getRitzEstimates() { return ritzEstimates; } + // Negate the field on all odd-parity (checkerboard) sites. + // Uses pickCheckerboard/setCheckerboard so the SIMD layout is handled + // correctly and the negation runs in parallel. + static void parityFlippedField(const Field& v, Field& out) { + GridCartesian* cgrid = dynamic_cast(v.Grid()); + assert(cgrid != nullptr); + GridRedBlackCartesian rbgrid(cgrid); + Field veven(&rbgrid), vodd(&rbgrid); + pickCheckerboard(Even, veven, v); + pickCheckerboard(Odd, vodd, v); + vodd = -vodd; + setCheckerboard(out, veven); + setCheckerboard(out, vodd); + } + //-------------------------------------------------------------------- // Verification: print H and B, check A V = V H + F B^dag explicitly //-------------------------------------------------------------------- diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 9afeebff4..f22de77d2 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -126,7 +126,10 @@ class HarmonicBlockKrylovSchur { public: std::vector evecs; - bool doEvalCheck = false; + bool doEvalCheck = false; + bool useParityFlip = false; + bool useGamma5 = false; + std::function gamma5Func; //-------------------------------------------------------------------- // Constructor @@ -153,10 +156,16 @@ public: Nstop = _Nstop; Nblock = _Nblock; - assert((int)v0.size() >= Nblock); + { + int divisor = 1; + if (useParityFlip) divisor *= 2; + if (useGamma5) divisor *= 2; + assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor); + } assert(Nm % Nblock == 0); assert(Nk % Nblock == 0); assert(Nk < Nm); + if (useGamma5) assert(gamma5Func && "useGamma5: gamma5Func must be set"); int N = Nm; @@ -170,8 +179,31 @@ public: H = CMat::Zero(N, N); B = CMat::Zero(N, Nblock); + int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1); int start = 0; - std::vector startBlock(v0.begin(), v0.begin() + Nblock); + std::vector startBlock; + startBlock.reserve(Nblock); + for (int i = 0; i < Nblock / divisor; i++) { + std::vector group; + group.push_back(v0[i]); + if (useParityFlip) { + int n = (int)group.size(); + for (int j = 0; j < n; j++) { + Field fp(Grid_); + parityFlippedField(group[j], fp); + group.push_back(std::move(fp)); + } + } + if (useGamma5) { + int n = (int)group.size(); + for (int j = 0; j < n; j++) { + Field g5v(Grid_); + gamma5Func(group[j], g5v); + group.push_back(std::move(g5v)); + } + } + for (auto& f : group) startBlock.push_back(std::move(f)); + } for (int iter = 0; iter < MaxIter; iter++) { std::cout << GridLogMessage @@ -280,6 +312,18 @@ public: CVec getEvals() { return evals; } std::vector getRitzEstimates() { return ritzEstimates; } + static void parityFlippedField(const Field& v, Field& out) { + GridCartesian* cgrid = dynamic_cast(v.Grid()); + assert(cgrid != nullptr); + GridRedBlackCartesian rbgrid(cgrid); + Field veven(&rbgrid), vodd(&rbgrid); + pickCheckerboard(Even, veven, v); + pickCheckerboard(Odd, vodd, v); + vodd = -vodd; + setCheckerboard(out, veven); + setCheckerboard(out, vodd); + } + //-------------------------------------------------------------------- // Verification: check A V = V H + F B^dag explicitly //-------------------------------------------------------------------- diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 8cadbb105..397f61365 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -334,12 +334,13 @@ int main (int argc, char ** argv) // Run KrylovSchur and Arnoldi on a Hermitian matrix RealD shift=LanParams.shift; -#if 1 +#if 0 std::cout << GridLogMessage << "Running Krylov Schur" << std::endl; KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrySchur(src[0], maxIter, Nm, Nk, Nstop); KrySchur.doEvalCheck=true; - KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); +// KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); + KrySchur(src[0], maxIter, Nm, Nk, Nstop); std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #else std::cout << GridLogMessage << "Running BlockKrylovSchur" << std::endl; @@ -347,9 +348,16 @@ int main (int argc, char ** argv) Nblock=LanParams.Nblock; bool if_verify=false; if(LanParams.verify) if_verify=true; - BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); -// HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall); +// BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); + HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall); KrySchur.doEvalCheck=true; + KrySchur.useParityFlip=true; std::cout << GridLogMessage << "useParityFlip= " <