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

Expanding debug output after restart in block KS

Browse Source
This commit is contained in:
Chulwoo Jung
2026-04-06 12:00:33 -04:00
parent 84707cc3a8
commit d4b685cf2d
2 changed files with 37 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/iterative/BlockKrylovSchur.h
+32 -11
View File
@@ -314,23 +314,28 @@ public:
return; return;
} }
// Full allocated dimension of H (= Nm); may be larger than nBasis after restart+truncation
int Nfull = (int)H.rows();
std::cout << GridLogMessage std::cout << GridLogMessage
<< "======== BlockKrylovSchur::verify [" << label << "] ========" << std::endl; << "======== BlockKrylovSchur::verify [" << label << "] ========" << std::endl;
std::cout << GridLogMessage std::cout << GridLogMessage
<< " nBasis = " << nBasis << " Nblock = " << Nblock << " nBasis = " << nBasis << " Nfull = " << Nfull
<< " nF = " << nF << std::endl; << " Nblock = " << Nblock << " nF = " << nF << std::endl;
// ---- Print H ---- // ---- Print H (full Nm x Nm matrix) ----
std::cout << GridLogMessage << "H (" << nBasis << " x " << nBasis << "):" << std::endl; // After restart+truncation nBasis < Nfull; entries outside [0:nBasis,0:nBasis]
for (int i = 0; i < nBasis; i++) { // should be ~0 — printing the full matrix makes this visible.
for (int j = 0; j < nBasis; j++) std::cout << GridLogMessage << "H (" << Nfull << " x " << Nfull << "):" << std::endl;
for (int i = 0; i < Nfull; i++) {
for (int j = 0; j < Nfull; j++)
std::cout << " " << std::setw(14) << H(i, j); std::cout << " " << std::setw(14) << H(i, j);
std::cout << std::endl; std::cout << std::endl;
} }
// ---- Print B ---- // ---- Print B (full Nm x Nblock matrix) ----
std::cout << GridLogMessage << "B (" << nBasis << " x " << nF << "):" << std::endl; std::cout << GridLogMessage << "B (" << Nfull << " x " << nF << "):" << std::endl;
for (int i = 0; i < nBasis; i++) { for (int i = 0; i < Nfull; i++) {
for (int t = 0; t < nF; t++) for (int t = 0; t < nF; t++)
std::cout << " " << std::setw(14) << B(i, t); std::cout << " " << std::setw(14) << B(i, t);
std::cout << std::endl; std::cout << std::endl;
@@ -352,13 +357,29 @@ public:
std::cout << std::endl; std::cout << std::endl;
} }
// ---- max |H - M| ---- // ---- max |H - M| (for the nBasis x nBasis Rayleigh-quotient block) ----
RealD maxHM = 0.0; RealD maxHM = 0.0;
for (int i = 0; i < nBasis; i++) for (int i = 0; i < nBasis; i++)
for (int j = 0; j < nBasis; j++) for (int j = 0; j < nBasis; j++)
maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j))); maxHM = std::max(maxHM, std::abs(H(i,j) - M(i,j)));
std::cout << GridLogMessage std::cout << GridLogMessage
<< " max |H[i,j] - M[i,j]| = " << maxHM << std::endl; << " max |H[i,j] - M[i,j]| (i,j < nBasis) = " << maxHM << std::endl;
// ---- Check entries of H OUTSIDE the nBasis x nBasis block ----
// These live in rows i >= nBasis or columns j >= nBasis.
// They must be ~0 after restart+truncation (before the coupling rows are
// re-filled by the next blockArnoldiIteration call).
// During a full Arnoldi run (nBasis == Nfull) this loop is a no-op.
RealD maxOutside = 0.0;
for (int i = 0; i < Nfull; i++)
for (int j = 0; j < Nfull; j++)
if (i >= nBasis || j >= nBasis)
maxOutside = std::max(maxOutside, std::abs(H(i,j)));
if (Nfull > nBasis)
std::cout << GridLogMessage
<< " max |H[i,j]| outside [0:" << nBasis << ",0:" << nBasis
<< "] block = " << maxOutside
<< " (should be ~0 after restart+truncation)" << std::endl;
// ---- Check orthonormality of basis ---- // ---- Check orthonormality of basis ----
CMat G = CMat::Zero(nBasis, nBasis); CMat G = CMat::Zero(nBasis, nBasis);
examples/Example_krylov_schur.cc
+5 -5
View File
@@ -347,7 +347,7 @@ int main (int argc, char ** argv)
// KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall); // KrylovSchur KrySchur (HermOp2, UGrid, resid,EvalNormSmall);
// Hacked, really EvalImagSmall // Hacked, really EvalImagSmall
RealD shift=1.5; RealD shift=1.5;
#if 1 #if 0
KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
KrySchur(src[0], maxIter, Nm, Nk, Nstop); KrySchur(src[0], maxIter, Nm, Nk, Nstop);
// KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift); // KrySchur(src[0], maxIter, Nm, Nk, Nstop,&shift);
@@ -356,10 +356,10 @@ int main (int argc, char ** argv)
Nblock=LanParams.Nblock; Nblock=LanParams.Nblock;
bool if_verify=false; bool if_verify=false;
if(LanParams.verify) if_verify=true; if(LanParams.verify) if_verify=true;
KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); // KrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,true); // KrySchur(src, maxIter, Nm, Nk, Nstop,true,if_verify);
// BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall); BlockKrylovSchur KrySchur (Dwilson, UGrid, resid,EvalImNormSmall);
// KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify); KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true,if_verify);
// HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall); // HarmonicBlockKrylovSchur KrySchur (Dwilson, UGrid, resid,shift,EvalImNormSmall);
// KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true); // KrySchur(src, maxIter, Nm, Nk, Nstop,Nblock,true);
#endif #endif