From 7752e7ead1a18e4f05fd8c5124bd60f76f8c9b1a Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 23 Apr 2026 16:28:48 -0400 Subject: [PATCH] Sorting eval/evec according to RitzFilter --- Grid/algorithms/iterative/BlockKrylovSchur.h | 20 +++++++++- .../iterative/HarmonicBlockKrylovSchur.h | 20 +++++++++- Grid/algorithms/iterative/KrylovSchur.h | 39 ++++++++++--------- examples/Example_krylov_schur.cc | 6 ++- examples/LanParams.xml | 2 + 5 files changed, 62 insertions(+), 25 deletions(-) diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index 6c599f612..968942126 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -30,6 +30,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #define GRID_BLOCKED_KRYLOV_SCHUR_H #include +#include NAMESPACE_BEGIN(Grid); @@ -697,8 +698,22 @@ private: { Eigen::ComplexEigenSolver es; es.compute(Hk); - evals = es.eigenvalues(); - littleEvecs = es.eigenvectors(); + + // Sort to match schurReorder ordering. + int n = es.eigenvalues().size(); + ComplexComparator cComp(ritzFilter); + std::vector idx(n); + std::iota(idx.begin(), idx.end(), 0); + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b))); + }); + + evals.resize(n); + littleEvecs.resize(n, n); + for (int k = 0; k < n; k++) { + evals(k) = es.eigenvalues()(idx[k]); + littleEvecs.col(k) = es.eigenvectors().col(idx[k]); + } evecs.clear(); for (int k = 0; k < Nkeep; k++) { @@ -735,6 +750,7 @@ private: std::cout << GridLogMessage << "BlockKrylovSchur: Ritz estimate[" << k << "] = " << res << " eval = " << evals[k] << std::endl; if (res < rtol) Nconv++; + else break; } return Nconv; } diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 1b405f975..9afeebff4 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -30,6 +30,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #define GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H #include +#include NAMESPACE_BEGIN(Grid); @@ -574,8 +575,22 @@ private: { Eigen::ComplexEigenSolver es; es.compute(Hk); - evals = es.eigenvalues(); - littleEvecs = es.eigenvectors(); + + // Sort to match schurReorder ordering. + int n = es.eigenvalues().size(); + ComplexComparator cComp(ritzFilter); + std::vector idx(n); + std::iota(idx.begin(), idx.end(), 0); + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b))); + }); + + evals.resize(n); + littleEvecs.resize(n, n); + for (int k = 0; k < n; k++) { + evals(k) = es.eigenvalues()(idx[k]); + littleEvecs.col(k) = es.eigenvectors().col(idx[k]); + } evecs.clear(); for (int k = 0; k < Nkeep; k++) { @@ -612,6 +627,7 @@ private: << "HarmonicBlockKrylovSchur: Ritz estimate[" << k << "] = " << res << " eval = " << evals[k] << std::endl; if (res < rtol) Nconv++; + else break; } return Nconv; } diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index bd97e8529..5504fd256 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -105,8 +105,8 @@ struct ComplexComparator bool operator()(std::complex z1, std::complex z2) { RealD tmp1=std::abs(std::imag(z1)); RealD tmp2=std::abs(std::imag(z2)); - if ( std::abs(std::real(z1)) >4.) tmp1 += 100.; - if ( std::abs(std::real(z2)) >4.) tmp2 += 100.; + if ( std::abs(std::real(z1)) >2.) tmp1 += 100.; + if ( std::abs(std::real(z2)) >2.) tmp2 += 100.; switch (RF) { case EvalNormSmall: return std::abs(z1) < std::abs(z2); @@ -738,35 +738,36 @@ if (!shift){ { std::cout << GridLogMessage << "Computing eigenvalues." << std::endl; - // 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; Eigen::ComplexEigenSolver es; - // es.compute(Rayleigh); es.compute(S); - evals = es.eigenvalues(); - littleEvecs = es.eigenvectors(); - // std::cout << GridLogDebug << "Little evecs: " << littleEvecs << std::endl; - // std::cout << "Rayleigh diag: " << S.diagonal() << std::endl; - // std::cout << "Rayleigh evals: " << evals << std::endl; + // Sort eigenvalues/evecs to match the schurReorder ordering. + int n = es.eigenvalues().size(); + ComplexComparator cComp(ritzFilter); + std::vector idx(n); + std::iota(idx.begin(), idx.end(), 0); + std::sort(idx.begin(), idx.end(), [&](int a, int b){ + return cComp(toStdCmplx(es.eigenvalues()(a)), toStdCmplx(es.eigenvalues()(b))); + }); + + evals.resize(n); + littleEvecs.resize(n, n); + for (int k = 0; k < n; k++) { + evals(k) = es.eigenvalues()(idx[k]); + littleEvecs.col(k) = es.eigenvectors().col(idx[k]); + } // Convert evecs to lattice fields for (int k = 0; k < n; k++) { Eigen::VectorXcd vec = littleEvecs.col(k); Field tmp (basis[0].Grid()); tmp = Zero(); - for (int j = 0; j < basis.size(); j++) { + for (int j = 0; j < (int)basis.size(); j++) { tmp = tmp + vec[j] * basis[j]; } evecs.push_back(tmp); - // evecs[k] = tmp; } } @@ -831,12 +832,12 @@ if (!shift){ 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++; + } else { + break; } - } // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) || // checkRitzEstimate(); diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 991b204ce..f2698f245 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -338,6 +338,7 @@ int main (int argc, char ** argv) #if 0 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrySchur(src[0], maxIter, Nm, Nk, Nstop); + KrySchur.doEvalCheck=true; KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #else @@ -345,8 +346,9 @@ 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(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); std::cout << GridLogMessage << "BlockKrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #endif diff --git a/examples/LanParams.xml b/examples/LanParams.xml index 41239a09c..ade8f411d 100644 --- a/examples/LanParams.xml +++ b/examples/LanParams.xml @@ -10,8 +10,10 @@ 100 0 1000 + 1 4 0 + 1.5 1e-10 1 100