From d4b685cf2d005ba04983cbd110a408eade2a45f3 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Mon, 6 Apr 2026 12:00:33 -0400 Subject: [PATCH 1/2] Expanding debug output after restart in block KS --- Grid/algorithms/iterative/BlockKrylovSchur.h | 43 +++++++++++++++----- examples/Example_krylov_schur.cc | 10 ++--- 2 files changed, 37 insertions(+), 16 deletions(-) diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index fa86e241e..c3835c88d 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -314,23 +314,28 @@ public: return; } + // Full allocated dimension of H (= Nm); may be larger than nBasis after restart+truncation + int Nfull = (int)H.rows(); + std::cout << GridLogMessage << "======== BlockKrylovSchur::verify [" << label << "] ========" << std::endl; std::cout << GridLogMessage - << " nBasis = " << nBasis << " Nblock = " << Nblock - << " nF = " << nF << std::endl; + << " nBasis = " << nBasis << " Nfull = " << Nfull + << " Nblock = " << Nblock << " nF = " << nF << std::endl; - // ---- Print H ---- - std::cout << GridLogMessage << "H (" << nBasis << " x " << nBasis << "):" << std::endl; - for (int i = 0; i < nBasis; i++) { - for (int j = 0; j < nBasis; j++) + // ---- Print H (full Nm x Nm matrix) ---- + // After restart+truncation nBasis < Nfull; entries outside [0:nBasis,0:nBasis] + // should be ~0 — printing the full matrix makes this visible. + std::cout << GridLogMessage << "H (" << Nfull << " x " << Nfull << "):" << std::endl; + for (int i = 0; i < Nfull; i++) { + for (int j = 0; j < Nfull; j++) std::cout << " " << std::setw(14) << H(i, j); std::cout << std::endl; } - // ---- Print B ---- - std::cout << GridLogMessage << "B (" << nBasis << " x " << nF << "):" << std::endl; - for (int i = 0; i < nBasis; i++) { + // ---- Print B (full Nm x Nblock matrix) ---- + std::cout << GridLogMessage << "B (" << Nfull << " x " << nF << "):" << std::endl; + for (int i = 0; i < Nfull; i++) { for (int t = 0; t < nF; t++) std::cout << " " << std::setw(14) << B(i, t); std::cout << std::endl; @@ -352,13 +357,29 @@ public: std::cout << std::endl; } - // ---- max |H - M| ---- + // ---- max |H - M| (for the nBasis x nBasis Rayleigh-quotient block) ---- RealD maxHM = 0.0; for (int i = 0; i < nBasis; i++) for (int j = 0; j < nBasis; j++) maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j))); std::cout << GridLogMessage - << " max |H[i,j] - M[i,j]| = " << maxHM << std::endl; + << " max |H[i,j] - M[i,j]| (i,j < nBasis) = " << maxHM << std::endl; + + // ---- Check entries of H OUTSIDE the nBasis x nBasis block ---- + // These live in rows i >= nBasis or columns j >= nBasis. + // They must be ~0 after restart+truncation (before the coupling rows are + // re-filled by the next blockArnoldiIteration call). + // During a full Arnoldi run (nBasis == Nfull) this loop is a no-op. + RealD maxOutside = 0.0; + for (int i = 0; i < Nfull; i++) + for (int j = 0; j < Nfull; j++) + if (i >= nBasis || j >= nBasis) + maxOutside = std::max(maxOutside, std::abs(H(i,j))); + if (Nfull > nBasis) + std::cout << GridLogMessage + << " max |H[i,j]| outside [0:" << nBasis << ",0:" << nBasis + << "] block = " << maxOutside + << " (should be ~0 after restart+truncation)" << std::endl; // ---- Check orthonormality of basis ---- CMat G = CMat::Zero(nBasis, nBasis); diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index c945f852d..11f8a9fc4 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -347,7 +347,7 @@ int main (int argc, char ** argv) // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // Hacked, really EvalImagSmall RealD shift=1.5; -#if 1 +#if 0 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrySchur(src[0], maxIter, Nm, Nk, Nstop); // KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); @@ -356,10 +356,10 @@ int main (int argc, char ** argv) Nblock=LanParams.Nblock; bool if_verify=false; if(LanParams.verify) if_verify=true; - KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); - KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,true); -// BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); -// KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); +// KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); +// KrySchur(src, maxIter, Nm, Nk, Nstop,true,if_verify); + BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); + KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); // HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall); // KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true); #endif From 82a8f1342e61c2f8a8739f77b118b46693cea1f9 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Mon, 6 Apr 2026 13:35:43 -0400 Subject: [PATCH 2/2] Adding Lanczos parameter file in examples directory --- examples/Example_krylov_schur.cc | 3 ++- examples/LanParams.xml | 20 ++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) create mode 100644 examples/LanParams.xml diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 11f8a9fc4..df6e06ede 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -351,6 +351,7 @@ int main (int argc, char ** argv) KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrySchur(src[0], maxIter, Nm, Nk, Nstop); // KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); + std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #else int Nblock=4; Nblock=LanParams.Nblock; @@ -362,8 +363,8 @@ int main (int argc, char ** argv) KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); // HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall); // KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true); + std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #endif - std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl; src[0]=KrySchur.evecs[0]; std::cout << GridLogMessage << "KrySchur.evecs= "<< KrySchur.evecs.size() < + + + 0 + -0.025 + 1.8 + 48 + 80 + 100 + 100 + 0 + 1000 + 4 + 0 + 1e-10 + 1 + 100 + 51 + +