portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-07-17 23:53:27 +01:00
Code Releases Activity

Adding (Harmonic) Block KS

Browse Source
This commit is contained in:
Chulwoo Jung
2026-03-31 12:36:18 -04:00
parent 09aa843984
commit 167f94e86c
6 changed files with 1673 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/Algorithms.h
+4
View File
@@ -85,8 +85,12 @@ NAMESPACE_CHECK(multigrid);
#include <Grid/algorithms/FFT.h> #include <Grid/algorithms/FFT.h>
#include <Grid/algorithms/iterative/KrylovSchur.h> #include <Grid/algorithms/iterative/KrylovSchur.h>
#include <Grid/algorithms/iterative/BlockedKrylovSchur.h>
#include <Grid/algorithms/iterative/HarmonicBlockedKrylovSchur.h>
#include <Grid/algorithms/iterative/Arnoldi.h> #include <Grid/algorithms/iterative/Arnoldi.h>
#include <Grid/algorithms/iterative/LanczosBidiagonalization.h> #include <Grid/algorithms/iterative/LanczosBidiagonalization.h>
#include <Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h> #include <Grid/algorithms/iterative/RestartedLanczosBidiagonalization.h>
#include <Grid/algorithms/iterative/GCR.h>
#include <Grid/algorithms/iterative/MultiSplittingPreconditionedCG.h>
#endif #endif
Grid/algorithms/iterative/BlockKrylovSchur.h
+701
View File
@@ -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 <paboyle@ph.ed.ac.uk>
Author: Chulwoo Jung <chulwoo@bnl.gov>
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 <iomanip>
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 Field>
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<Field>& 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<Field> 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<Field> 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<RealD> ritzEstimates;
public:
std::vector<Field> evecs;
//--------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------
BlockKrylovSchur(LinearOperatorBase<Field>& _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<Field>& 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<Field> 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<Field> 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<Field> 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<Field> getEvecs() { return evecs; }
CVec getEvals() { return evals; }
std::vector<RealD> 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] = <basis[i] | A basis[j]>
*
* and compared column-by-column against H. Separately, the nBasis × Nblock
* residual coupling matrix R is computed:
*
* R[j, t] = <basis[j] | F[t]> * ||F[t]|| (scaled by F-block norms)
*
* but since F is already normalised, R[j,t] = <basis[j] | F[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 <V | A V>)
* - 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] = <basis[i] | A basis[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 = <V|AV> (" << 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 |<V_i|V_j> - 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 |<V_i | F_t>| (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<Field>& 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<Field> 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<Field> 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] = <basis[i] | W[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<Field>& 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<Field>& 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<Field>& UR, std::vector<Field>& 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<CMat> 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
Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h
+700
View File
@@ -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 <paboyle@ph.ed.ac.uk>
Author: Chulwoo Jung <chulwoo@bnl.gov>
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 <iomanip>
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<Field> hbks(LinOp, Grid, tol, shift, EvalNormSmall);
* std::vector<Field> 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 Field>
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<Field>& Linop;
GridBase* Grid_;
RitzFilter ritzFilter;
std::vector<Field> basis; // Nm*Nblock flat basis
CMat H; // (Nm*Nblock)² block-Hessenberg Rayleigh quotient
std::vector<Field> F; // Nblock residual vectors
CMat B; // (Nm*Nblock) × Nblock coupling matrix
RealD beta_k;
RealD rtol;
CVec evals;
CMat littleEvecs;
std::vector<RealD> ritzEstimates;
public:
std::vector<Field> evecs;
//--------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------
HarmonicBlockKrylovSchur(LinearOperatorBase<Field>& _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<Field>& 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<Field> 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<Field> 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<Field> 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<Field> getEvecs() { return evecs; }
CVec getEvals() { return evals; }
std::vector<RealD> 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 <V | A V>)
* - max |H[i,j] - M[i,j]| (should be O(machine epsilon))
* - max |<V_i|V_j> - delta_ij| (orthonormality check)
* - for each basis column j: || A v_j - V H[:,j] - F B[j,:]^* ||
* - max |<V_i | F_t>| (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] = <basis[i] | A basis[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 = <V|AV> (" << 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 |<V_i|V_j> - 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 |<V_i | F_t>| (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<Field>& 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<Field> 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<Field> 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] = <basis[kBase+t] | A basis[j]>
// = conj(<basis[j] | A basis[kBase+t]>) = 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<Field>& 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<Field>& 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<Field>& UR, std::vector<Field>& 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<CMat> 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
examples/Example_krylov_schur.cc
+15 -4
View File
@@ -52,6 +52,8 @@ struct LanczosParameters: Serializable {
Integer, Np, Integer, Np,
Integer, ReadEvec, Integer, ReadEvec,
Integer, maxIter, Integer, maxIter,
Integer, Nblock,
Integer, verify,
RealD, resid, RealD, resid,
RealD, ChebyLow, RealD, ChebyLow,
RealD, ChebyHigh, RealD, ChebyHigh,
@@ -333,18 +335,27 @@ int main (int argc, char ** argv)
// KrySchur(src, maxIter, Nm, Nk, Nstop); // KrySchur(src, maxIter, Nm, Nk, Nstop);
// KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
// Hacked, really EvalImagSmall // Hacked, really EvalImagSmall
#if 1
RealD shift=1.5; RealD shift=1.5;
#if 0
KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
#else #else
KrylovSchur KrySchur (Iwilson, UGrid, resid,EvalImNormSmall); int Nblock=4;
KrySchur(src[0], maxIter, Nm, Nk, Nstop); 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 #endif
std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl; std::cout << GridLogMessage << "evec.size= " << KrySchur.evecs.size()<< std::endl;
src[0]=KrySchur.evecs[0]; src[0]=KrySchur.evecs[0];
for (int i=1;i<Nstop;i++) src[0]+=KrySchur.evecs[i]; std::cout << GridLogMessage << "KrySchur.evecs= "<< KrySchur.evecs.size() <<std::endl;
for (int i=1;i<Nk;i++) src[0]+=KrySchur.evecs[i];
for (int i=0;i<Nstop;i++) for (int i=0;i<Nstop;i++)
{ {
std::string evfile ("./evec_"+std::to_string(mass)+"_"+std::to_string(i)); std::string evfile ("./evec_"+std::to_string(mass)+"_"+std::to_string(i));
tests/lanczos/Test_blocked_krylov_schur.cc
+99
View File
@@ -0,0 +1,99 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Test for BlockedKrylovSchur: verifies the KS decomposition A V = V H + F B^dag
by explicit operator applies, before and after each restart.
Uses DumbOperator (diagonal, real, Hermitian) from Test_synthetic_lanczos.
Tests Nblock=1 (scalar, regression) and Nblock=2,3 (exercises the B^H fix).
*************************************************************************************/
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
// Diagonal real Hermitian operator (eigenvalues = scale lattice sites)
template<class Field>
class DumbOperator : public LinearOperatorBase<Field> {
public:
LatticeComplex scale;
DumbOperator(GridBase* grid) : scale(grid) {
GridParallelRNG pRNG(grid);
std::vector<int> 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<Field>& 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<Field> 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<Field> bks(op, grid, 1e-6, EvalReSmall);
std::vector<Field> 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;
}
tests/lanczos/Test_harmonic_vs_krylov_schur.cc
+154
View File
@@ -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 <Grid/Grid.h>
using namespace std;
using namespace Grid;
// Diagonal real Hermitian operator out = scale * in
template<class Field>
class DumbOperator : public LinearOperatorBase<Field> {
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<Field>& 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<Field> 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<Field> v0(Nblock, Field(grid));
std::vector<Field> 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<Field> 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<Field> 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<ComplexD> bvec(bks_evals.data(),
bks_evals.data() + bks_evals.size());
std::vector<ComplexD> 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;
}