From 84707cc3a8dbca1444ef771bbf213b8bd6607058 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 2 Apr 2026 15:22:21 -0400 Subject: [PATCH] Changing definition of Nm and Nk in BlockKrylovSchur --- Grid/algorithms/iterative/BlockKrylovSchur.h | 60 ++++++++++--------- .../iterative/HarmonicBlockKrylovSchur.h | 38 ++++++------ tests/lanczos/Test_blocked_krylov_schur.cc | 10 ++-- .../lanczos/Test_harmonic_vs_krylov_schur.cc | 12 ++-- 4 files changed, 62 insertions(+), 58 deletions(-) diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index e12604a70..fa86e241e 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -44,10 +44,10 @@ NAMESPACE_BEGIN(Grid); * A V_k = V_k H_k + F_k B_k^dag * * where - * V_k = Nm*Nblock orthonormal basis vectors (stored flat in basis[]) - * H_k = (Nm*Nblock) x (Nm*Nblock) upper block-Hessenberg Rayleigh quotient + * V_k = Nm orthonormal basis vectors (stored flat in basis[]) + * H_k = Nm x Nm upper block-Hessenberg Rayleigh quotient * F_k = Nblock residual vectors (the next block beyond V_k) - * B_k = (Nm*Nblock) x Nblock coupling matrix (non-zero only in last Nblock rows) + * B_k = Nm x Nblock coupling matrix (non-zero only in last Nblock rows) * * Each block Arnoldi step applies A to each of the Nblock vectors in the * current block, orthogonalises against all previous basis vectors, and @@ -56,15 +56,15 @@ NAMESPACE_BEGIN(Grid); * * The restart is a thick restart via the Schur decomposition of H_k: * H_k = Q^dag S Q - * The leading Nk*Nblock Schur vectors (chosen by RitzFilter) are retained, + * The leading Nk Schur vectors (chosen by RitzFilter) are retained, * the basis and Rayleigh quotient are truncated, and block Arnoldi continues - * from the Nk-th block. + * from the (Nk/Nblock)-th block. * * Parameters * ---------- * Nblock : block size p - * Nm : number of block steps (total Krylov dimension = Nm * Nblock) - * Nk : number of block steps to keep after each restart (Nk < Nm) + * Nm : total Krylov dimension (must be divisible by Nblock) + * Nk : total vectors to keep after each restart (must be divisible by Nblock, Nk < Nm) * Nstop : declare convergence when this many eigenpairs have converged * MaxIter : maximum number of outer (restart) iterations * Tolerance : relative convergence tolerance (||r|| < Tolerance * |lambda_max|) @@ -83,8 +83,8 @@ protected: // Parameters (set by operator()) //-------------------------------------------------------------------- int Nblock; // block size - int Nm; // block steps (total dim = Nm * Nblock) - int Nk; // blocks retained after restart + int Nm; // total Krylov dimension (multiple of Nblock) + int Nk; // total vectors retained after restart (multiple of Nblock) int Nstop; int MaxIter; RealD Tolerance; @@ -97,19 +97,19 @@ protected: RitzFilter ritzFilter; // Flat storage: basis[s*Nblock + t] is the t-th vector of block s - // After construction: basis has Nm*Nblock entries + // After construction: basis has Nm entries std::vector basis; - // Rayleigh quotient (Nm*Nblock) x (Nm*Nblock) + // Rayleigh quotient Nm x Nm CMat H; - // Residual block: Nblock vectors (the (Nm+1)-th block, unnormalised before + // Residual block: Nblock vectors (the (Nm/Nblock+1)-th block, unnormalised before // QR; normalised and orthogonalised as part of block Arnoldi) std::vector F; - // Coupling matrix B: (Nm*Nblock) x Nblock. + // Coupling matrix B: Nm x Nblock. // In exact arithmetic only the last Nblock rows are non-zero: - // B(Nm*Nblock - Nblock + t, s) = H_{Nm+1, Nm}(t, s) (the subdiagonal block) + // B(Nm - Nblock + t, s) = H_{Nm/Nblock+1, Nm/Nblock}(t, s) (the subdiagonal block) // We keep it as a full matrix for generality after restarts. CMat B; @@ -118,7 +118,7 @@ protected: // Output CVec evals; - CMat littleEvecs; // Nm*Nblock columns + CMat littleEvecs; // Nm columns std::vector ritzEstimates; public: @@ -144,8 +144,8 @@ public: * ---------- * v0 : block of Nblock starting vectors (size >= Nblock) * _maxIter : maximum outer (restart) iterations - * _Nm : number of block steps per cycle - * _Nk : number of block steps to keep after restart (Nk < Nm) + * _Nm : total Krylov dimension (must be divisible by _Nblock) + * _Nk : total vectors to keep after restart (must be divisible by _Nblock, _Nk < _Nm) * _Nstop : stop after _Nstop eigenvalues converged * _Nblock : block size */ @@ -160,10 +160,12 @@ public: Nblock = _Nblock; assert((int)v0.size() >= Nblock); + assert(Nm % Nblock == 0); + assert(Nk % Nblock == 0); assert(Nk < Nm); preRun(); // hook: derived classes add parameter assertions here - int N = Nm * Nblock; // total Krylov dimension + int N = Nm; // total Krylov dimension // Approximate largest eigenvalue for tolerance normalisation RealD approxLambdaMax = approxMaxEval(v0[0]); @@ -181,11 +183,11 @@ public: for (int iter = 0; iter < MaxIter; iter++) { std::cout << GridLogMessage << "BlockKrylovSchur: restart iteration " << iter << std::endl; - // ---- Block Arnoldi: extend from block start to block Nm ---- - blockArnoldiIteration(startBlock, Nm, start, doubleOrthog); + // ---- Block Arnoldi: extend from block start to block Nm/Nblock ---- + blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog); - // After first full cycle start from block Nk - start = Nk; + // After first full cycle start from block Nk/Nblock + start = Nk/Nblock; if (doVerify) { std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi"; @@ -196,8 +198,8 @@ public: ComplexSchurDecomposition schur(H, false, ritzFilter); std::cout << GridLogMessage << "BlockKrylovSchur: Schur decomposed." << std::endl; - // Reorder: bring wanted Nk*Nblock Schur values to top-left - schur.schurReorder(Nk * Nblock); + // Reorder: bring wanted Nk Schur values to top-left + schur.schurReorder(Nk); std::cout << GridLogMessage << "BlockKrylovSchur: Schur reordered." << std::endl; CMat Q = schur.getMatrixQ(); @@ -212,8 +214,8 @@ public: B = Q * B; H = schur.getMatrixS(); - // ---- Truncate to Nk*Nblock ---- - int Nkeep = Nk * Nblock; + // ---- Truncate to Nk ---- + int Nkeep = Nk; CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); H = CMat::Zero(N, N); @@ -449,7 +451,7 @@ private: void blockArnoldiIteration(std::vector& startBlock, int endBlock, int startIdx, bool doubleOrthog) { - int N = Nm * Nblock; + int N = Nm; if (startIdx == 0) { basis.clear(); @@ -523,7 +525,7 @@ private: { int kBase = k * Nblock; // first flat index of current block int prevN = kBase + Nblock; // number of basis vectors so far after this step - int N = Nm * Nblock; + int N = Nm; // W[t] = op(basis[kBase + t]) — dispatches to applyBlock() virtual std::vector W(Nblock, Field(Grid_)); @@ -547,7 +549,7 @@ private: // Store residual block F F = W; - if (k == Nm - 1) { + if (k == Nm/Nblock - 1) { // Last block: compute coupling matrix B for KS decomp. // // blockQR modifies F in-place (F → Q orthonormal) and returns R diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 78f660910..6b7daca67 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -49,9 +49,9 @@ NAMESPACE_BEGIN(Grid); * * A V = V H + F B^dag (1) * - * with V orthonormal (Nm*Nblock columns), H the (Nm*Nblock)² block + * with V orthonormal (Nm columns), H the Nm² block * upper-Hessenberg Rayleigh quotient, F the Nblock residual vectors and B - * the (Nm*Nblock)×Nblock coupling matrix, the harmonic Rayleigh quotient + * the Nm×Nblock coupling matrix, the harmonic Rayleigh quotient * relative to shift σ is * * Hhat = H + (H - σI)^{-H} B B^H (2) @@ -94,8 +94,8 @@ NAMESPACE_BEGIN(Grid); * ---------- * shift : target shift σ (default 0.0) * Nblock : block size p - * Nm : number of block steps (total dim = Nm * Nblock) - * Nk : blocks to retain after each restart (Nk < Nm) + * Nm : total Krylov dimension (must be divisible by Nblock) + * Nk : total vectors to retain after each restart (must be divisible by Nblock, Nk < Nm) * Nstop : stop when this many eigenpairs converge * MaxIter : maximum outer (restart) iterations * Tolerance: relative convergence tolerance @@ -119,8 +119,8 @@ class HarmonicBlockKrylovSchur { // Parameters //-------------------------------------------------------------------- int Nblock; - int Nm; - int Nk; + int Nm; // total Krylov dimension (multiple of Nblock) + int Nk; // total vectors retained after restart (multiple of Nblock) int Nstop; int MaxIter; RealD Tolerance; @@ -173,9 +173,11 @@ public: Nblock = _Nblock; assert((int)v0.size() >= Nblock); + assert(Nm % Nblock == 0); + assert(Nk % Nblock == 0); assert(Nk < Nm); - int N = Nm * Nblock; + int N = Nm; RealD approxLambdaMax = approxMaxEval(v0[0]); rtol = Tolerance * approxLambdaMax; @@ -194,9 +196,9 @@ public: std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: restart iteration " << iter << std::endl; - // ---- Block Arnoldi: extend from block 'start' to block Nm ---- - blockArnoldiIteration(startBlock, Nm, start, doubleOrthog); - start = Nk; + // ---- Block Arnoldi: extend from block 'start' to block Nm/Nblock ---- + blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog); + start = Nk/Nblock; if (doVerify) { std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi"; @@ -209,12 +211,12 @@ public: // ---- Schur decompose Hhat ---- ComplexSchurDecomposition schur(Hhat, false, ritzFilter); - schur.schurReorder(Nk * Nblock); + schur.schurReorder(Nk); std::cout << GridLogMessage - << "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk*Nblock):" << std::endl; + << "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk):" << std::endl; CMat S = schur.getMatrixS(); - for (int i = 0; i < Nk * Nblock; i++) + for (int i = 0; i < Nk; i++) std::cout << GridLogMessage << " [" << i << "] " << S(i, i) << std::endl; CMat Q = schur.getMatrixQ(); @@ -229,8 +231,8 @@ public: H = Q * H * Qt; B = Q * B; - // ---- Truncate to Nk*Nblock ---- - int Nkeep = Nk * Nblock; + // ---- Truncate to Nk ---- + int Nkeep = Nk; CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); H = CMat::Zero(N, N); @@ -440,7 +442,7 @@ private: void blockArnoldiIteration(std::vector& startBlock, int endBlock, int startIdx, bool doubleOrthog) { - int N = Nm * Nblock; + int N = Nm; if (startIdx == 0) { basis.clear(); @@ -492,7 +494,7 @@ private: { int kBase = k * Nblock; int prevN = kBase + Nblock; - int N = Nm * Nblock; + int N = Nm; std::vector W(Nblock, Field(Grid_)); for (int t = 0; t < Nblock; t++) @@ -514,7 +516,7 @@ private: F = W; - if (k == Nm - 1) { + if (k == Nm/Nblock - 1) { // Last block: record coupling in B as R^H (Hermitian conjugate of QR factor) // KS relation requires B[kBase+t, s] = conj(R[s,t]) CMat R = blockQR(F); diff --git a/tests/lanczos/Test_blocked_krylov_schur.cc b/tests/lanczos/Test_blocked_krylov_schur.cc index a277a1185..1e716d085 100644 --- a/tests/lanczos/Test_blocked_krylov_schur.cc +++ b/tests/lanczos/Test_blocked_krylov_schur.cc @@ -68,7 +68,7 @@ int main(int argc, char** argv) int maxIter, int Nstop) { std::cout << GridLogMessage << "===== " << label << " =====" << std::endl; - BlockedKrylovSchur bks(op, grid, 1e-6, EvalReSmall); + BlockKrylovSchur bks(op, grid, 1e-6, EvalReSmall); std::vector v0(Nblock, Field(grid)); for (int t = 0; t < Nblock; t++) random(RNG, v0[t]); @@ -79,17 +79,17 @@ int main(int argc, char** argv) std::cout << GridLogMessage << label << " done." << std::endl; }; - // Test 1: Nblock=1 — scalar case, regression + // Test 1: Nblock=1 — scalar case, regression (Nm,Nk now total vectors) runTest("Nblock=1 Nm=10 Nk=5 maxIter=3", 1, 10, 5, 3, 5); // Test 2: Nblock=2 — exercises the B^H fix for off-diagonal elements - runTest("Nblock=2 Nm=8 Nk=4 maxIter=3", 2, 8, 4, 3, 4); + runTest("Nblock=2 Nm=16 Nk=8 maxIter=3", 2, 16, 8, 3, 4); // Test 3: Nblock=3 — further stress-test the B^H fix - runTest("Nblock=3 Nm=9 Nk=3 maxIter=3", 3, 9, 3, 3, 3); + runTest("Nblock=3 Nm=27 Nk=9 maxIter=3", 3, 27, 9, 3, 3); // Test 4: Nblock=2, larger cycle — more restarts - runTest("Nblock=2 Nm=12 Nk=6 maxIter=5", 2, 12, 6, 5, 6); + runTest("Nblock=2 Nm=24 Nk=12 maxIter=5", 2, 24, 12, 5, 6); if (nFail == 0) std::cout << GridLogMessage << "All BlockedKrylovSchur tests completed." << std::endl; diff --git a/tests/lanczos/Test_harmonic_vs_krylov_schur.cc b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc index 98bd3b12c..da2a35058 100644 --- a/tests/lanczos/Test_harmonic_vs_krylov_schur.cc +++ b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc @@ -62,8 +62,8 @@ int main(int argc, char** argv) // Parameters (kept small so output is readable) //---------------------------------------------------------------------- const int Nblock = 2; - const int Nm = 6; - const int Nk = 3; + const int Nm = 12; // total vectors (= 6 blocks * Nblock=2) + const int Nk = 6; // total kept (= 3 blocks * Nblock=2) const int Nstop = 2; const int maxIter = 4; const RealD tol = 1e-6; @@ -82,12 +82,12 @@ int main(int argc, char** argv) std::cout << GridLogMessage << "\n========================================" << std::endl; std::cout << GridLogMessage - << " BlockedKrylovSchur (Nblock=" << Nblock + << " BlockKrylovSchur (Nblock=" << Nblock << " Nm=" << Nm << " Nk=" << Nk << ")" << std::endl; std::cout << GridLogMessage << "========================================\n" << std::endl; - BlockedKrylovSchur bks(op, grid, tol, EvalReSmall); + BlockKrylovSchur bks(op, grid, tol, EvalReSmall); bks(v0, maxIter, Nm, Nk, Nstop, Nblock, /*doubleOrthog=*/true, /*doVerify=*/true); @@ -103,12 +103,12 @@ int main(int argc, char** argv) std::cout << GridLogMessage << "\n========================================" << std::endl; std::cout << GridLogMessage - << " HarmonicBlockedKrylovSchur (Nblock=" << Nblock + << " HarmonicBlockKrylovSchur (Nblock=" << Nblock << " Nm=" << Nm << " Nk=" << Nk << " shift=0)" << std::endl; std::cout << GridLogMessage << "========================================\n" << std::endl; - HarmonicBlockedKrylovSchur hbks(op, grid, tol, 0.0, EvalNormSmall); + HarmonicBlockKrylovSchur hbks(op, grid, tol, 0.0, EvalNormSmall); hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock, /*doubleOrthog=*/true, /*doVerify=*/true);