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

Changing definition of Nm and Nk in BlockKrylovSchur

Browse Source
This commit is contained in:
Chulwoo Jung
2026-04-02 15:22:21 -04:00
parent 88ea24687f
commit 84707cc3a8
4 changed files with 62 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/iterative/BlockKrylovSchur.h
+31 -29
View File
@@ -44,10 +44,10 @@ NAMESPACE_BEGIN(Grid);
* 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
* V_k = Nm orthonormal basis vectors (stored flat in basis[])
* H_k = Nm x Nm 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)
* B_k = Nm 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
@@ -56,15 +56,15 @@ NAMESPACE_BEGIN(Grid);
*
* 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 leading Nk Schur vectors (chosen by RitzFilter) are retained,
* the basis and Rayleigh quotient are truncated, and block Arnoldi continues
* from the Nk-th block.
* from the (Nk/Nblock)-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)
* Nm : total Krylov dimension (must be divisible by Nblock)
* Nk : total vectors to keep after each restart (must be divisible by Nblock, 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|)
@@ -83,8 +83,8 @@ protected:
// Parameters (set by operator())
//--------------------------------------------------------------------
int Nblock; // block size
int Nm; // block steps (total dim = Nm * Nblock)
int Nk; // blocks retained after restart
int Nm; // total Krylov dimension (multiple of Nblock)
int Nk; // total vectors retained after restart (multiple of Nblock)
int Nstop;
int MaxIter;
RealD Tolerance;
@@ -97,19 +97,19 @@ protected:
RitzFilter ritzFilter;
// Flat storage: basis[s*Nblock + t] is the t-th vector of block s
// After construction: basis has Nm*Nblock entries
// After construction: basis has Nm entries
std::vector<Field> basis;
// Rayleigh quotient (Nm*Nblock) x (Nm*Nblock)
// Rayleigh quotient Nm x Nm
CMat H;
// Residual block: Nblock vectors (the (Nm+1)-th block, unnormalised before
// Residual block: Nblock vectors (the (Nm/Nblock+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.
// Coupling matrix B: Nm 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)
// B(Nm - Nblock + t, s) = H_{Nm/Nblock+1, Nm/Nblock}(t, s) (the subdiagonal block)
// We keep it as a full matrix for generality after restarts.
CMat B;
@@ -118,7 +118,7 @@ protected:
// Output
CVec evals;
CMat littleEvecs; // Nm*Nblock columns
CMat littleEvecs; // Nm columns
std::vector<RealD> ritzEstimates;
public:
@@ -144,8 +144,8 @@ public:
* ----------
* 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)
* _Nm : total Krylov dimension (must be divisible by _Nblock)
* _Nk : total vectors to keep after restart (must be divisible by _Nblock, _Nk < _Nm)
* _Nstop : stop after _Nstop eigenvalues converged
* _Nblock : block size
*/
@@ -160,10 +160,12 @@ public:
Nblock = _Nblock;
assert((int)v0.size() >= Nblock);
assert(Nm % Nblock == 0);
assert(Nk % Nblock == 0);
assert(Nk < Nm);
preRun(); // hook: derived classes add parameter assertions here
int N = Nm * Nblock; // total Krylov dimension
int N = Nm; // total Krylov dimension
// Approximate largest eigenvalue for tolerance normalisation
RealD approxLambdaMax = approxMaxEval(v0[0]);
@@ -181,11 +183,11 @@ public:
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);
// ---- Block Arnoldi: extend from block start to block Nm/Nblock ----
blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog);
// After first full cycle start from block Nk
start = Nk;
// After first full cycle start from block Nk/Nblock
start = Nk/Nblock;
if (doVerify) {
std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi";
@@ -196,8 +198,8 @@ public:
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);
// Reorder: bring wanted Nk Schur values to top-left
schur.schurReorder(Nk);
std::cout << GridLogMessage << "BlockKrylovSchur: Schur reordered." << std::endl;
CMat Q = schur.getMatrixQ();
@@ -212,8 +214,8 @@ public:
B = Q * B;
H = schur.getMatrixS();
// ---- Truncate to Nk*Nblock ----
int Nkeep = Nk * Nblock;
// ---- Truncate to Nk ----
int Nkeep = Nk;
CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
H = CMat::Zero(N, N);
@@ -449,7 +451,7 @@ private:
void blockArnoldiIteration(std::vector<Field>& startBlock, int endBlock,
int startIdx, bool doubleOrthog)
{
int N = Nm * Nblock;
int N = Nm;
if (startIdx == 0) {
basis.clear();
@@ -523,7 +525,7 @@ private:
{
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;
int N = Nm;
// W[t] = op(basis[kBase + t]) — dispatches to applyBlock() virtual
std::vector<Field> W(Nblock, Field(Grid_));
@@ -547,7 +549,7 @@ private:
// Store residual block F
F = W;
if (k == Nm - 1) {
if (k == Nm/Nblock - 1) {
// Last block: compute coupling matrix B for KS decomp.
//
// blockQR modifies F in-place (F → Q orthonormal) and returns R
Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h
+20 -18
View File
@@ -49,9 +49,9 @@ NAMESPACE_BEGIN(Grid);
*
* A V = V H + F B^dag (1)
*
* with V orthonormal (Nm*Nblock columns), H the (Nm*Nblock)² block
* with V orthonormal (Nm columns), H the Nm² block
* upper-Hessenberg Rayleigh quotient, F the Nblock residual vectors and B
* the (Nm*Nblock)×Nblock coupling matrix, the harmonic Rayleigh quotient
* the Nm×Nblock coupling matrix, the harmonic Rayleigh quotient
* relative to shift σ is
*
* Hhat = H + (H - σI)^{-H} B B^H (2)
@@ -94,8 +94,8 @@ NAMESPACE_BEGIN(Grid);
* ----------
* 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)
* Nm : total Krylov dimension (must be divisible by Nblock)
* Nk : total vectors to retain after each restart (must be divisible by Nblock, Nk < Nm)
* Nstop : stop when this many eigenpairs converge
* MaxIter : maximum outer (restart) iterations
* Tolerance: relative convergence tolerance
@@ -119,8 +119,8 @@ class HarmonicBlockKrylovSchur {
// Parameters
//--------------------------------------------------------------------
int Nblock;
int Nm;
int Nk;
int Nm; // total Krylov dimension (multiple of Nblock)
int Nk; // total vectors retained after restart (multiple of Nblock)
int Nstop;
int MaxIter;
RealD Tolerance;
@@ -173,9 +173,11 @@ public:
Nblock = _Nblock;
assert((int)v0.size() >= Nblock);
assert(Nm % Nblock == 0);
assert(Nk % Nblock == 0);
assert(Nk < Nm);
int N = Nm * Nblock;
int N = Nm;
RealD approxLambdaMax = approxMaxEval(v0[0]);
rtol = Tolerance * approxLambdaMax;
@@ -194,9 +196,9 @@ public:
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;
// ---- Block Arnoldi: extend from block 'start' to block Nm/Nblock ----
blockArnoldiIteration(startBlock, Nm/Nblock, start, doubleOrthog);
start = Nk/Nblock;
if (doVerify) {
std::string lbl = "iter " + std::to_string(iter) + " after Arnoldi";
@@ -209,12 +211,12 @@ public:
// ---- Schur decompose Hhat ----
ComplexSchurDecomposition schur(Hhat, false, ritzFilter);
schur.schurReorder(Nk * Nblock);
schur.schurReorder(Nk);
std::cout << GridLogMessage
<< "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk*Nblock):" << std::endl;
<< "HarmonicBlockKrylovSchur: harmonic Ritz values (first Nk):" << std::endl;
CMat S = schur.getMatrixS();
for (int i = 0; i < Nk * Nblock; i++)
for (int i = 0; i < Nk; i++)
std::cout << GridLogMessage << " [" << i << "] " << S(i, i) << std::endl;
CMat Q = schur.getMatrixQ();
@@ -229,8 +231,8 @@ public:
H = Q * H * Qt;
B = Q * B;
// ---- Truncate to Nk*Nblock ----
int Nkeep = Nk * Nblock;
// ---- Truncate to Nk ----
int Nkeep = Nk;
CMat Htmp = H(Eigen::seqN(0, Nkeep), Eigen::seqN(0, Nkeep));
H = CMat::Zero(N, N);
@@ -440,7 +442,7 @@ private:
void blockArnoldiIteration(std::vector<Field>& startBlock, int endBlock,
int startIdx, bool doubleOrthog)
{
int N = Nm * Nblock;
int N = Nm;
if (startIdx == 0) {
basis.clear();
@@ -492,7 +494,7 @@ private:
{
int kBase = k * Nblock;
int prevN = kBase + Nblock;
int N = Nm * Nblock;
int N = Nm;
std::vector<Field> W(Nblock, Field(Grid_));
for (int t = 0; t < Nblock; t++)
@@ -514,7 +516,7 @@ private:
F = W;
if (k == Nm - 1) {
if (k == Nm/Nblock - 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);
tests/lanczos/Test_blocked_krylov_schur.cc
+5 -5
View File
@@ -68,7 +68,7 @@ int main(int argc, char** argv)
int maxIter, int Nstop) {
std::cout << GridLogMessage << "===== " << label << " =====" << std::endl;
BlockedKrylovSchur<Field> bks(op, grid, 1e-6, EvalReSmall);
BlockKrylovSchur<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]);
@@ -79,17 +79,17 @@ int main(int argc, char** argv)
std::cout << GridLogMessage << label << " done." << std::endl;
};
// Test 1: Nblock=1 — scalar case, regression
// Test 1: Nblock=1 — scalar case, regression (Nm,Nk now total vectors)
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);
runTest("Nblock=2 Nm=16 Nk=8 maxIter=3", 2, 16, 8, 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);
runTest("Nblock=3 Nm=27 Nk=9 maxIter=3", 3, 27, 9, 3, 3);
// Test 4: Nblock=2, larger cycle — more restarts
runTest("Nblock=2 Nm=12 Nk=6 maxIter=5", 2, 12, 6, 5, 6);
runTest("Nblock=2 Nm=24 Nk=12 maxIter=5", 2, 24, 12, 5, 6);
if (nFail == 0)
std::cout << GridLogMessage << "All BlockedKrylovSchur tests completed." << std::endl;
tests/lanczos/Test_harmonic_vs_krylov_schur.cc
+6 -6
View File
@@ -62,8 +62,8 @@ int main(int argc, char** argv)
// Parameters (kept small so output is readable)
//----------------------------------------------------------------------
const int Nblock = 2;
const int Nm = 6;
const int Nk = 3;
const int Nm = 12; // total vectors (= 6 blocks * Nblock=2)
const int Nk = 6; // total kept (= 3 blocks * Nblock=2)
const int Nstop = 2;
const int maxIter = 4;
const RealD tol = 1e-6;
@@ -82,12 +82,12 @@ int main(int argc, char** argv)
std::cout << GridLogMessage
<< "\n========================================" << std::endl;
std::cout << GridLogMessage
<< " BlockedKrylovSchur (Nblock=" << Nblock
<< " BlockKrylovSchur (Nblock=" << Nblock
<< " Nm=" << Nm << " Nk=" << Nk << ")" << std::endl;
std::cout << GridLogMessage
<< "========================================\n" << std::endl;
BlockedKrylovSchur<Field> bks(op, grid, tol, EvalReSmall);
BlockKrylovSchur<Field> bks(op, grid, tol, EvalReSmall);
bks(v0, maxIter, Nm, Nk, Nstop, Nblock,
/*doubleOrthog=*/true, /*doVerify=*/true);
@@ -103,12 +103,12 @@ int main(int argc, char** argv)
std::cout << GridLogMessage
<< "\n========================================" << std::endl;
std::cout << GridLogMessage
<< " HarmonicBlockedKrylovSchur (Nblock=" << Nblock
<< " HarmonicBlockKrylovSchur (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);
HarmonicBlockKrylovSchur<Field> hbks(op, grid, tol, 0.0, EvalNormSmall);
hbks(v0b, maxIter, Nm, Nk, Nstop, Nblock,
/*doubleOrthog=*/true, /*doVerify=*/true);