diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index e7e1f043f..6c30beaad 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -85,8 +85,12 @@ NAMESPACE_CHECK(multigrid); #include #include +#include +#include #include #include #include +#include +#include #endif diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h new file mode 100644 index 000000000..8881bce86 --- /dev/null +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -0,0 +1,701 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/BlockKrylovSchur.h + +Copyright (C) 2015 + +Author: Peter Boyle +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 */ +#ifndef GRID_BLOCKED_KRYLOV_SCHUR_H +#define GRID_BLOCKED_KRYLOV_SCHUR_H + +#include + +NAMESPACE_BEGIN(Grid); + +/** + * Block (block-Arnoldi) restarted Krylov-Schur eigensolver for + * general non-Hermitian operators. + * + * Algorithm + * --------- + * Uses a block Arnoldi factorisation of block size Nblock: + * + * 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 + * 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) + * + * Each block Arnoldi step applies A to each of the Nblock vectors in the + * current block, orthogonalises against all previous basis vectors, and + * reduces the residual block to upper-triangular form via Householder QR + * (implemented here as modified Gram-Schmidt within the block). + * + * 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 basis and Rayleigh quotient are truncated, and block Arnoldi continues + * from the Nk-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) + * Nstop : declare convergence when this many eigenpairs have converged + * MaxIter : maximum number of outer (restart) iterations + * Tolerance : relative convergence tolerance (||r|| < Tolerance * |lambda_max|) + */ +template +class BlockKrylovSchur { + + //-------------------------------------------------------------------- + // Types + //-------------------------------------------------------------------- + typedef Eigen::MatrixXcd CMat; + typedef Eigen::VectorXcd CVec; + + //-------------------------------------------------------------------- + // Parameters (set by operator()) + //-------------------------------------------------------------------- + int Nblock; // block size + int Nm; // block steps (total dim = Nm * Nblock) + int Nk; // blocks retained after restart + int Nstop; + int MaxIter; + RealD Tolerance; + + //-------------------------------------------------------------------- + // Internal state + //-------------------------------------------------------------------- + LinearOperatorBase& Linop; + GridBase* Grid_; + RitzFilter ritzFilter; + + // Flat storage: basis[s*Nblock + t] is the t-th vector of block s + // After construction: basis has Nm*Nblock entries + std::vector basis; + + // Rayleigh quotient (Nm*Nblock) x (Nm*Nblock) + CMat H; + + // Residual block: Nblock vectors (the (Nm+1)-th block, unnormalised before + // QR; normalised and orthogonalised as part of block Arnoldi) + std::vector F; + + // Coupling matrix B: (Nm*Nblock) 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) + // We keep it as a full matrix for generality after restarts. + CMat B; + + RealD beta_k; // Frobenius norm of the last subdiagonal block + RealD rtol; // absolute tolerance = Tolerance * approxLambdaMax + + // Output + CVec evals; + CMat littleEvecs; // Nm*Nblock columns + std::vector ritzEstimates; + +public: + std::vector evecs; + + //-------------------------------------------------------------------- + // Constructor + //-------------------------------------------------------------------- + BlockKrylovSchur(LinearOperatorBase& _Linop, GridBase* _Grid, + RealD _Tolerance, RitzFilter _rf = EvalReSmall) + : Linop(_Linop), Grid_(_Grid), Tolerance(_Tolerance), ritzFilter(_rf), + Nblock(-1), Nm(-1), Nk(-1), Nstop(-1), MaxIter(-1), + beta_k(0.0), rtol(0.0) + {} + + //-------------------------------------------------------------------- + // Main entry point + //-------------------------------------------------------------------- + /** + * Run the blocked Krylov-Schur algorithm. + * + * Parameters + * ---------- + * 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) + * _Nstop : stop after _Nstop eigenvalues converged + * _Nblock : block size + */ + 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; + Nk = _Nk; + Nstop = _Nstop; + Nblock = _Nblock; + + assert((int)v0.size() >= Nblock); + assert(Nk < Nm); + + int N = Nm * Nblock; // total Krylov dimension + + // Approximate largest eigenvalue for tolerance normalisation + RealD approxLambdaMax = approxMaxEval(v0[0]); + rtol = Tolerance * approxLambdaMax; + std::cout << GridLogMessage << "BlockKrylovSchur: approx max eval = " + << approxLambdaMax << ", rtol = " << rtol << std::endl; + + // Initialise + H = CMat::Zero(N, N); + B = CMat::Zero(N, Nblock); + + int start = 0; + std::vector startBlock(v0.begin(), v0.begin() + Nblock); + + 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); + + // After first full cycle start from block Nk + start = Nk; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi"; + verify(lbl); + } + + // ---- Schur decompose H ---- + 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); + std::cout << GridLogMessage << "BlockKrylovSchur: Schur reordered." << std::endl; + + CMat Q = schur.getMatrixQ(); + CMat Qt = Q.adjoint(); + + // Rotate Krylov basis: basis_new[i] = sum_j basis[j] * Qt(j,i) + std::vector basis2; + constructUR(basis2, basis, Qt, N); + basis = basis2; + + // Update b and H + B = Q * B; + H = schur.getMatrixS(); + + // ---- Truncate to Nk*Nblock ---- + int Nkeep = Nk * Nblock; + + 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; + + std::vector basisTmp(basis.begin(), basis.begin() + Nkeep); + basis = basisTmp; + + CMat Btmp = B(Eigen::seqN(0, Nkeep), Eigen::all); + B = CMat::Zero(N, Nblock); + B(Eigen::seqN(0, Nkeep), Eigen::all) = Btmp; + + // beta_k = Frobenius norm of the effective coupling + beta_k = Btmp.norm(); + std::cout << GridLogMessage << "BlockKrylovSchur: beta_k = " << beta_k << std::endl; + + // Restart: the new starting block is F (the residual block from Arnoldi) + startBlock = F; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation"; + verify(lbl); + } + + // ---- Compute eigensystem of truncated H for convergence check ---- + CMat Hk = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); + computeEigensystem(Hk, Nkeep); + + int Nconv = converged(Nkeep); + std::cout << GridLogMessage << "BlockKrylovSchur: converged " << Nconv + << " / " << Nstop << std::endl; + + if (Nconv >= Nstop || iter == MaxIter - 1) { + std::cout << GridLogMessage << "BlockKrylovSchur: done after " << iter + << " restarts, " << Nconv << " converged." << std::endl; + std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl; + return; + } + } + } + + // Accessors + std::vector getEvecs() { return evecs; } + CVec getEvals() { return evals; } + std::vector getRitzEstimates() { return ritzEstimates; } + + //-------------------------------------------------------------------- + // Verification: print H and B, check A V = V H + F B^dag explicitly + //-------------------------------------------------------------------- + /** + * Checks the block Arnoldi / Krylov-Schur decomposition + * + * A V = V H + F B^dag (KS) + * + * by explicit operator applications. For each basis vector j: + * + * w_j = A basis[j] + * + * The nBasis × nBasis matrix M of inner products is computed: + * + * M[i, j] = + * + * and compared column-by-column against H. Separately, the nBasis × Nblock + * residual coupling matrix R is computed: + * + * R[j, t] = * ||F[t]|| (scaled by F-block norms) + * + * but since F is already normalised, R[j,t] = . + * + * The KS relation for column j reads: + * w_j = sum_i basis[i] H[i,j] + sum_t F[t] B[j,t]* + * so the deviation in column j is + * dev_j = w_j - sum_i basis[i] M[i,j] (should be zero for exact arithmetic) + * augmented by the F B^dag term in the last block. + * + * Prints: + * - H (current Rayleigh quotient, nBasis × nBasis) + * - B (coupling matrix, nBasis × Nblock) + * - M (explicit inner product matrix ) + * - max |H[i,j] - M[i,j]| (should be O(machine epsilon)) + * - for each basis column j: || A v_j - V H[:,j] - F B[j,:]^* || + * + * Parameters + * ---------- + * label : string printed at the start (e.g. "after restart 2") + */ + void verify(const std::string& label = "") + { + int nBasis = (int)basis.size(); + int nF = (int)F.size(); + + if (nBasis == 0) { + std::cout << GridLogMessage + << "BlockKrylovSchur::verify [" << label + << "]: basis is empty." << std::endl; + return; + } + + std::cout << GridLogMessage + << "======== BlockKrylovSchur::verify [" << label << "] ========" << std::endl; + std::cout << GridLogMessage + << " nBasis = " << nBasis << " Nblock = " << Nblock + << " nF = " << nF << std::endl; + + // ---- Print H ---- + 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::endl; + } + + // ---- Print B ---- + 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::endl; + } + + // ---- Compute M[i,j] = ---- + CMat M = CMat::Zero(nBasis, nBasis); + Field w(Grid_); + for (int j = 0; j < nBasis; j++) { + Linop.Op(basis[j], w); + for (int i = 0; i < nBasis; i++) + M(i, j) = toStdCmplx(innerProduct(basis[i], w)); + } + + 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::endl; + } + + // ---- max |H - M| ---- + RealD maxHM = 0.0; + for (int i = 0; i < nBasis; i++) + for (int j = 0; j < nBasis; j++) + maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j))); + std::cout << GridLogMessage + << " max |H[i,j] - M[i,j]| = " << maxHM << std::endl; + + // ---- Check orthonormality of basis ---- + CMat G = CMat::Zero(nBasis, nBasis); + for (int i = 0; i < nBasis; i++) + for (int j = 0; j < nBasis; j++) + G(i, j) = toStdCmplx(innerProduct(basis[i], basis[j])); + CMat Gerr = G - CMat::Identity(nBasis, nBasis); + std::cout << GridLogMessage + << " max | - delta_ij| = " << Gerr.cwiseAbs().maxCoeff() << std::endl; + + // ---- Per-column residual: || A v_j - V H[:,j] - F B[j,:]^* || ---- + // For each basis vector j, compute A v_j then subtract V H[:,j] and F B[j,:]^* + RealD maxColDev = 0.0; + for (int j = 0; j < nBasis; j++) { + Linop.Op(basis[j], w); + + // subtract V H[:,j] + for (int i = 0; i < nBasis; i++) + w -= basis[i] * H(i, j); + + // subtract F B[j,:]^* (F[t] * conj(B[j,t])) + for (int t = 0; t < nF; t++) + w -= F[t] * std::conj(B(j, t)); + + RealD dev = std::sqrt(norm2(w)); + std::cout << GridLogMessage + << " || A v[" << j << "] - V H[:,j] - F B[j,:]* || = " << dev << std::endl; + maxColDev = std::max(maxColDev, dev); + } + std::cout << GridLogMessage + << " max column deviation = " << maxColDev << std::endl; + + // ---- Check F block orthogonality against basis ---- + if (nF > 0) { + RealD maxFV = 0.0; + for (int t = 0; t < nF; t++) + for (int i = 0; i < nBasis; i++) { + RealD ip = std::abs(toStdCmplx(innerProduct(basis[i], F[t]))); + maxFV = std::max(maxFV, ip); + } + std::cout << GridLogMessage + << " max || (should be ~0) = " << maxFV << std::endl; + } + + std::cout << GridLogMessage + << "======== end verify ========" << std::endl; + } + +private: + + //-------------------------------------------------------------------- + // Block Arnoldi iteration + //-------------------------------------------------------------------- + /** + * Extends the block Arnoldi factorisation from block index 'start' to + * block index 'Nm'. + * + * On entry (start > 0): basis[0..start*Nblock-1] already set, + * H[0..start*Nblock-1, 0..start*Nblock-1] already set, + * B[start*Nblock-1, :] set (coupling from prior residual block). + * startBlock = the normalised residual block F from the previous cycle. + * + * On entry (start == 0): initialises everything from startBlock. + */ + void blockArnoldiIteration(std::vector& startBlock, int endBlock, + int startIdx, bool doubleOrthog) + { + int N = Nm * Nblock; + + if (startIdx == 0) { + basis.clear(); + F.clear(); + H = CMat::Zero(N, N); + B = CMat::Zero(N, Nblock); + + // Orthonormalise starting block via modified Gram-Schmidt + std::vector V0 = startBlock; + 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: + // + // A V_k = V_k S_k + F_old B_old^dag (*) + // + // 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 entries are the "restart coupling rows" that connect the new + // block to all retained Schur vectors and must be set before Arnoldi + // continues, otherwise A V_k = V_k H[:,0:Nkeep] would be missing the + // F_old B_old^dag term for those columns. + + int Nkeep = startIdx * Nblock; + for (auto& v : startBlock) basis.push_back(v); + + // Fill restart coupling rows into H + for (int t = 0; t < Nblock; t++) + for (int j = 0; j < Nkeep; j++) + H(Nkeep + t, j) = std::conj(B(j, t)); + + // Zero out B for the retained columns now that the coupling is in H + for (int j = 0; j < Nkeep; j++) + for (int t = 0; t < Nblock; t++) + B(j, t) = 0.0; + } + + // Main block Arnoldi loop + for (int k = startIdx; k < endBlock; k++) { + blockArnoldiStep(k, doubleOrthog); + } + } + + //-------------------------------------------------------------------- + /** + * One block Arnoldi step: extends by one block (Nblock vectors). + * + * Computes block column k of H and the next basis block V_{k+1}. + * + * Layout of basis (flat): + * basis[j*Nblock + t] = t-th vector of j-th block, j = 0..k + * + * After this call: + * H[i, k*Nblock : (k+1)*Nblock] filled for i = 0..(k+1)*Nblock - 1 + * basis[k*Nblock .. (k+1)*Nblock - 1] normalised (already set on entry) + * F = residual block (to become V_{k+1} after this step if k < Nm-1) + * + * If k < Nm-1, also: + * H[(k+1)*Nblock : (k+2)*Nblock, k*Nblock : (k+1)*Nblock] = subdiag block (from QR of residual) + * basis extended by Nblock (the normalised residual vectors) + */ + void blockArnoldiStep(int k, bool doubleOrthog) + { + 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; + + // W[t] = A * basis[kBase + t] + std::vector W(Nblock, Field(Grid_)); + for (int t = 0; t < Nblock; t++) { + Linop.Op(basis[kBase + t], W[t]); + } + + // Orthogonalise W against all current basis vectors (full reorthogonalisation) + // H[i, kBase + t] = + for (int pass = 0; pass < (doubleOrthog ? 2 : 1); pass++) { + for (int i = 0; i < prevN; i++) { + for (int t = 0; t < Nblock; t++) { + ComplexD coeff = innerProduct(basis[i], W[t]); + if (pass == 0) + H(i, kBase + t) = toStdCmplx(coeff); + else + H(i, kBase + t) += toStdCmplx(coeff); + W[t] -= coeff * basis[i]; + } + } + } + + // Store residual block F + F = W; + + if (k == Nm - 1) { + // Last block: compute coupling matrix B for KS decomp. + // + // blockQR modifies F in-place (F → Q orthonormal) and returns R + // such that W_orig[t] = sum_s F[s] * R[s,t] (W_orig = F_after * R). + // + // The KS relation for column j = kBase+t requires the coefficient of F[s] + // to be (B†)[s,j] = conj(B[j,s]). Matching with R[s,t]: + // conj(B[kBase+t, s]) = R[s,t] → B[kBase+t, s] = conj(R[s,t]) + // + // Equivalently the last Nblock rows of B are R^H (Hermitian conjugate of R). + // Note: for Nblock=1, R is scalar real positive, so this reduces to B = R. ✓ + CMat R = blockQR(F); // F is modified in-place to become Q; returns R + for (int t = 0; t < Nblock; t++) + for (int s = 0; s < Nblock; s++) + B(kBase + t, s) = std::conj(R(s, t)); // B_block = R^H + + beta_k = R.norm(); + return; + } + + // Not last block: QR-decompose residual to get V_{k+1} + CMat R = blockQR(F); // F orthonormalised in-place, R is upper triangular + + // Subdiagonal block of H: H[(k+1)*Nblock : (k+2)*Nblock, kBase : kBase+Nblock] = R + int nextBase = (k + 1) * Nblock; + for (int i = 0; i < Nblock; i++) + for (int j = 0; j < Nblock; j++) + H(nextBase + i, kBase + j) = R(i, j); + + // Append normalised residual block to basis + for (int t = 0; t < Nblock; t++) + basis.push_back(F[t]); + } + + //-------------------------------------------------------------------- + // Block QR via modified Gram-Schmidt within the block + //-------------------------------------------------------------------- + /** + * Given a block of Nblock vectors W (not necessarily orthonormal), + * orthonormalises them in-place and returns the upper-triangular R + * such that W_in = W_out * R. + * + * Handles (near-)linear dependence by zeroing vectors below threshold. + */ + CMat blockQR(std::vector& W) + { + CMat R = CMat::Zero(Nblock, Nblock); + const RealD deflThresh = 1e-14; + + for (int j = 0; j < Nblock; j++) { + // Orthogonalise W[j] against W[0..j-1] + for (int i = 0; i < j; i++) { + ComplexD coeff = innerProduct(W[i], W[j]); + R(i, j) = toStdCmplx(coeff); + W[j] -= coeff * W[i]; + } + RealD nrm = std::sqrt(norm2(W[j])); + R(j, j) = nrm; + if (nrm > deflThresh) { + W[j] *= (1.0 / nrm); + } else { + // deflation: zero this vector + W[j] = Zero(); + std::cout << GridLogMessage + << "BlockKrylovSchur: deflation at block column " << j + << " (norm = " << nrm << ")" << std::endl; + } + } + return R; + } + + //-------------------------------------------------------------------- + // Orthonormalise a block against itself (no prior basis) + //-------------------------------------------------------------------- + void blockOrthonormalise(std::vector& V) + { + for (int j = 0; j < (int)V.size(); j++) { + for (int i = 0; i < j; i++) { + ComplexD c = innerProduct(V[i], V[j]); + V[j] -= c * V[i]; + } + RealD nrm = std::sqrt(norm2(V[j])); + assert(nrm > 1e-14); + V[j] *= (1.0 / nrm); + } + } + + //-------------------------------------------------------------------- + // Basis rotation: UR[i] = sum_j U[j] * R(j, i) + //-------------------------------------------------------------------- + void constructUR(std::vector& UR, std::vector& U, + CMat& R, int N) + { + UR.clear(); + Field tmp(Grid_); + for (int i = 0; i < N; i++) { + tmp = Zero(); + for (int j = 0; j < N; j++) + tmp += U[j] * R(j, i); + UR.push_back(tmp); + } + } + + //-------------------------------------------------------------------- + // Eigensystem of the truncated Rayleigh quotient + //-------------------------------------------------------------------- + void computeEigensystem(CMat& Hk, int Nkeep) + { + Eigen::ComplexEigenSolver es; + es.compute(Hk); + evals = es.eigenvalues(); + littleEvecs = es.eigenvectors(); + + evecs.clear(); + for (int k = 0; k < Nkeep; k++) { + CVec vec = littleEvecs.col(k); + Field tmp(Grid_); + tmp = Zero(); + for (int j = 0; j < (int)basis.size() && j < Nkeep; j++) + tmp += vec[j] * basis[j]; + evecs.push_back(tmp); + } + } + + //-------------------------------------------------------------------- + // Convergence check + //-------------------------------------------------------------------- + /** + * An eigenpair (lambda, y) is converged if the Ritz estimate + * r = || B^dag y || + * satisfies r < rtol. Here B is the (Nkeep x Nblock) coupling matrix + * and y is the little eigenvector (Nkeep-vector) of H. + */ + int converged(int Nkeep) + { + ritzEstimates.clear(); + int Nconv = 0; + + CMat Bk = B(Eigen::seqN(0, Nkeep), Eigen::all); // Nkeep x Nblock + + for (int k = 0; k < Nkeep; k++) { + CVec yk = littleEvecs.col(k); // Nkeep-vector + CVec Bty = Bk.adjoint() * yk; // Nblock-vector + RealD res = Bty.norm(); + ritzEstimates.push_back(res); + std::cout << GridLogMessage << "BlockKrylovSchur: Ritz estimate[" << k + << "] = " << res << " eval = " << evals[k] << std::endl; + if (res < rtol) Nconv++; + } + return Nconv; + } + + //-------------------------------------------------------------------- + // Approximate maximum eigenvalue via power iteration + //-------------------------------------------------------------------- + RealD approxMaxEval(const Field& v0, int MAX_ITER = 50) + { + assert(norm2(v0) > 1e-8); + RealD lam = 0.0, denom = std::sqrt(norm2(v0)); + Field vcur(Grid_), vtmp(Grid_); + vcur = v0; + for (int i = 0; i < MAX_ITER; i++) { + Linop.Op(vcur, vtmp); + vcur = vtmp; + RealD num = std::sqrt(norm2(vcur)); + lam = num / denom; + denom = num; + } + return lam; + } + +}; + +NAMESPACE_END(Grid); + +#endif // GRID_BLOCKED_KRYLOV_SCHUR_H diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h new file mode 100644 index 000000000..78f660910 --- /dev/null +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -0,0 +1,700 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/HarmonicBlockKrylovSchur.h + +Copyright (C) 2015 + +Author: Peter Boyle +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 */ +#ifndef GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H +#define GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H + +#include + +NAMESPACE_BEGIN(Grid); + +/** + * Block harmonic restarted 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 + * + * A V = V H + F B^dag (1) + * + * with V orthonormal (Nm*Nblock columns), H the (Nm*Nblock)² block + * upper-Hessenberg Rayleigh quotient, F the Nblock residual vectors and B + * the (Nm*Nblock)×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 σ. + * + * 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: + * + * 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] + * + * Block Arnoldi then resumes from block Nk, restoring H to full size + * as new columns are appended. + * + * 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. + * + * Parameters + * ---------- + * 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) + * Nstop : stop when this many eigenpairs converge + * MaxIter : maximum outer (restart) iterations + * Tolerance: relative convergence tolerance + * + * Usage + * ----- + * HarmonicBlockKrylovSchur hbks(LinOp, Grid, tol, shift, EvalNormSmall); + * std::vector v0(Nblock, Field(Grid)); + * // fill v0 with random starting vectors + * hbks(v0, maxIter, Nm, Nk, Nstop, Nblock); + * auto evals = hbks.getEvals(); + * auto evecs = hbks.getEvecs(); + */ +template +class HarmonicBlockKrylovSchur { + + typedef Eigen::MatrixXcd CMat; + typedef Eigen::VectorXcd CVec; + + //-------------------------------------------------------------------- + // Parameters + //-------------------------------------------------------------------- + int Nblock; + int Nm; + int Nk; + int Nstop; + int MaxIter; + RealD Tolerance; + ComplexD shift; // target shift σ + + //-------------------------------------------------------------------- + // Internal state + //-------------------------------------------------------------------- + LinearOperatorBase& Linop; + GridBase* Grid_; + RitzFilter ritzFilter; + + std::vector basis; // Nm*Nblock flat basis + CMat H; // (Nm*Nblock)² block-Hessenberg Rayleigh quotient + std::vector F; // Nblock residual vectors + CMat B; // (Nm*Nblock) × Nblock coupling matrix + RealD beta_k; + RealD rtol; + + CVec evals; + CMat littleEvecs; + std::vector ritzEstimates; + +public: + std::vector evecs; + + //-------------------------------------------------------------------- + // Constructor + //-------------------------------------------------------------------- + HarmonicBlockKrylovSchur(LinearOperatorBase& _Linop, GridBase* _Grid, + RealD _Tolerance, ComplexD _shift = 0.0, + RitzFilter _rf = EvalNormSmall) + : Linop(_Linop), Grid_(_Grid), Tolerance(_Tolerance), shift(_shift), + ritzFilter(_rf), + Nblock(-1), Nm(-1), Nk(-1), Nstop(-1), MaxIter(-1), + beta_k(0.0), rtol(0.0) + {} + + //-------------------------------------------------------------------- + // 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) + { + MaxIter = _maxIter; + Nm = _Nm; + Nk = _Nk; + Nstop = _Nstop; + Nblock = _Nblock; + + assert((int)v0.size() >= Nblock); + assert(Nk < Nm); + + int N = Nm * Nblock; + + RealD approxLambdaMax = approxMaxEval(v0[0]); + rtol = Tolerance * approxLambdaMax; + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: 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(v0.begin(), v0.begin() + Nblock); + + for (int iter = 0; iter < MaxIter; iter++) { + 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; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi"; + 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.schurReorder(Nk * Nblock); + + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk*Nblock):" << std::endl; + CMat S = schur.getMatrixS(); + for (int i = 0; i < Nk * Nblock; i++) + std::cout << GridLogMessage << " [" << i << "] " << S(i, i) << std::endl; + + CMat Q = schur.getMatrixQ(); + CMat Qt = Q.adjoint(); + + // ---- Rotate Krylov basis using Q from Hhat ---- + std::vector basis2; + constructUR(basis2, basis, Qt, N); + basis = basis2; + + // ---- Update H and B (rotate H, not Hhat) ---- + H = Q * H * Qt; + B = Q * B; + + // ---- Truncate to Nk*Nblock ---- + int Nkeep = Nk * Nblock; + + 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; + + std::vector basisTmp(basis.begin(), basis.begin() + Nkeep); + basis = basisTmp; + + CMat Btmp = B(Eigen::seqN(0, Nkeep), Eigen::all); + B = CMat::Zero(N, Nblock); + B(Eigen::seqN(0, Nkeep), Eigen::all) = Btmp; + + beta_k = Btmp.norm(); + 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; + + if (doVerify) { + std::string lbl = "iter " + std::to_string(iter) + " after restart+truncation"; + verify(lbl); + } + + // ---- Eigensystem of truncated H for convergence ---- + CMat Hk = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep)); + computeEigensystem(Hk, Nkeep); + + int Nconv = converged(Nkeep); + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: converged " << Nconv + << " / " << Nstop << std::endl; + + if (Nconv >= Nstop || iter == MaxIter - 1) { + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: done after " << iter + << " restarts, " << Nconv << " converged." << std::endl; + std::cout << GridLogMessage << "Eigenvalues: " << evals.transpose() << std::endl; + return; + } + } + } + + // Accessors + std::vector getEvecs() { return evecs; } + CVec getEvals() { return evals; } + std::vector getRitzEstimates() { return ritzEstimates; } + + //-------------------------------------------------------------------- + // Verification: check A V = V H + F B^dag explicitly + //-------------------------------------------------------------------- + /** + * Checks the block Arnoldi / Krylov-Schur decomposition + * + * 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. + * + * Prints: + * - H (current Rayleigh quotient, nBasis × nBasis) + * - B (coupling matrix, nBasis × Nblock) + * - M (explicit inner product matrix ) + * - max |H[i,j] - M[i,j]| (should be O(machine epsilon)) + * - max | - delta_ij| (orthonormality check) + * - for each basis column j: || A v_j - V H[:,j] - F B[j,:]^* || + * - max || (F orthogonal to basis) + */ + void verify(const std::string& label = "") + { + int nBasis = (int)basis.size(); + int nF = (int)F.size(); + + if (nBasis == 0) { + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur::verify [" << label + << "]: basis is empty." << std::endl; + return; + } + + std::cout << GridLogMessage + << "======== HarmonicBlockKrylovSchur::verify [" << label << "] ========" << std::endl; + std::cout << GridLogMessage + << " nBasis = " << nBasis << " Nblock = " << Nblock + << " nF = " << nF << std::endl; + + // ---- Print H ---- + 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::endl; + } + + // ---- Print B ---- + 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::endl; + } + + // ---- Compute M[i,j] = ---- + CMat M = CMat::Zero(nBasis, nBasis); + Field w(Grid_); + for (int j = 0; j < nBasis; j++) { + Linop.Op(basis[j], w); + for (int i = 0; i < nBasis; i++) + M(i, j) = toStdCmplx(innerProduct(basis[i], w)); + } + + 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::endl; + } + + // ---- max |H - M| ---- + RealD maxHM = 0.0; + for (int i = 0; i < nBasis; i++) + for (int j = 0; j < nBasis; j++) + maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j))); + std::cout << GridLogMessage + << " max |H[i,j] - M[i,j]| = " << maxHM << std::endl; + + // ---- Check orthonormality of basis ---- + CMat G = CMat::Zero(nBasis, nBasis); + for (int i = 0; i < nBasis; i++) + for (int j = 0; j < nBasis; j++) + G(i, j) = toStdCmplx(innerProduct(basis[i], basis[j])); + CMat Gerr = G - CMat::Identity(nBasis, nBasis); + std::cout << GridLogMessage + << " max | - delta_ij| = " << Gerr.cwiseAbs().maxCoeff() << std::endl; + + // ---- Per-column residual: || A v_j - V H[:,j] - F B[j,:]^* || ---- + RealD maxColDev = 0.0; + for (int j = 0; j < nBasis; j++) { + Linop.Op(basis[j], w); + + // subtract V H[:,j] + for (int i = 0; i < nBasis; i++) + w -= basis[i] * H(i, j); + + // subtract F B[j,:]^* (F[t] * conj(B[j,t])) + for (int t = 0; t < nF; t++) + w -= F[t] * std::conj(B(j, t)); + + RealD dev = std::sqrt(norm2(w)); + std::cout << GridLogMessage + << " || A v[" << j << "] - V H[:,j] - F B[j,:]* || = " << dev << std::endl; + maxColDev = std::max(maxColDev, dev); + } + std::cout << GridLogMessage + << " max column deviation = " << maxColDev << std::endl; + + // ---- Check F block orthogonality against basis ---- + if (nF > 0) { + RealD maxFV = 0.0; + for (int t = 0; t < nF; t++) + for (int i = 0; i < nBasis; i++) { + RealD ip = std::abs(toStdCmplx(innerProduct(basis[i], F[t]))); + maxFV = std::max(maxFV, ip); + } + std::cout << GridLogMessage + << " max || (should be ~0) = " << maxFV << std::endl; + } + + std::cout << GridLogMessage + << "======== end verify ========" << std::endl; + } + +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 + //-------------------------------------------------------------------- + void blockArnoldiIteration(std::vector& startBlock, int endBlock, + int startIdx, bool doubleOrthog) + { + int N = Nm * Nblock; + + if (startIdx == 0) { + basis.clear(); + F.clear(); + H = CMat::Zero(N, N); + B = CMat::Zero(N, Nblock); + + std::vector V0 = startBlock; + 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: + // + // A V_k = V_k S_k + F_old B_old^dag (*) + // + // 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. + int Nkeep = startIdx * Nblock; + for (auto& v : startBlock) basis.push_back(v); + + // Fill restart coupling rows into H + for (int t = 0; t < Nblock; t++) + for (int j = 0; j < Nkeep; j++) + H(Nkeep + t, j) = std::conj(B(j, t)); + + // Zero out B for the retained columns now that the coupling is in H + for (int j = 0; j < Nkeep; j++) + for (int t = 0; t < Nblock; t++) + B(j, t) = 0.0; + } + + for (int k = startIdx; k < endBlock; k++) + blockArnoldiStep(k, doubleOrthog); + } + + //-------------------------------------------------------------------- + // One block Arnoldi step + //-------------------------------------------------------------------- + void blockArnoldiStep(int k, bool doubleOrthog) + { + int kBase = k * Nblock; + int prevN = kBase + Nblock; + int N = Nm * Nblock; + + std::vector W(Nblock, Field(Grid_)); + for (int t = 0; t < Nblock; t++) + Linop.Op(basis[kBase + t], W[t]); + + // Full reorthogonalisation against all current basis vectors + for (int pass = 0; pass < (doubleOrthog ? 2 : 1); pass++) { + for (int i = 0; i < prevN; i++) { + for (int t = 0; t < Nblock; t++) { + ComplexD coeff = innerProduct(basis[i], W[t]); + if (pass == 0) + H(i, kBase + t) = toStdCmplx(coeff); + else + H(i, kBase + t) += toStdCmplx(coeff); + W[t] -= coeff * basis[i]; + } + } + } + + F = W; + + if (k == Nm - 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); + for (int t = 0; t < Nblock; t++) + 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; + } + + // Not last: QR the residual, extend basis + CMat R = blockQR(F); + + int nextBase = (k + 1) * Nblock; + for (int i = 0; i < Nblock; i++) + for (int j = 0; j < Nblock; j++) + H(nextBase + i, kBase + j) = R(i, j); + + 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)); + } + + //-------------------------------------------------------------------- + // Block QR (modified Gram-Schmidt within the block) + //-------------------------------------------------------------------- + CMat blockQR(std::vector& W) + { + CMat R = CMat::Zero(Nblock, Nblock); + const RealD deflThresh = 1e-14; + + for (int j = 0; j < Nblock; j++) { + for (int i = 0; i < j; i++) { + ComplexD coeff = innerProduct(W[i], W[j]); + R(i, j) = toStdCmplx(coeff); + W[j] -= coeff * W[i]; + } + RealD nrm = std::sqrt(norm2(W[j])); + R(j, j) = nrm; + if (nrm > deflThresh) { + W[j] *= (1.0 / nrm); + } else { + W[j] = Zero(); + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: deflation at block column " << j + << " (norm = " << nrm << ")" << std::endl; + } + } + return R; + } + + //-------------------------------------------------------------------- + // Orthonormalise starting block + //-------------------------------------------------------------------- + void blockOrthonormalise(std::vector& V) + { + for (int j = 0; j < (int)V.size(); j++) { + for (int i = 0; i < j; i++) { + ComplexD c = innerProduct(V[i], V[j]); + V[j] -= c * V[i]; + } + RealD nrm = std::sqrt(norm2(V[j])); + assert(nrm > 1e-14); + V[j] *= (1.0 / nrm); + } + } + + //-------------------------------------------------------------------- + // Basis rotation: UR[i] = sum_j U[j] * R(j,i) + //-------------------------------------------------------------------- + void constructUR(std::vector& UR, std::vector& U, + CMat& R, int N) + { + UR.clear(); + Field tmp(Grid_); + for (int i = 0; i < N; i++) { + tmp = Zero(); + for (int j = 0; j < N; j++) + tmp += U[j] * R(j, i); + UR.push_back(tmp); + } + } + + //-------------------------------------------------------------------- + // Eigensystem of the truncated H (not Hhat) + //-------------------------------------------------------------------- + /** + * Eigenvalues of H_k are the standard Ritz values in the retained + * subspace. After convergence has been declared via harmonic estimates, + * the final reported eigenvalues and vectors come from H_k (not Hhat_k), + * since H_k contains the true projected operator. + */ + void computeEigensystem(CMat& Hk, int Nkeep) + { + Eigen::ComplexEigenSolver es; + es.compute(Hk); + evals = es.eigenvalues(); + littleEvecs = es.eigenvectors(); + + evecs.clear(); + for (int k = 0; k < Nkeep; k++) { + CVec vec = littleEvecs.col(k); + Field tmp(Grid_); + tmp = Zero(); + for (int j = 0; j < (int)basis.size() && j < Nkeep; j++) + tmp += vec[j] * basis[j]; + evecs.push_back(tmp); + } + } + + //-------------------------------------------------------------------- + // Convergence check + //-------------------------------------------------------------------- + /** + * Ritz estimate for eigenpair k: || B^H y_k || + * where y_k is the k-th eigenvector of the truncated H. + * The same bound applies whether using Ritz or harmonic Ritz restart. + */ + int converged(int Nkeep) + { + ritzEstimates.clear(); + int Nconv = 0; + + CMat Bk = B(Eigen::seqN(0, Nkeep), Eigen::all); + + for (int k = 0; k < Nkeep; k++) { + CVec yk = littleEvecs.col(k); + CVec Bty = Bk.adjoint() * yk; + RealD res = Bty.norm(); + ritzEstimates.push_back(res); + std::cout << GridLogMessage + << "HarmonicBlockKrylovSchur: Ritz estimate[" << k + << "] = " << res << " eval = " << evals[k] << std::endl; + if (res < rtol) Nconv++; + } + return Nconv; + } + + //-------------------------------------------------------------------- + // Approximate maximum eigenvalue (power iteration) + //-------------------------------------------------------------------- + RealD approxMaxEval(const Field& v0, int MAX_ITER = 50) + { + assert(norm2(v0) > 1e-8); + RealD lam = 0.0, denom = std::sqrt(norm2(v0)); + Field vcur(Grid_), vtmp(Grid_); + vcur = v0; + for (int i = 0; i < MAX_ITER; i++) { + Linop.Op(vcur, vtmp); + vcur = vtmp; + RealD num = std::sqrt(norm2(vcur)); + lam = num / denom; + denom = num; + } + return lam; + } + +}; + +NAMESPACE_END(Grid); + +#endif // GRID_HARMONIC_BLOCKED_KRYLOV_SCHUR_H diff --git a/examples/Example_krylov_schur.cc b/examples/Example_krylov_schur.cc index 48967fda1..078588975 100644 --- a/examples/Example_krylov_schur.cc +++ b/examples/Example_krylov_schur.cc @@ -52,6 +52,8 @@ struct LanczosParameters: Serializable { Integer, Np, Integer, ReadEvec, Integer, maxIter, + Integer, Nblock, + Integer, verify, RealD, resid, RealD, ChebyLow, RealD, ChebyHigh, @@ -333,18 +335,27 @@ int main (int argc, char ** argv) // KrySchur(src, maxIter, Nm, Nk, Nstop); // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // Hacked, really EvalImagSmall -#if 1 RealD shift=1.5; +#if 0 KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); #else - KrylovSchur KrySchur (Iwilson, UGrid, resid,EvalImNormSmall); - KrySchur(src[0], maxIter, Nm, Nk, Nstop); + int Nblock=4; + 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,Nblock,true,true); +// BlockedKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); +// KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); + HarmonicBlockedKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall); + KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true); #endif std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl; src[0]=KrySchur.evecs[0]; - for (int i=1;i + +using namespace std; +using namespace Grid; + +// Diagonal real Hermitian operator (eigenvalues = scale lattice sites) +template +class DumbOperator : public LinearOperatorBase { +public: + LatticeComplex scale; + + DumbOperator(GridBase* grid) : scale(grid) { + GridParallelRNG pRNG(grid); + std::vector seeds({5,6,7,8}); + pRNG.SeedFixedIntegers(seeds); + random(pRNG, scale); + scale = exp(-Grid::real(scale) * 3.0); + } + + void OpDirAll(const Field& in, std::vector& out) {} + void OpDiag(const Field& in, Field& out) {} + void OpDir(const Field& in, Field& out, int dir, int disp) {} + + void Op(const Field& in, Field& out) { out = scale * in; } + void AdjOp(const Field& in, Field& out) { out = scale * in; } + void HermOp(const Field& in, Field& out) { out = scale * in; } + void HermOpAndNorm(const Field& in, Field& out, double& n1, double& n2) { + out = scale * in; + ComplexD d = innerProduct(in, out); n1 = real(d); + d = innerProduct(out, out); n2 = real(d); + } +}; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), + GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + + GridParallelRNG RNG(grid); + RNG.SeedFixedIntegers({1,2,3,4}); + + typedef LatticeComplex Field; + + DumbOperator op(grid); + + int nFail = 0; + + //-------------------------------------------------------------------- + // Helper lambda: run BKS with doVerify and check it doesn't crash + //-------------------------------------------------------------------- + auto runTest = [&](const std::string& label, int Nblock, int Nm, int Nk, + int maxIter, int Nstop) { + std::cout << GridLogMessage << "===== " << label << " =====" << std::endl; + + BlockedKrylovSchur bks(op, grid, 1e-6, EvalReSmall); + + std::vector v0(Nblock, Field(grid)); + for (int t = 0; t < Nblock; t++) random(RNG, v0[t]); + + bks(v0, maxIter, Nm, Nk, Nstop, Nblock, + /*doubleOrthog=*/true, /*doVerify=*/true); + + std::cout << GridLogMessage << label << " done." << std::endl; + }; + + // Test 1: Nblock=1 — scalar case, regression + 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); + + // 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); + + // Test 4: Nblock=2, larger cycle — more restarts + runTest("Nblock=2 Nm=12 Nk=6 maxIter=5", 2, 12, 6, 5, 6); + + if (nFail == 0) + std::cout << GridLogMessage << "All BlockedKrylovSchur tests completed." << std::endl; + + Grid_finalize(); + return nFail; +} diff --git a/tests/lanczos/Test_harmonic_vs_krylov_schur.cc b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc new file mode 100644 index 000000000..98bd3b12c --- /dev/null +++ b/tests/lanczos/Test_harmonic_vs_krylov_schur.cc @@ -0,0 +1,154 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Comparison test: HarmonicBlockedKrylovSchur vs BlockedKrylovSchur. + + Both algorithms are run on the same diagonal Hermitian operator and the + resulting eigenvalues are compared. doVerify=true is used so the KS + decomposition check max|H-M| and the per-column residuals are printed + at each step. For BKS these should be O(machine epsilon) at all times. + For HBKS they should be O(machine epsilon) AFTER Arnoldi, but may show + large per-column deviations AFTER restart+truncation (because the rotation + Q from Schur(Hhat) does not give an upper-triangular H_new, so the + truncated KS relation is only approximate). + +*************************************************************************************/ +#include + +using namespace std; +using namespace Grid; + +// Diagonal real Hermitian operator out = scale * in +template +class DumbOperator : public LinearOperatorBase { +public: + LatticeComplex scale; + DumbOperator(GridBase* grid) : scale(grid) { + GridParallelRNG pRNG(grid); + pRNG.SeedFixedIntegers({5,6,7,8}); + random(pRNG, scale); + scale = exp(-Grid::real(scale) * 3.0); + } + void OpDirAll(const Field& in, std::vector& out) {} + void OpDiag(const Field& in, Field& out) {} + void OpDir(const Field& in, Field& out, int dir, int disp) {} + void Op(const Field& in, Field& out) { out = scale * in; } + void AdjOp(const Field& in, Field& out) { out = scale * in; } + void HermOp(const Field& in, Field& out) { out = scale * in; } + void HermOpAndNorm(const Field& in, Field& out, double& n1, double& n2) { + out = scale * in; + ComplexD d = innerProduct(in, out); n1 = real(d); + d = innerProduct(out, out); n2 = real(d); + } +}; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), + GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + + GridParallelRNG RNG(grid); + RNG.SeedFixedIntegers({1,2,3,4}); + + typedef LatticeComplex Field; + DumbOperator op(grid); + + //---------------------------------------------------------------------- + // Parameters (kept small so output is readable) + //---------------------------------------------------------------------- + const int Nblock = 2; + const int Nm = 6; + const int Nk = 3; + const int Nstop = 2; + const int maxIter = 4; + const RealD tol = 1e-6; + + // Two identical starting blocks + std::vector v0(Nblock, Field(grid)); + std::vector v0b(Nblock, Field(grid)); + for (int t = 0; t < Nblock; t++) { + random(RNG, v0[t]); + v0b[t] = v0[t]; + } + + //---------------------------------------------------------------------- + // Run BlockedKrylovSchur with doVerify=true + //---------------------------------------------------------------------- + std::cout << GridLogMessage + << "\n========================================" << std::endl; + std::cout << GridLogMessage + << " BlockedKrylovSchur (Nblock=" << Nblock + << " Nm=" << Nm << " Nk=" << Nk << ")" << std::endl; + std::cout << GridLogMessage + << "========================================\n" << std::endl; + + BlockedKrylovSchur bks(op, grid, tol, EvalReSmall); + bks(v0, maxIter, Nm, Nk, Nstop, Nblock, + /*doubleOrthog=*/true, /*doVerify=*/true); + + auto bks_evals = bks.getEvals(); + std::cout << GridLogMessage + << "BKS eigenvalues (" << bks_evals.size() << "):" << std::endl; + for (int k = 0; k < (int)bks_evals.size(); k++) + std::cout << GridLogMessage << " [" << k << "] " << bks_evals[k] << std::endl; + + //---------------------------------------------------------------------- + // Run HarmonicBlockedKrylovSchur with doVerify=true + //---------------------------------------------------------------------- + std::cout << GridLogMessage + << "\n========================================" << std::endl; + std::cout << GridLogMessage + << " HarmonicBlockedKrylovSchur (Nblock=" << Nblock + << " Nm=" << Nm << " Nk=" << Nk << " shift=0)" << std::endl; + std::cout << GridLogMessage + << "========================================\n" << std::endl; + + HarmonicBlockedKrylovSchur hbks(op, grid, tol, 0.0, EvalNormSmall); + hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock, + /*doubleOrthog=*/true, /*doVerify=*/true); + + auto hbks_evals = hbks.getEvals(); + std::cout << GridLogMessage + << "HBKS eigenvalues (" << hbks_evals.size() << "):" << std::endl; + for (int k = 0; k < (int)hbks_evals.size(); k++) + std::cout << GridLogMessage << " [" << k << "] " << hbks_evals[k] << std::endl; + + //---------------------------------------------------------------------- + // Compare + //---------------------------------------------------------------------- + std::cout << GridLogMessage + << "\n========================================" << std::endl; + std::cout << GridLogMessage << " Eigenvalue comparison" << std::endl; + std::cout << GridLogMessage + << "========================================" << std::endl; + + // Sort both sets by real part for comparison + std::vector bvec(bks_evals.data(), + bks_evals.data() + bks_evals.size()); + std::vector hvec(hbks_evals.data(), + hbks_evals.data() + hbks_evals.size()); + auto cmpRe = [](const ComplexD& a, const ComplexD& b){ return a.real() < b.real(); }; + std::sort(bvec.begin(), bvec.end(), cmpRe); + std::sort(hvec.begin(), hvec.end(), cmpRe); + + int nCmp = std::min(bvec.size(), hvec.size()); + double maxDiff = 0.0; + for (int k = 0; k < nCmp; k++) { + double diff = std::abs(bvec[k].real() - hvec[k].real()) + std::abs(bvec[k].imag() - hvec[k].imag()); + maxDiff = std::max(maxDiff, diff); + std::cout << GridLogMessage + << " k=" << k + << " BKS=" << bvec[k] + << " HBKS=" << hvec[k] + << " |diff|=" << diff << std::endl; + } + std::cout << GridLogMessage << " max |BKS - HBKS| = " << maxDiff << std::endl; + + Grid_finalize(); + return 0; +}