diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 435d60bef..2f7c4b1e9 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -89,12 +89,13 @@ NAMESPACE_CHECK(multigrid); #include #include #include +#include #include //#include #include #include #include //#include -//#include +#include #endif diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index f22de77d2..1407b3714 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -92,6 +92,7 @@ NAMESPACE_BEGIN(Grid); template class HarmonicBlockKrylovSchur { +protected: typedef Eigen::MatrixXcd CMat; typedef Eigen::VectorXcd CVec; @@ -143,12 +144,14 @@ public: beta_k(0.0), rtol(0.0) {} + virtual ~HarmonicBlockKrylovSchur() = default; + //-------------------------------------------------------------------- // Main entry point //-------------------------------------------------------------------- - void operator()(const std::vector& v0, int _maxIter, int _Nm, int _Nk, - int _Nstop, int _Nblock = 1, bool doubleOrthog = true, - bool doVerify = false) + virtual void operator()(const std::vector& v0, int _maxIter, int _Nm, int _Nk, + int _Nstop, int _Nblock = 1, bool doubleOrthog = true, + bool doVerify = false) { MaxIter = _maxIter; Nm = _Nm; @@ -161,6 +164,8 @@ public: if (useParityFlip) divisor *= 2; if (useGamma5) divisor *= 2; assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor); + std::cout << GridLogMessage << "divisor= " << divisor << std::endl; + } assert(Nm % Nblock == 0); assert(Nk % Nblock == 0); @@ -179,31 +184,8 @@ public: H = CMat::Zero(N, N); B = CMat::Zero(N, Nblock); - int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1); int start = 0; - std::vector startBlock; - startBlock.reserve(Nblock); - for (int i = 0; i < Nblock / divisor; i++) { - std::vector group; - group.push_back(v0[i]); - if (useParityFlip) { - int n = (int)group.size(); - for (int j = 0; j < n; j++) { - Field fp(Grid_); - parityFlippedField(group[j], fp); - group.push_back(std::move(fp)); - } - } - if (useGamma5) { - int n = (int)group.size(); - for (int j = 0; j < n; j++) { - Field g5v(Grid_); - gamma5Func(group[j], g5v); - group.push_back(std::move(g5v)); - } - } - for (auto& f : group) startBlock.push_back(std::move(f)); - } + std::vector startBlock = expandStartBlock(v0); for (int iter = 0; iter < MaxIter; iter++) { std::cout << GridLogMessage @@ -448,7 +430,42 @@ public: << "======== end verify ========" << std::endl; } -private: +protected: + + //-------------------------------------------------------------------- + // Starting-block expansion (parity flip / gamma5 partners) + //-------------------------------------------------------------------- + // Expands Nblock/divisor seed vectors into the full Nblock starting + // block, pairing each seed with its parity-flipped and/or gamma5 + // partners. Requires Nblock, Grid_ and the flags to be set. + std::vector expandStartBlock(const std::vector& v0) + { + int divisor = (useParityFlip ? 2 : 1) * (useGamma5 ? 2 : 1); + std::vector startBlock; + startBlock.reserve(Nblock); + for (int i = 0; i < Nblock / divisor; i++) { + std::vector group; + group.push_back(v0[i]); + if (useParityFlip) { + int n = (int)group.size(); + for (int j = 0; j < n; j++) { + Field fp(Grid_); + parityFlippedField(group[j], fp); + group.push_back(std::move(fp)); + } + } + if (useGamma5) { + int n = (int)group.size(); + for (int j = 0; j < n; j++) { + Field g5v(Grid_); + gamma5Func(group[j], g5v); + group.push_back(std::move(g5v)); + } + } + for (auto& f : group) startBlock.push_back(std::move(f)); + } + return startBlock; + } //-------------------------------------------------------------------- // Block Arnoldi iteration diff --git a/Grid/algorithms/iterative/KrylovSchur.h b/Grid/algorithms/iterative/KrylovSchur.h index a5ed84c41..6a8b821ba 100644 --- a/Grid/algorithms/iterative/KrylovSchur.h +++ b/Grid/algorithms/iterative/KrylovSchur.h @@ -451,18 +451,6 @@ if(0){ Eigen::MatrixXcd Q_s = schurS.getMatrixQ(); Eigen::MatrixXcd Qt_s = Q_s.adjoint(); // TODO should Q be real? -#if 0 - std::cout << GridLogMessage << "Q_s" << Q_s < b = Q*b std::vector basis2; - constructUR(basis2, basis, Qt, Nm); + constructUR(basis2, basis, Qt, Nm,Nm); basis = basis2; if(0){ Field w(Grid); @@ -861,7 +851,7 @@ if (!shift){ // std::cout << GridLogDebug << "Rayleigh in KSDecomposition: " << std::endl << Rayleigh << std::endl; std::vector rotated = basis; - constructUR(rotated, basis, Rayleigh, k); // manually rotate + constructUR(rotated, basis, Rayleigh, k,k); // manually rotate // Eigen::MatrixXcd Rt = Rayleigh.adjoint(); // basisRotate(rotated, Rt, 0, k, 0, k, k); // UR @@ -957,7 +947,7 @@ if (!shift){ * Note that I believe this is equivalent to basisRotate(U, R.adjoint(), 0, N, 0, N, N), but I'm * not 100% sure (this will be slower and unoptimized though). */ - void constructUR(std::vector& UR, std::vector &U, Eigen::MatrixXcd& R, int N) { + void constructUR(std::vector& UR, std::vector &U, Eigen::MatrixXcd& R, int N, int N2) { Field tmp (Grid); UR.clear(); @@ -967,6 +957,7 @@ if (!shift){ for (int i = 0; i < N; i++) { tmp = Zero(); + if (i < N2) for (int j = 0; j < N; j++) { std::cout << GridLogDebug << "Adding R("< + +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 */ +#ifndef GRID_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H +#define GRID_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H + +#include + +NAMESPACE_BEGIN(Grid); + +/** + * True harmonic-Ritz block Krylov-Schur eigensolver. + * + * Difference from the base class + * ------------------------------ + * HarmonicBlockKrylovSchur is a *shift-targeted* Krylov-Schur: it Schur- + * decomposes (H - sigma I) and sorts the standard Ritz values by |lambda - + * sigma|. For interior eigenvalues of Hermitian indefinite operators the + * standard Ritz values can be spurious ("ghosts"), so the sort can select + * junk. This variant performs genuine harmonic Ritz extraction instead. + * + * Harmonic extraction + * ------------------- + * From the exact block Arnoldi relation + * + * A V = V H + F B^dag (1) + * + * harmonic Ritz pairs (theta, g) w.r.t. shift sigma satisfy the Petrov- + * Galerkin condition (A - sigma)Vg - (theta - sigma)Vg perp (A - sigma)V, + * which reduces (with Hs = H - sigma I) to the small eigenproblem + * + * Hhat g = (theta - sigma) g, Hhat = Hs + Hs^{-H} B B^dag (2) + * + * Exact thick restart + * ------------------- + * From (2), H g = theta g - Hs^{-H} B (B^dag g), so the residual of every + * harmonic pair lies in the Nblock-dimensional column space of the single + * block + * + * W = F - V Hs^{-H} B : + * A (V g) - theta (V g) = W (B^dag g) (3) + * + * Hence {V g_1 .. g_Nk} + cols(W) admits an EXACT truncated Krylov-Schur + * relation A Y = Y H_k + What Bc^dag with dense H_k, and block Arnoldi + * resumes cleanly from What. Krylov-Schur (Stewart, SIMAX 23(3), 2001) + * does not require H to be triangular after restart -- only that the + * decomposition holds exactly, which (3) guarantees. + * + * This is the block analogue of the GMRES-DR restart (Morgan, SISC 24(1), + * 2002; block version: Morgan, Appl. Numer. Math. 54, 2005). The exact + * composition "block harmonic Krylov-Schur" is a derivation from these + * ingredients rather than a published algorithm; the relation (1) is + * checked explicitly at every restart when doVerify=true, and identity (3) + * was validated to machine precision in a standalone Eigen prototype. + * + * Caveats + * ------- + * - If sigma is (numerically) an eigenvalue of H, Hs^{-H} is singular; + * an assert fires -- perturb the shift slightly. + * - Reported eigenvalues come from the eigensystem of the retained H_k + * (the true projected Rayleigh quotient Y^dag A Y), as in the base + * class, so convergence tests and reporting are inherited unchanged. + * + * Usage: identical to HarmonicBlockKrylovSchur. + */ +template +class TrueHarmonicBlockKrylovSchur : public HarmonicBlockKrylovSchur { + + typedef HarmonicBlockKrylovSchur Base; + typedef typename Base::CMat CMat; + typedef typename Base::CVec CVec; + + // Members of the dependent base class need explicit scope. + using Base::Nblock; + using Base::Nm; + using Base::Nk; + using Base::Nstop; + using Base::MaxIter; + using Base::Tolerance; + using Base::shift; + using Base::Linop; + using Base::Grid_; + using Base::ritzFilter; + using Base::basis; + using Base::H; + using Base::F; + using Base::B; + using Base::beta_k; + using Base::rtol; + using Base::evals; + using Base::littleEvecs; + using Base::ritzEstimates; + + using Base::expandStartBlock; + using Base::blockArnoldiIteration; + using Base::blockQR; + using Base::computeEigensystem; + using Base::converged; + using Base::approxMaxEval; + +public: + using Base::evecs; + using Base::doEvalCheck; + using Base::useParityFlip; + using Base::useGamma5; + using Base::gamma5Func; + using Base::verify; + + TrueHarmonicBlockKrylovSchur(LinearOperatorBase& _Linop, GridBase* _Grid, + RealD _Tolerance, ComplexD _shift = 0.0, + RitzFilter _rf = EvalNormSmall) + : Base(_Linop, _Grid, _Tolerance, _shift, _rf) + {} + + //-------------------------------------------------------------------- + // Main entry point (same signature and outer structure as the base; + // the Schur-sort restart is replaced by harmonicRestart()). + //-------------------------------------------------------------------- + virtual void operator()(const std::vector& v0, int _maxIter, int _Nm, int _Nk, + int _Nstop, int _Nblock = 1, bool doubleOrthog = true, + bool doVerify = false) override + { + MaxIter = _maxIter; + Nm = _Nm; + Nk = _Nk; + Nstop = _Nstop; + Nblock = _Nblock; + + { + int divisor = 1; + if (useParityFlip) divisor *= 2; + if (useGamma5) divisor *= 2; + assert(Nblock % divisor == 0 && (int)v0.size() >= Nblock / divisor); + } + assert(Nm % Nblock == 0); + assert(Nk % Nblock == 0); + assert(Nk < Nm); + if (useGamma5) assert(gamma5Func && "useGamma5: gamma5Func must be set"); + + int N = Nm; + + RealD approxLambdaMax = approxMaxEval(v0[0]); + rtol = Tolerance * approxLambdaMax; + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: approx max eval = " << approxLambdaMax + << ", rtol = " << rtol + << ", shift = " << shift << std::endl; + + H = CMat::Zero(N, N); + B = CMat::Zero(N, Nblock); + + int start = 0; + std::vector startBlock = expandStartBlock(v0); + + for (int iter = 0; iter < MaxIter; iter++) { + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: restart iteration " << iter << std::endl; + + // ---- Block Arnoldi: extend from block 'start' to block Nm/Nblock ---- + blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog); + std::cout << GridLogMessage << "blockArnoldiIteration done " << std::endl; + start = Nk/Nblock; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi"; + verify(lbl); + } + + // ---- Harmonic extraction + exact thick restart ---- + harmonicRestart(); + std::cout << GridLogMessage << "harmonicRestart done " << std::endl; + + // Restart from the residual block W (exact by identity (3)) + startBlock = F; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after harmonic restart"; + verify(lbl); + } + + // ---- Eigensystem of retained H_k = Y^dag A Y for convergence ---- + CMat Hk = H(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)); + computeEigensystem(Hk, Nk); + + int Nconv = converged(Nk); + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: converged " << Nconv + << " / " << Nstop << std::endl; + + if (Nconv >= Nstop || iter == MaxIter - 1) { + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: 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 << "TrueHarmonicBlockKrylovSchur: evec[" << k << "]" + << " eval_reported = " << evals[k] + << " eval_est = " << eval_est + << " || A v - eval_est * v || = " << res << std::endl; + } + } + + return; + } + } + } + +protected: + + //-------------------------------------------------------------------- + // Harmonic thick restart + //-------------------------------------------------------------------- + // On entry: A V = V H + F B^dag exact, V = basis (Nm fields), F Nblock + // residual fields, H/B the Nm-sized small matrices. + // On exit: basis = Y (Nk orthonormal fields spanning the harmonic + // Ritz space), F = What (Nblock orthonormal fields, perp Y), + // H top-left Nk x Nk = Y^dag A Y (dense), B top Nk rows the + // new coupling; relation A Y = Y H_k + F B^dag again exact. + void harmonicRestart() + { + int N = Nm; + + // ---- Small harmonic eigenproblem: Hhat g = (theta - sigma) g ---- + CMat Hs = H - shift * CMat::Identity(N, N); + Eigen::FullPivLU lu(Hs.adjoint()); + assert(lu.isInvertible() && + "TrueHarmonicBlockKrylovSchur: H - shift*I singular; perturb the shift"); + CMat M = lu.solve(B); // N x Nblock : Hs^{-H} B + CMat Hhat = Hs + M * B.adjoint(); // N x N + + Eigen::ComplexEigenSolver es(Hhat); + + // Sort harmonic values (theta - sigma) by the Ritz filter + // (EvalNormSmall -> closest to the shift). + 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))); + }); + + CMat G(N, Nk); + CVec theta(Nk); + for (int k = 0; k < Nk; k++) { + G.col(k) = es.eigenvectors().col(idx[k]); + theta(k) = es.eigenvalues()(idx[k]) + toStdCmplx(shift); + } + + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: harmonic Ritz values nearest shift:" << std::endl; + for (int k = 0; k < Nk; k++) + std::cout << GridLogMessage << " [" << k << "] " << theta(k) << std::endl; + + // ---- Orthonormalise the little vectors: G = Gt Rg ---- + Eigen::HouseholderQR qr(G); + CMat Gt = qr.householderQ() * CMat::Identity(N, Nk); + CMat Rg = qr.matrixQR().topLeftCorner(Nk, Nk).template triangularView(); + for (int k = 0; k < Nk; k++) + assert(std::abs(Rg(k, k)) > 1e-14 && + "TrueHarmonicBlockKrylovSchur: harmonic vectors linearly dependent"); + CMat RgInv = Rg.inverse(); + + // A (V G) = (V G) diag(theta) + W (B^dag G), W = F - V M + // => A (V Gt) = (V Gt) T + W Ct + CMat T = Rg * theta.asDiagonal() * RgInv; // Nk x Nk + CMat Ct = (B.adjoint() * G) * RgInv; // Nblock x Nk + + // ---- Large-field work ---- + // Y = V Gt : new orthonormal basis (Nk fields) + std::vector Y; + constructThin(Y, basis, Gt, Nk); + + // W = F - V M : Nblock residual directions containing every harmonic + // residual (identity (3) in the class comment) + std::vector W = F; + for (int t = 0; t < Nblock; t++) + for (int j = 0; j < N; j++) + W[t] -= basis[j] * M(j, t); + + // Orthogonalise W against Y (two passes), recording P = Y^dag W + CMat P = CMat::Zero(Nk, Nblock); + for (int pass = 0; pass < 2; pass++) { + for (int j = 0; j < Nk; j++) { + for (int t = 0; t < Nblock; t++) { + ComplexD c = innerProduct(Y[j], W[t]); + P(j, t) += toStdCmplx(c); + W[t] -= c * Y[j]; + } + } + } + + // Block QR of the projected residual: W <- What, returns Rw + CMat Rw = blockQR(W); + + // ---- Reassemble the exact Krylov-Schur relation ---- + // A Y = Y (T + P Ct) + What (Rw Ct) + CMat Hk = T + P * Ct; // Nk x Nk, dense + CMat Bc = Rw * Ct; // Nblock x Nk + + H = CMat::Zero(N, N); + H(Eigen::seqN(0, Nk), Eigen::seqN(0, Nk)) = Hk; + B = CMat::Zero(N, Nblock); + for (int j = 0; j < Nk; j++) + for (int t = 0; t < Nblock; t++) + B(j, t) = std::conj(Bc(t, j)); + + basis = Y; + F = W; + beta_k = Bc.norm(); + + std::cout << GridLogMessage + << "TrueHarmonicBlockKrylovSchur: beta_k = " << beta_k << std::endl; + } + + //-------------------------------------------------------------------- + // Thin basis combination: out[i] = sum_j U[j] * R(j,i), i < ncol + //-------------------------------------------------------------------- + void constructThin(std::vector& out, const std::vector& U, + const CMat& R, int ncol) + { + out.clear(); + Field tmp(Grid_); + for (int i = 0; i < ncol; i++) { + tmp = Zero(); + for (int j = 0; j < (int)U.size(); j++) + tmp += U[j] * R(j, i); + out.push_back(tmp); + } + } + +}; + +NAMESPACE_END(Grid); + +#endif // GRID_TRUE_HARMONIC_BLOCKED_KRYLOV_SCHUR_H diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 397f61365..208343106 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -334,13 +334,13 @@ int main (int argc, char ** argv) // Run KrylovSchur and Arnoldi on a Hermitian matrix RealD shift=LanParams.shift; -#if 0 +#if 1 std::cout << GridLogMessage << "Running Krylov Schur" << std::endl; KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrySchur(src[0], maxIter, Nm, Nk, Nstop); KrySchur.doEvalCheck=true; -// KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); - KrySchur(src[0], maxIter, Nm, Nk, Nstop); + KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); +// KrySchur(src[0], maxIter, Nm, Nk, Nstop); std::cout << GridLogMessage << "KrylovSchur evec.size= " << KrySchur.evecs.size()<< std::endl; #else std::cout << GridLogMessage << "Running BlockKrylovSchur" << std::endl; @@ -349,10 +349,19 @@ int main (int argc, char ** argv) bool if_verify=false; if(LanParams.verify) if_verify=true; // BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); - HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalNormSmall); + bool useTrueHarmonic = GridCmdOptionExists(argv, argv+argc, std::string("--true-harmonic")); + HarmonicBlockKrylovSchur KrySchurShift (Dwilson, UGrid, resid,shift,EvalImNormSmall); + TrueHarmonicBlockKrylovSchur KrySchurTrue (Dwilson, UGrid, resid,shift,EvalImNormSmall); + HarmonicBlockKrylovSchur& KrySchur = useTrueHarmonic + ? static_cast&>(KrySchurTrue) + : KrySchurShift; + std::cout << GridLogMessage + << (useTrueHarmonic ? "Using TrueHarmonicBlockKrylovSchur (harmonic Ritz)" + : "Using HarmonicBlockKrylovSchur (shift-sorted Ritz)") + << std::endl; KrySchur.doEvalCheck=true; - KrySchur.useParityFlip=true; std::cout << GridLogMessage << "useParityFlip= " <