From af1fb41d55921d64aa6fa7c5bac7086cff537b9b Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Wed, 8 Apr 2026 22:50:37 -0400 Subject: [PATCH 1/5] Precision reduction in Block KS verify() for better readability Adding Plaquette stat testing --- Grid/algorithms/iterative/BlockKrylovSchur.h | 6 +- .../iterative/HarmonicBlockKrylovSchur.h | 6 +- tests/core/Test_plaquette_stats.cc | 211 ++++++++++++++++++ 3 files changed, 217 insertions(+), 6 deletions(-) create mode 100644 tests/core/Test_plaquette_stats.cc diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index c3835c88d..c48c8b360 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -329,7 +329,7 @@ public: 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::setprecision(4) << std::setw(14) << H(i, j); std::cout << std::endl; } @@ -337,7 +337,7 @@ public: 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::setprecision(4) << std::setw(14) << B(i, t); std::cout << std::endl; } @@ -353,7 +353,7 @@ public: std::cout << GridLogMessage << "M = (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << M(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << M(i, j); std::cout << std::endl; } diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 6b7daca67..117ef6a07 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -326,7 +326,7 @@ public: std::cout << GridLogMessage << "H (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << H(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << H(i, j); std::cout << std::endl; } @@ -334,7 +334,7 @@ public: std::cout << GridLogMessage << "B (" << nBasis << " x " << nF << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int t = 0; t < nF; t++) - std::cout << " " << std::setw(14) << B(i, t); + std::cout << " " << std::setprecision(4) << std::setw(14) << B(i, t); std::cout << std::endl; } @@ -350,7 +350,7 @@ public: std::cout << GridLogMessage << "M = (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << M(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << M(i, j); std::cout << std::endl; } diff --git a/tests/core/Test_plaquette_stats.cc b/tests/core/Test_plaquette_stats.cc new file mode 100644 index 000000000..10b33535c --- /dev/null +++ b/tests/core/Test_plaquette_stats.cc @@ -0,0 +1,211 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: tests/core/Test_plaquette_stats.cc + +Copyright (C) 2015 + +Author: Chulwoo Jung + +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 */ + +/** + * Test_plaquette_stats + * + * Measure every plaquette Re Tr[U_mu(x) U_nu(x+mu) U_mu†(x+nu) U_nu†(x)] / Nc + * across all sites and all (mu,nu) planes and report max, min, and average. + * + * Usage: + * ./Test_plaquette_stats [Grid options] [--file ] [--hot] + * + * --file Read gauge field from NERSC-format file + * --hot Use random (hot) SU(3) start (default: cold/unit start) + * + * Grid size defaults to 8^3 x 16; override with --grid (e.g. --grid 4.4.4.8). + */ + +#include + +using namespace std; +using namespace Grid; + +// Return the plane label string for output +static std::string planeName(int mu, int nu) +{ + const char dirs[] = "xyzt"; + std::string s; + s += dirs[mu]; + s += dirs[nu]; + return s; +} + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + // ---- Lattice geometry ---- + Coordinate latt_size = GridDefaultLatt(); + if (latt_size.size() == 0) { + latt_size = {8, 8, 8, 16}; + } + Coordinate simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + GridCartesian grid(latt_size, simd_layout, mpi_layout); + + std::cout << GridLogMessage << "Lattice: "; + for (int d = 0; d < Nd; d++) + std::cout << latt_size[d] << (d < Nd-1 ? "x" : "\n"); + + // ---- Gauge field ---- + LatticeGaugeField Umu(&grid); + + // Check for --file argument + std::string config_file = ""; + for (int i = 1; i < argc - 1; i++) { + if (std::string(argv[i]) == "--file") { + config_file = argv[i+1]; + break; + } + } + + bool doHot = false; + for (int i = 1; i < argc; i++) { + if (std::string(argv[i]) == "--hot") { doHot = true; break; } + } + + if (!config_file.empty()) { + std::cout << GridLogMessage << "Reading gauge field from " << config_file << std::endl; + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, config_file); + } else { + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&grid); + pRNG.SeedFixedIntegers(seeds); + if (doHot) { + std::cout << GridLogMessage << "Generating hot (random SU(3)) start" << std::endl; + SU::HotConfiguration(pRNG, Umu); + } else { + std::cout << GridLogMessage << "Using cold (unit) gauge start" << std::endl; + SU::ColdConfiguration(pRNG, Umu); + } + } + + // ---- Extract link matrices ---- + std::vector U(Nd, &grid); + for (int mu = 0; mu < Nd; mu++) + U[mu] = PeekIndex(Umu, mu); + + // ---- Per-plane plaquette statistics ---- + // + // For each (mu, nu) plane (mu > nu) compute + // P_munu(x) = Re Tr[plaquette] / Nc + // then report max, min, mean over all sites. + // + // WilsonLoops::traceDirPlaquette gives Tr[U_mu U_nu U_mu† U_nu†] (complex). + + double vol = grid.gSites(); + + // Accumulate site-average plaquette (sum over planes / Nplanes / Nc) + LatticeComplex plaq_all(&grid); + plaq_all = Zero(); + int Nplanes = 0; + + std::cout << GridLogMessage + << "======== Per-plane plaquette statistics ========" << std::endl; + std::cout << GridLogMessage + << std::setw(6) << "plane" + << std::setw(20) << "max" + << std::setw(20) << "min" + << std::setw(20) << "average" + << std::endl; + + for (int mu = 1; mu < Nd; mu++) { + for (int nu = 0; nu < mu; nu++) { + + // Per-site trace of plaquette in (mu,nu) plane + LatticeComplex sitePlaq(&grid); + ColourWilsonLoops::traceDirPlaquette(sitePlaq, U, mu, nu); + + plaq_all = plaq_all + sitePlaq; + Nplanes++; + + // --- global average via sum() --- + TComplex Tsum = sum(sitePlaq); + ComplexD csum = TensorRemove(Tsum); + RealD avg = csum.real() / vol / Nc; + + // --- global max and min via unvectorize + GlobalMax/GlobalMin --- + std::vector sv; + unvectorizeToLexOrdArray(sv, sitePlaq); + + RealD local_max = -1e38, local_min = 1e38; + for (auto& tc : sv) { + RealD r = TensorRemove(tc).real() / Nc; + if (r > local_max) local_max = r; + if (r < local_min) local_min = r; + } + + // Reduce across MPI ranks + grid.GlobalMax(local_max); + grid.GlobalMin(local_min); + + std::cout << GridLogMessage + << std::setw(6) << planeName(mu, nu) + << std::setw(20) << std::setprecision(10) << local_max + << std::setw(20) << std::setprecision(10) << local_min + << std::setw(20) << std::setprecision(10) << avg + << std::endl; + } + } + + // ---- Overall (averaged over all planes) statistics ---- + // plaq_all = sum of Tr[...] over all 6 (mu,nu) planes + // Normalise to Re Tr / Nc per plane + plaq_all = plaq_all * (1.0 / Nc / Nplanes); + + TComplex Tsum_all = sum(plaq_all); + RealD avg_all = TensorRemove(Tsum_all).real() / vol; + + std::vector sv_all; + unvectorizeToLexOrdArray(sv_all, plaq_all); + + RealD max_all = -1e38, min_all = 1e38; + for (auto& tc : sv_all) { + RealD r = TensorRemove(tc).real(); + if (r > max_all) max_all = r; + if (r < min_all) min_all = r; + } + grid.GlobalMax(max_all); + grid.GlobalMin(min_all); + + // Cross-check with built-in avgPlaquette + RealD avg_builtin = ColourWilsonLoops::avgPlaquette(Umu); + + std::cout << GridLogMessage + << "======== Overall plaquette statistics (all planes) ========" << std::endl; + std::cout << GridLogMessage << " max = " << max_all << std::endl; + std::cout << GridLogMessage << " min = " << min_all << std::endl; + std::cout << GridLogMessage << " average = " << avg_all << std::endl; + std::cout << GridLogMessage << " avgPlaquette (builtin check) = " << avg_builtin << std::endl; + + Grid_finalize(); + return 0; +} From 99be36c891b52376b0f86cf9f04bb990734b5cf2 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Wed, 8 Apr 2026 23:37:26 -0400 Subject: [PATCH 2/5] Fixing tests/core/Test_plaquette_stats --- tests/core/Test_plaquette_stats.cc | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/tests/core/Test_plaquette_stats.cc b/tests/core/Test_plaquette_stats.cc index 10b33535c..2ba4c57be 100644 --- a/tests/core/Test_plaquette_stats.cc +++ b/tests/core/Test_plaquette_stats.cc @@ -63,7 +63,7 @@ int main(int argc, char** argv) // ---- Lattice geometry ---- Coordinate latt_size = GridDefaultLatt(); if (latt_size.size() == 0) { - latt_size = {8, 8, 8, 16}; + latt_size = Coordinate(std::vector{8, 8, 8, 16}); } Coordinate simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); Coordinate mpi_layout = GridDefaultMpi(); @@ -163,9 +163,11 @@ int main(int argc, char** argv) if (r < local_min) local_min = r; } - // Reduce across MPI ranks + // Reduce across MPI ranks (no GlobalMin; negate to use GlobalMax) grid.GlobalMax(local_max); - grid.GlobalMin(local_min); + local_min = -local_min; + grid.GlobalMax(local_min); + local_min = -local_min; std::cout << GridLogMessage << std::setw(6) << planeName(mu, nu) @@ -194,7 +196,9 @@ int main(int argc, char** argv) if (r < min_all) min_all = r; } grid.GlobalMax(max_all); - grid.GlobalMin(min_all); + min_all = -min_all; + grid.GlobalMax(min_all); + min_all = -min_all; // Cross-check with built-in avgPlaquette RealD avg_builtin = ColourWilsonLoops::avgPlaquette(Umu); From aecc50869c7162d7c36df8c68a493652fa18e389 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 9 Apr 2026 21:13:41 -0400 Subject: [PATCH 3/5] Another try at block Harmonic KS. Still not working --- .../iterative/HarmonicBlockKrylovSchur.h | 72 ++++++++++++++----- examples/Example_krylov_schur.cc | 7 +- 2 files changed, 57 insertions(+), 22 deletions(-) diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 117ef6a07..bbfe9ebd6 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -234,6 +234,41 @@ public: // ---- Truncate to Nk ---- int Nkeep = Nk; + // ---- Option B: corrected restart starting block ------------------------- + // + // In standard Krylov-Schur, Q diagonalises H itself, so Q H Q^H = S is + // upper triangular and the off-diagonal coupling block H_dk = S[Nkeep:,:Nkeep] = 0. + // Truncation is therefore exact. + // + // Here Q diagonalises Hhat (not H), so H_new = Q H Q^H is generally dense. + // The off-diagonal block + // + // H_dk = H_new[Nkeep:, :Nkeep] (size (N-Nkeep) x Nkeep) + // + // is non-zero, and the true KS relation after truncation is: + // + // A V_k = V_k H_k + V_disc H_dk + F B_k^H + // + // If we restart from F only, the V_disc H_dk term is lost. + // + // Fix: include the V_disc coupling in the new starting block: + // + // G[t] = F[t] + sum_{s=Nkeep}^{N-1} basis[s] * H(s, t), t = 0..Nblock-1 + // + // G[:,t] is the t-th column of (V_disc H_dk + F B_k^H) restricted to the + // first Nblock columns of H_dk. Since F ⊥ V_k and V_disc ⊥ V_k, G is + // automatically orthogonal to the retained subspace. + // + // This must be computed BEFORE basis[Nkeep:] and H[Nkeep:, :] are discarded. + + std::vector G(Nblock, Field(Grid_)); + for (int t = 0; t < Nblock; t++) { + G[t] = F[t]; + for (int s = Nkeep; s < N; s++) + G[t] += basis[s] * H(s, t); + } + blockQR(G); // orthonormalise within G (G is already ⊥ to basis[0:Nkeep]) + CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); H = CMat::Zero(N, N); H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)) = Htmp; @@ -249,10 +284,10 @@ public: std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl; - // Restart from the residual block F (unchanged from last Arnoldi step). - // Note: for a Hermitian operator the correct H rows H[i,j] for i >= Nkeep+Nblock, - // j < Nkeep are filled via Hermitian symmetry inside blockArnoldiStep. - startBlock = F; + // Use corrected starting block G (not bare F). + // G encodes the coupling from the discarded basis vectors V_disc through + // H_dk[:,0:Nblock], restoring the exact KS relation to first order. + startBlock = G; if (doVerify) { std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation"; @@ -454,25 +489,26 @@ private: blockOrthonormalise(V0); for (auto& v : V0) basis.push_back(v); } else { - // Append residual block (startBlock = F_old) to basis. - // The truncated KS relation after restart is: + // Append the new starting block to the retained basis. // - // A V_k = V_k S_k + F_old B_old^dag (*) + // Standard KS (startBlock = F): + // The exact truncated relation is A V_k = V_k H_k + F B_k^dag, + // so the coupling rows are H[Nkeep+t, j] = conj(B_k[j,t]). // - // where V_k = basis[0:Nkeep], S_k is stored in H[0:Nkeep,0:Nkeep], - // B_old = B[0:Nkeep,:], F_old = startBlock. - // - // Once F_old is appended as basis[Nkeep:Nkeep+Nblock], (*) becomes - // a statement about the extended H matrix: - // - // H[Nkeep+t, j] = (B_old^dag)[t,j] = conj(B_old[j,t]) - // for t=0..Nblock-1, j=0..Nkeep-1 - // - // These "restart coupling rows" must be set before Arnoldi continues. + // Harmonic KS Option B (startBlock = G, where G = F + V_disc H_dk[:,0:Nblock]): + // The exact coupling rows are H[Nkeep+t, j] = , + // which differs from conj(B_k[j,t]) because G ≠ F. + // For a Hermitian operator these preset rows are overwritten by the + // Hermitian symmetry fill inside blockArnoldiStep (via explicit inner + // products), so the approximate preset below does no harm. + // For a non-Hermitian operator the preset is approximate; a more + // expensive fix would compute explicitly here. int Nkeep = startIdx * Nblock; for (auto& v : startBlock) basis.push_back(v); - // Fill restart coupling rows into H + // Fill restart coupling rows into H (exact for standard KS; approximate + // for harmonic KS with Option-B starting block, but overwritten by + // Hermitian symmetry fill for Hermitian operators). for (int t = 0; t < Nblock; t++) for (int j = 0; j < Nkeep; j++) H(Nkeep + t, j) = std::conj(B(j, t)); diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index df6e06ede..cad93a4d5 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -349,8 +349,8 @@ int main (int argc, char ** argv) RealD shift=1.5; #if 0 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); - KrySchur(src[0], maxIter, Nm, Nk, Nstop); -// KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); +// 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; @@ -360,9 +360,8 @@ int main (int argc, char ** argv) // 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); + KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #endif From a696953485c8d879ce2d55739154e3665850626b Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Tue, 14 Apr 2026 14:40:55 -0400 Subject: [PATCH 4/5] Finally Block Harmonic(?) KS working --- .../iterative/HarmonicBlockKrylovSchur.h | 206 ++++-------------- examples/Example_krylov_schur.cc | 22 +- .../lanczos/Test_harmonic_vs_krylov_schur.cc | 16 +- 3 files changed, 59 insertions(+), 185 deletions(-) diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index bbfe9ebd6..9a179c1ff 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -34,65 +34,44 @@ See the full license in the file "LICENSE" in the top level distribution directo NAMESPACE_BEGIN(Grid); /** - * Block harmonic restarted Krylov-Schur eigensolver. + * Block shift-targeted Krylov-Schur eigensolver. * - * Harmonic Ritz values - * -------------------- - * Standard Ritz values of A in a Krylov space K_m minimise the residual - * in a Galerkin sense; they are good approximations to eigenvalues at the - * *exterior* of the spectrum. For eigenvalues *near* a target shift σ - * (e.g. the smallest eigenvalues when σ=0) harmonic Ritz values are - * better-suited: they are obtained by a Petrov-Galerkin condition that - * requires the residual to be orthogonal to (A-σI)K_m instead of K_m. - * - * Given the block Arnoldi factorisation + * Algorithm + * --------- + * Uses a block Arnoldi factorisation: * * A V = V H + F B^dag (1) * - * with V orthonormal (Nm columns), H the Nm² block - * upper-Hessenberg Rayleigh quotient, F the Nblock residual vectors and B - * the Nm×Nblock coupling matrix, the harmonic Rayleigh quotient - * relative to shift σ is - * - * Hhat = H + (H - σI)^{-H} B B^H (2) - * - * Derivation: the harmonic Ritz condition (A-σI)Vy ⊥ (A-σI)V leads to - * - * [ (H-σI)^H (H-σI) + B B^H ] y = μ (H-σI)^H y - * - * Left-multiplying by (H-σI)^{-H} and setting θ = μ + σ gives the - * standard eigenvalue problem Hhat y = θ y with Hhat as in (2). - * - * The harmonic Ritz values θ_j are eigenvalues of Hhat; among these, - * the ones closest to σ (smallest |θ_j - σ|) are the best approximations - * to the eigenvalues of A near σ. + * with V orthonormal (Nm columns), H the Nm×Nm block upper-Hessenberg + * Rayleigh quotient, F the Nblock residual vectors and B the Nm×Nblock + * coupling matrix. * * Thick restart * ------------- - * The Schur decomposition Hhat = Q^dag S Q is computed and the - * leading Nk*Nblock Schur values (sorted by the RitzFilter) are kept. - * The same unitary rotation Q is applied to both the Krylov basis and - * to the *original* Rayleigh quotient H (not Hhat) for the restart: + * To target eigenvalues near shift σ, the Schur decomposition is computed + * for the shifted Rayleigh quotient: * - * V_new = V Q^dag [first Nk*Nblock columns] - * H_new = Q H Q^dag [truncated Nk*Nblock × Nk*Nblock] - * B_new = Q B [truncated Nk*Nblock × Nblock] + * (H - σI) = Q^dag S Q (2) * - * Block Arnoldi then resumes from block Nk, restoring H to full size - * as new columns are appended. + * Sorting the Schur values of (H - σI) by smallest |S(i,i)| = |λ - σ| and + * retaining the leading Nk is equivalent to selecting the Ritz values of H + * closest to σ. Since Q diagonalises H - σI (and hence H itself), the + * rotated Rayleigh quotient is exactly upper triangular: + * + * H_new = Q H Q^dag = S + σI (upper triangular) (3) + * + * Truncation to Nk is therefore exact: the off-diagonal coupling block + * H_new[Nk:, :Nk] = 0 by triangularity. The Krylov-Schur relation after + * restart is exact and block Arnoldi resumes cleanly from F. * * Convergence * ----------- - * For a harmonic Ritz pair (θ, y) the true Ritz residual bound is - * - * || (A - θI) V y || ≤ || B^H y || - * - * (same as for standard Ritz, because B captures the full coupling). - * Convergence is declared when || B^H y || < Tolerance * approxLambdaMax. + * Convergence is declared when || B^H y_k || < Tolerance * approxLambdaMax + * for each Ritz pair (λ_k, y_k) of the truncated H. * * Parameters * ---------- - * shift : target shift σ (default 0.0) + * shift : target shift σ (default 0.0); Schur values sorted by |λ - σ| * Nblock : block size p * Nm : total Krylov dimension (must be divisible by Nblock) * Nk : total vectors to retain after each restart (must be divisible by Nblock, Nk < Nm) @@ -205,70 +184,38 @@ public: verify(lbl); } - // ---- Form harmonic Rayleigh quotient ---- - // Hhat = H + (H - σI)^{-H} * B * B^H - CMat Hhat = harmonicRayleigh(H, B, N); - - // ---- Schur decompose Hhat ---- - ComplexSchurDecomposition schur(Hhat, false, ritzFilter); + // ---- Schur decompose (H - σI) to select Schur vectors closest to σ ---- + // Sorting the Schur values of (H - σI) by |S(i,i)| = |λ - σ| gives the + // Ritz values of H nearest the target shift without any matrix inversion. + // Because Q diagonalises (H - σI), it also diagonalises H: + // H_new = Q H Q^dag = S + σI (upper triangular) + // Truncation is therefore exact. + CMat Hshift = H - shift * CMat::Identity(N, N); + ComplexSchurDecomposition schur(Hshift, false, ritzFilter); schur.schurReorder(Nk); std::cout << GridLogMessage - << "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk):" << std::endl; + << "HarmonicBlockKrylovSchur: Ritz values nearest shift (first Nk):" << std::endl; CMat S = schur.getMatrixS(); for (int i = 0; i < Nk; i++) - std::cout << GridLogMessage << " [" << i << "] " << S(i, i) << std::endl; + std::cout << GridLogMessage << " [" << i << "] " << S(i, i) + shift << std::endl; CMat Q = schur.getMatrixQ(); CMat Qt = Q.adjoint(); - // ---- Rotate Krylov basis using Q from Hhat ---- + // ---- Rotate Krylov basis ---- std::vector basis2; constructUR(basis2, basis, Qt, N); basis = basis2; - // ---- Update H and B (rotate H, not Hhat) ---- - H = Q * H * Qt; + // ---- Update H and B ---- + // H_new = S + σI is upper triangular; off-diagonal block H_new[Nk:,:Nk] = 0 + H = S + shift * CMat::Identity(N, N); B = Q * B; - // ---- Truncate to Nk ---- + // ---- Truncate to Nk (exact: H upper triangular) ---- int Nkeep = Nk; - // ---- Option B: corrected restart starting block ------------------------- - // - // In standard Krylov-Schur, Q diagonalises H itself, so Q H Q^H = S is - // upper triangular and the off-diagonal coupling block H_dk = S[Nkeep:,:Nkeep] = 0. - // Truncation is therefore exact. - // - // Here Q diagonalises Hhat (not H), so H_new = Q H Q^H is generally dense. - // The off-diagonal block - // - // H_dk = H_new[Nkeep:, :Nkeep] (size (N-Nkeep) x Nkeep) - // - // is non-zero, and the true KS relation after truncation is: - // - // A V_k = V_k H_k + V_disc H_dk + F B_k^H - // - // If we restart from F only, the V_disc H_dk term is lost. - // - // Fix: include the V_disc coupling in the new starting block: - // - // G[t] = F[t] + sum_{s=Nkeep}^{N-1} basis[s] * H(s, t), t = 0..Nblock-1 - // - // G[:,t] is the t-th column of (V_disc H_dk + F B_k^H) restricted to the - // first Nblock columns of H_dk. Since F ⊥ V_k and V_disc ⊥ V_k, G is - // automatically orthogonal to the retained subspace. - // - // This must be computed BEFORE basis[Nkeep:] and H[Nkeep:, :] are discarded. - - std::vector G(Nblock, Field(Grid_)); - for (int t = 0; t < Nblock; t++) { - G[t] = F[t]; - for (int s = Nkeep; s < N; s++) - G[t] += basis[s] * H(s, t); - } - blockQR(G); // orthonormalise within G (G is already ⊥ to basis[0:Nkeep]) - CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); H = CMat::Zero(N, N); H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)) = Htmp; @@ -284,10 +231,8 @@ public: std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl; - // Use corrected starting block G (not bare F). - // G encodes the coupling from the discarded basis vectors V_disc through - // H_dk[:,0:Nblock], restoring the exact KS relation to first order. - startBlock = G; + // Restart from F (exact: no discarded-basis correction needed) + startBlock = F; if (doVerify) { std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation"; @@ -327,8 +272,7 @@ public: * A V = V H + F B^dag (KS) * * by explicit operator applications. H here is the standard Rayleigh - * quotient (not Hhat), so the KS relation is the same as for - * BlockKrylovSchur. + * quotient, so the KS relation is the same as for BlockKrylovSchur. * * Prints: * - H (current Rayleigh quotient, nBasis × nBasis) @@ -445,32 +389,6 @@ public: private: - //-------------------------------------------------------------------- - // Harmonic Rayleigh quotient - //-------------------------------------------------------------------- - /** - * Forms the harmonic Rayleigh quotient relative to shift σ: - * - * Hhat = H + (H - σI)^{-H} * B * B^H - * - * where H is the N×N block-Hessenberg, B is the N×Nblock coupling matrix. - * - * The N×N solve (H - σI)^H X = B B^H is done via Eigen's LU - * factorisation. If H - σI is (nearly) singular the result is - * ill-conditioned; in that case σ should be perturbed slightly. - */ - CMat harmonicRayleigh(const CMat& H_, const CMat& B_, int N) - { - CMat K = H_ - shift * CMat::Identity(N, N); - CMat KH = K.adjoint(); // (H - σI)^H - - // Solve KH * X = B B^H → X = KH^{-1} B B^H = (H-σI)^{-H} B B^H - CMat BBH = B_ * B_.adjoint(); // N × N - CMat X = KH.lu().solve(BBH); // N × N - - return H_ + X; - } - //-------------------------------------------------------------------- // Block Arnoldi iteration //-------------------------------------------------------------------- @@ -491,24 +409,15 @@ private: } else { // Append the new starting block to the retained basis. // - // Standard KS (startBlock = F): - // The exact truncated relation is A V_k = V_k H_k + F B_k^dag, - // so the coupling rows are H[Nkeep+t, j] = conj(B_k[j,t]). - // - // Harmonic KS Option B (startBlock = G, where G = F + V_disc H_dk[:,0:Nblock]): - // The exact coupling rows are H[Nkeep+t, j] = , - // which differs from conj(B_k[j,t]) because G ≠ F. - // For a Hermitian operator these preset rows are overwritten by the - // Hermitian symmetry fill inside blockArnoldiStep (via explicit inner - // products), so the approximate preset below does no harm. - // For a non-Hermitian operator the preset is approximate; a more - // expensive fix would compute explicitly here. + // The exact truncated KS relation is A V_k = V_k H_k + F B_k^dag, + // so the coupling rows are H[Nkeep+t, j] = conj(B_k[j,t]). + // Since H_new = S + σI is upper triangular, the off-diagonal block + // H_new[Nkeep:, :Nkeep] = 0 and the restart from F is exact. int Nkeep = startIdx * Nblock; for (auto& v : startBlock) basis.push_back(v); - // Fill restart coupling rows into H (exact for standard KS; approximate - // for harmonic KS with Option-B starting block, but overwritten by - // Hermitian symmetry fill for Hermitian operators). + // Fill restart coupling rows into H (exact: H_new is upper triangular, + // so B encodes the only non-zero coupling to the new block). for (int t = 0; t < Nblock; t++) for (int j = 0; j < Nkeep; j++) H(Nkeep + t, j) = std::conj(B(j, t)); @@ -560,10 +469,6 @@ private: for (int s = 0; s < Nblock; s++) B(kBase + t, s) = std::conj(R(s, t)); // B_block = R^H beta_k = R.norm(); - // Hermitian symmetry fill for last block (same as non-last path below) - for (int t = 0; t < Nblock; t++) - for (int j = 0; j < kBase; j++) - H(kBase + t, j) = std::conj(H(j, kBase + t)); return; } @@ -578,23 +483,6 @@ private: for (int t = 0; t < Nblock; t++) basis.push_back(F[t]); - // Hermitian symmetry fill: H[kBase+t, j] = conj(H[j, kBase+t]) for j < kBase. - // - // In a fresh block Arnoldi the Krylov structure forces H[kBase+t, j] = 0 for - // j < kBase-Nblock (sub-subdiagonal), so this is a no-op. - // - // After a non-Schur restart (e.g. harmonic restart where H_new = Q H Q^dag is - // a full matrix), A v_k_j for j < Nkeep has components in ALL new extended - // vectors, making these elements non-zero. The Arnoldi step fills column - // kBase+t (H[j, kBase+t] for j < prevN) via inner products, but never fills - // the corresponding row. For a Hermitian operator the two are related by - // H[kBase+t, j] = - // = conj() = conj(H[j, kBase+t]) - // Filling these ensures H = H^dag and fixes the M != H discrepancy that - // corrupts subsequent Arnoldi steps after a harmonic restart. - for (int t = 0; t < Nblock; t++) - for (int j = 0; j < kBase; j++) - H(kBase + t, j) = std::conj(H(j, kBase + t)); } //-------------------------------------------------------------------- diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index cad93a4d5..991b204ce 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -56,6 +56,7 @@ struct LanczosParameters: Serializable { Integer, maxIter, Integer, Nblock, Integer, verify, + RealD, shift , RealD, resid, RealD, ChebyLow, RealD, ChebyHigh, @@ -117,12 +118,9 @@ public: InvertNonHermitianLinearOperator(Matrix &Mat,RealD stp=1e-8): _Mat(Mat),_stp(stp){}; // Support for coarsening to a multigrid void OpDiag (const Field &in, Field &out) { -// _Mat.Mdiag(in,out); -// out = out + shift*in; assert(0); } void OpDir (const Field &in, Field &out,int dir,int disp) { -// _Mat.Mdir(in,out,dir,disp); assert(0); } void OpDirAll (const Field &in, std::vector &out){ @@ -131,11 +129,6 @@ public: }; void Op (const Field &in, Field &out){ Field tmp(in.Grid()); -// _Mat.M(in,out); -// RealD mass=-shift; -// WilsonCloverFermionD Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); -// NonHermitianLinearOperator HermOp(_Mat); -// BiCGSTAB CG(_stp,10000); _Mat.Mdag(in,tmp); MdagMLinearOperator HermOp(_Mat); ConjugateGradient CG(_stp,10000); @@ -341,12 +334,7 @@ int main (int argc, char ** argv) // Run KrylovSchur and Arnoldi on a Hermitian matrix std::cout << GridLogMessage << "Running Krylov Schur" << std::endl; - // KrylovSchur KrySchur (Dsq, FGrid, 1e-8, EvalNormLarge); -// KrylovSchur KrySchur (Dsq, FGrid, 1e-8,EvalImNormSmall); -// KrySchur(src, maxIter, Nm, Nk, Nstop); -// KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); -// Hacked, really EvalImagSmall - RealD shift=1.5; + RealD shift=LanParams.shift; #if 0 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrySchur(src[0], maxIter, Nm, Nk, Nstop); @@ -357,10 +345,8 @@ 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,true,if_verify); - BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); -// HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall); +// BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); + HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall); KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #endif diff --git a/tests/lanczos/Test_harmonic_vs_krylov_schur.cc b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc index da2a35058..0eebb342f 100644 --- a/tests/lanczos/Test_harmonic_vs_krylov_schur.cc +++ b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc @@ -61,11 +61,11 @@ int main(int argc, char** argv) //---------------------------------------------------------------------- // Parameters (kept small so output is readable) //---------------------------------------------------------------------- - const int Nblock = 2; - 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 int Nblock = 4; + const int Nm = 20; // total vectors (= 5 blocks * Nblock=4) + const int Nk = 8; // total kept (= 2 blocks * Nblock=4) + const int Nstop = 4; + const int maxIter = 8; const RealD tol = 1e-6; // Two identical starting blocks @@ -73,7 +73,7 @@ int main(int argc, char** argv) std::vector v0b(Nblock, Field(grid)); for (int t = 0; t < Nblock; t++) { random(RNG, v0[t]); - v0b[t] = v0[t]; + v0b[t] = v0[t]; // identical start for fair comparison } //---------------------------------------------------------------------- @@ -87,7 +87,7 @@ int main(int argc, char** argv) std::cout << GridLogMessage << "========================================\n" << std::endl; - BlockKrylovSchur bks(op, grid, tol, EvalReSmall); + BlockKrylovSchur bks(op, grid, tol, EvalImNormSmall); bks(v0, maxIter, Nm, Nk, Nstop, Nblock, /*doubleOrthog=*/true, /*doVerify=*/true); @@ -108,7 +108,7 @@ int main(int argc, char** argv) std::cout << GridLogMessage << "========================================\n" << std::endl; - HarmonicBlockKrylovSchur hbks(op, grid, tol, 0.0, EvalNormSmall); + HarmonicBlockKrylovSchur hbks(op, grid, tol, 0.0, EvalImNormSmall); hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock, /*doubleOrthog=*/true, /*doVerify=*/true); From 3b28e3c98ffaa1dd0933fed3bb7a61af9d248e44 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Fri, 17 Apr 2026 11:08:34 -0400 Subject: [PATCH 5/5] Added flags for explicit eval checks at the end of operator() --- Grid/algorithms/iterative/BlockKrylovSchur.h | 16 ++++++++++++++++ .../iterative/HarmonicBlockKrylovSchur.h | 16 ++++++++++++++++ Grid/algorithms/iterative/KrylovSchur.h | 19 +++++++++++++++++-- examples/Example_krylov_schur.cc | 2 +- 4 files changed, 50 insertions(+), 3 deletions(-) diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index c48c8b360..6c599f612 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -123,6 +123,7 @@ protected: public: std::vector evecs; + bool doEvalCheck = false; //-------------------------------------------------------------------- // Constructor @@ -252,6 +253,21 @@ public: std::cout << GridLogMessage << "BlockKrylovSchur: done after " << iter << " restarts, " << Nconv << " converged." << std::endl; std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl; + + if (doEvalCheck) { + Field w(Grid_); + for (int k = 0; k < (int)evecs.size(); k++) { + Linop.Op(evecs[k], w); + ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w)); + w -= eval_est * evecs[k]; + RealD res = std::sqrt(norm2(w)); + std::cout << GridLogMessage << "BlockKrylovSchur: evec[" << k << "]" + << " eval_reported = " << evals[k] + << " eval_est = " << eval_est + << " || A v - eval_est * v || = " << res << std::endl; + } + } + return; } } diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 9a179c1ff..1b405f975 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -125,6 +125,7 @@ class HarmonicBlockKrylovSchur { public: std::vector evecs; + bool doEvalCheck = false; //-------------------------------------------------------------------- // Constructor @@ -253,6 +254,21 @@ public: << "HarmonicBlockKrylovSchur: done after " << iter << " restarts, " << Nconv << " converged." << std::endl; std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl; + + if (doEvalCheck) { + Field w(Grid_); + for (int k = 0; k < (int)evecs.size(); k++) { + Linop.Op(evecs[k], w); + ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w)); + w -= eval_est * evecs[k]; + RealD res = std::sqrt(norm2(w)); + std::cout << GridLogMessage << "HarmonicBlockKrylovSchur: evec[" << k << "]" + << " eval_reported = " << evals[k] + << " eval_est = " << eval_est + << " || A v - eval_est * v || = " << res << std::endl; + } + } + return; } } diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 0ebf29d24..bd97e8529 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -336,8 +336,9 @@ class KrylovSchur { RitzFilter ritzFilter; // how to sort evals - public: + public: RealD *shift; + bool doEvalCheck = false; 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), @@ -601,11 +602,25 @@ if (!shift){ int Nconv = converged(); std::cout << GridLogMessage << "Number of evecs converged: " << Nconv << std::endl; if (Nconv >= Nstop || i == MaxIter - 1) { - std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " + std::cout << GridLogMessage << "Converged with " << Nconv << " / " << Nstop << " eigenvectors on iteration " << i << "." << std::endl; // basisRotate(evecs, Qt, 0, Nk, 0, Nk, Nm); // Think this might have been the issue std::cout << GridLogMessage << "Eigenvalues: " << evals << std::endl; + if (doEvalCheck) { + Field w(Grid); + for (int k = 0; k < (int)evecs.size(); k++) { + Linop.Op(evecs[k], w); + ComplexD eval_est = toStdCmplx(innerProduct(evecs[k], w)); + w -= eval_est * evecs[k]; + RealD res = std::sqrt(norm2(w)); + std::cout << GridLogMessage << "KrylovSchur: evec[" << k << "]" + << " eval_reported = " << evals[k] + << " eval_est = " << eval_est + << " || A v - eval_est * v || = " << res << std::endl; + } + } + // writeEigensystem(path); return; diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 991b204ce..08a1e3b04 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -335,7 +335,7 @@ int main (int argc, char ** argv) // Run KrylovSchur and Arnoldi on a Hermitian matrix std::cout << GridLogMessage << "Running Krylov Schur" << std::endl; RealD shift=LanParams.shift; -#if 0 +#if 1 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrySchur(src[0], maxIter, Nm, Nk, Nstop); KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);