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

More debugging outputs

Browse Source
This commit is contained in:
Chulwoo Jung
2026-05-01 11:35:45 -04:00
parent 64ee6e90c9
commit 5754c965d5
2 changed files with 46 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/iterative/Gamma5BlockLanczos.h
+42 -2
View File
@@ -112,6 +112,9 @@ public:
assert(nrm > 1e-14);
v *= (1.0 / nrm);
applyGamma5(v, g5v);
// nrm = std::sqrt(norm2(g5v));
// assert(nrm > 1e-14);
// g5v *= (1.0 / nrm);
CMat2 G1 = gramMatrix(v, g5v);
ComplexD detG1 = G1(0,0)*G1(1,1) - G1(0,1)*G1(1,0);
@@ -400,7 +403,9 @@ public:
H_hess.block(0, 0, Nk, Nk) = H_L2;
// Trim basis to Nk (Q_{m+1} replaced by F_vec below)
basis.resize(Nk);
// Use erase instead of resize: resize instantiates _M_default_append
// which requires Lattice's default ctor (nonexistent).
basis.erase(basis.begin() + Nk, basis.end());
first_iter = false;
} else {
@@ -439,7 +444,7 @@ public:
H_hess.block(0, 0, Nk, Nk) = schur.getMatrixS().block(0, 0, Nk, Nk);
// Trim basis back to Nk (extension will re-grow it)
basis.resize(Nk);
basis.erase(basis.begin() + Nk, basis.end());
}
// ── L2-Arnoldi extension: add Np = Nmax-Nk vectors from F_vec ────────
@@ -1164,6 +1169,11 @@ private:
CMat2 Bk = B_blocks[step - 1];
Ck = invert2x2(Gkm1) * Bk.adjoint() * Gk;
}
ComplexD detG1 = Ck(0,0)*Ck(1,1) - Ck(0,1)*Ck(1,0);
std::cout << GridLogMessage
<< "Gamma5BlockLanczos: C "<<step<<"= \n" << Ck
<< "\n det = " << detG1 << std::endl;
C_blocks.push_back(Ck);
// (v) Residual R̂_k = P_k - Q_k A_k - Q_{k-1} C_k
@@ -1245,8 +1255,24 @@ private:
qnew1 = (r1 * U(0,0) + r2 * U(1,0)) * (1.0 / sqd0);
qnew2 = (r1 * U(0,1) + r2 * U(1,1)) * (1.0 / sqd1);
detG1 = Gkp1(0,0)*Gkp1(1,1) - Gkp1(0,1)*Gkp1(1,0);
std::cout << GridLogMessage
<< "Gamma5BlockLanczos: G "<<step<<"= \n" << Gkp1
<< "\n det G = " << detG1 << std::endl;
G_blocks.push_back(Gkp1);
detG1 = Bkp1(0,0)*Bkp1(1,1) - Bkp1(0,1)*Bkp1(1,0);
std::cout << GridLogMessage
<< "Gamma5BlockLanczos: B "<<step<<"= \n" << Bkp1
<< "\n det B = " << detG1 << std::endl;
B_blocks.push_back(Bkp1);
CMat2 G1 = gramMatrix(basis[0], basis[1],qnew1,qnew2);
detG1 = G1(0,0)*G1(1,1) - G1(0,1)*G1(1,0);
RealD absdetG1 = std::sqrt(detG1.real()*detG1.real() + detG1.imag()*detG1.imag());
std::cout << GridLogMessage
<< "Gamma5BlockLanczos: eta "<<step<<" = \n" << G1
<< "\n det = " << detG1 << std::endl;
basis.push_back(qnew1);
basis.push_back(qnew2);
@@ -1351,6 +1377,20 @@ private:
return G;
}
// γ5-Gram matrix of block [q1, q2]: G[i,j] = q^(i)† γ5 q^(j)
CMat2 gramMatrix(const Field& q11, const Field& q12, const Field& q21, const Field& q22)
{
Field g5q21(Grid_), g5q22(Grid_);
applyGamma5(q21, g5q21);
applyGamma5(q22, g5q22);
CMat2 G;
G(0,0) = toStdCmplx(innerProduct(q11, g5q21));
G(0,1) = toStdCmplx(innerProduct(q11, g5q22));
G(1,0) = toStdCmplx(innerProduct(q12, g5q21));
G(1,1) = toStdCmplx(innerProduct(q12, g5q22));
return G;
}
// M = Q†γ5 P for Q=[q1,q2], P=[p1,p2]: M[i,j] = q^(i)† γ5 p^(j)
CMat2 g5InnerBlock(const Field& q1, const Field& q2,
const Field& p1, const Field& p2)
examples/LanParams.xml
+4 -4
View File
@@ -5,11 +5,11 @@
<mstep>-0.025</mstep>
<M5>1.8</M5>
<Ls>48</Ls>
<Nstop>4</Nstop>
<Nk>16</Nk>
<Np>8</Np>
<Nstop>10</Nstop>
<Nk>10</Nk>
<Np>10</Np>
<ReadEvec>0</ReadEvec>
<maxIter>20</maxIter>
<maxIter>1000</maxIter>
<reorthog>1</reorthog>
<Nblock>4</Nblock>
<verify>1</verify>