From dcda74f924123eb52aafb138c8775b3564b3307d Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Thu, 18 Dec 2025 18:23:50 +0000 Subject: [PATCH] Timing info for schurReorder,etc --- Grid/algorithms/iterative/KrylovSchur.h | 49 +++++++++++++++++-------- examples/Example_krylov_schur.cc | 3 +- 2 files changed, 35 insertions(+), 17 deletions(-) diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index 86700083..fc3b8968 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -40,6 +40,7 @@ enum RitzFilter { EvalNormLarge, // Keep evals with largest norm EvalReSmall, // Keep evals with smallest real part EvalReLarge, // Keep evals with largest real part + EvalReNormLarge, // Keep evals with largest real part EvalImSmall, // Keep evals with smallest imaginary part EvalImLarge, // Keep evals with largest imaginary part EvalImNormSmall, // Keep evals with smallest |imaginary| part @@ -52,6 +53,7 @@ inline RitzFilter selectRitzFilter(std::string s) { if (s == "EvalNormLarge") { return EvalNormLarge; } else if (s == "EvalReSmall") { return EvalReSmall; } else if (s == "EvalReLarge") { return EvalReLarge; } else + if (s == "EvalReNormLarge") { return EvalReNormLarge; } else if (s == "EvalImSmall") { return EvalImSmall; } else if (s == "EvalImLarge") { return EvalImLarge; } else if (s == "EvalImNormSmall") { return EvalImNormSmall; } else @@ -70,6 +72,8 @@ inline std::string rfToString(RitzFilter RF) { return "EvalReSmall"; case EvalReLarge: return "EvalReLarge"; + case EvalReNormLarge: + return "EvalReNormLarge"; case EvalImSmall: return "EvalImSmall"; case EvalImLarge: @@ -89,6 +93,10 @@ struct ComplexComparator RitzFilter RF; ComplexComparator (RitzFilter _rf) : RF(_rf) {} 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.; switch (RF) { case EvalNormSmall: return std::abs(z1) < std::abs(z2); @@ -98,12 +106,15 @@ struct ComplexComparator return std::real(z1) < std::real(z2); // DELETE THE ABS HERE!!! case EvalReLarge: return std::real(z1) > std::real(z2); + case EvalReNormLarge: + return std::abs(std::real(z1)) < std::abs(std::real(z2)); case EvalImSmall: return std::imag(z1) < std::imag(z2); case EvalImLarge: return std::imag(z1) > std::imag(z2); case EvalImNormSmall: - return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); +// return std::abs(std::imag(z1)) < std::abs(std::imag(z2)); + return tmp1 < tmp2 ; case EvalImNormLarge: return std::abs(std::imag(z1)) > std::abs(std::imag(z2)); default: @@ -224,6 +235,8 @@ class ComplexSchurDecomposition { */ // void schurReorder(int Nk, std::function compare) { void schurReorder(int Nk) { + std::cout << GridLogMessage << "schurReorder start " << std::endl; + int nswap=0; for (int i = 0; i < Nk; i++) { for (int k = 0; k <= Nm - 2; k++) { int idx = Nm - 2 - k; @@ -232,9 +245,14 @@ class ComplexSchurDecomposition { // if (std::real(S(idx, idx)) > std::real(S(idx+1, idx+1))) { // sort by smallest real eigenvalue if ( cCompare(S(idx+1, idx+1), S(idx, idx)) ) { // sort by largest modulus of eigenvalue swapEvals(idx); + nswap++; } } } + std::cout << GridLogMessage << "schurReorder end "<< nswap << std::endl; + for (int i = 0; i < Nk; i++) { + std::cout << GridLogMessage << "schurR "< b = Q*b - // basisRotate(basis, Q, 0, Nm, 0, Nm, Nm); - // basisRotate(evecs, Q, 0, Nm, 0, Nm, Nm); std::vector basis2; - // basis2.reserve(Nm); - // for (int i = start; i < Nm; i++) { - // basis2.emplace_back(Grid); - // } constructUR(basis2, basis, Qt, Nm); basis = basis2; if(0){ @@ -513,10 +528,11 @@ if(0){ } } +} std::cout << GridLogMessage << "*** TRUNCATING FOR RESTART *** " << std::endl; - +if (!shift){ std::cout << GridLogDebug << "Rayleigh before truncation: " << std::endl << Rayleigh << std::endl; Eigen::MatrixXcd RayTmp = Rayleigh(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); @@ -536,6 +552,7 @@ if(0){ // Compute eigensystem of Rayleigh. Note the eigenvectors correspond to the sorted eigenvalues. computeEigensystem(Rayleigh); std::cout << GridLogMessage << "Eigenvalues (first Nk sorted): " << std::endl << evals << std::endl; +} if(shift){ Rayleigh = Btilde; @@ -627,7 +644,7 @@ if(0){ } if (doubleOrthog) { - std::cout << GridLogMessage << "Double orthogonalizing." << std::endl; + std::cout << GridLogDebug << "Double orthogonalizing." << std::endl; for (int j = 0; j < basis.size(); j++) { coeff = innerProduct(basis[j], w); // see if there is any residual component in basis[j] direction Rayleigh(j, i) += coeff; // if coeff is non-zero, adjust Rayleigh @@ -775,7 +792,7 @@ if(0){ } // Check that Ritz estimate is explicitly || D (Uy) - lambda (Uy) || - // checkRitzEstimate(); +// checkRitzEstimate(); return Nconv; } diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 38b5fd2d..04211cd0 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -540,8 +540,9 @@ int main (int argc, char ** argv) // Hacked, really EvalImagSmall KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // BlockKrylovSchur KrySchur (HermOp2, UGrid, Nu, resid,EvalNormSmall); - RealD shift=1.; + RealD shift=2.5; KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); +// KrySchur(src[0], maxIter, Nm, Nk, Nstop); std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl; src[0]=KrySchur.evecs[0];