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

Adding threshold in Plaquette testing fft5d bugfix

Browse Source
This commit is contained in:
Chulwoo Jung
2026-05-08 18:53:30 -04:00
parent d5cab295b2
commit bb99123e6d
3 changed files with 276 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Grid/algorithms/iterative/Gamma5ScalarLanczos.h
+230
View File
@@ -0,0 +1,230 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
γ5-Scalar Lanczos algorithm for γ5-Hermitian operators (Wilson Dirac).
Single-vector Lanczos using the γ5-inner product (u,v)_{γ5} = u†γ5v.
Builds the Krylov space {q, D_W q, D_W² q, ...} and projects D_W onto it,
producing a real tridiagonal matrix T_m whose eigenvalues approximate those
of D_W. Complex conjugate pairs of D_W eigenvalues appear as complex
eigenvalues of the real non-symmetric T_m.
Recurrence (G_k = sign(q_k†γ5 q_k) = ±1):
r_k = D_W q_k − α_k q_k − β_sup[k-1] q_{k-1}
β_sub[k] = √|r_k†γ5 r_k| (lower subdiag T_{k+1,k})
β_sup[k] = β_sub[k] G_{k+1}/G_k (upper superdiag T_{k,k+1})
q_{k+1} = r_k / β_sub[k]
T_m is real but not symmetric in general; eigenvalues via EigenSolver.
*************************************************************************************/
#ifndef GRID_GAMMA5_SCALAR_LANCZOS_H
#define GRID_GAMMA5_SCALAR_LANCZOS_H
#include <functional>
#include <numeric>
#include <iomanip>
NAMESPACE_BEGIN(Grid);
template<class Field>
class Gamma5ScalarLanczos {
public:
using Gamma5Func = std::function<void(const Field&, Field&)>;
private:
typedef Eigen::MatrixXcd CMat;
typedef Eigen::VectorXcd CVec;
typedef Eigen::MatrixXd RMat;
LinearOperatorBase<Field>& Linop;
GridBase* Grid_;
Gamma5Func applyGamma5;
RealD Tolerance;
int nSteps;
std::vector<Field> basis;
std::vector<RealD> alpha_; // diagonal of T_m
std::vector<RealD> beta_sub_; // lower subdiagonal: T_{k+1,k}
std::vector<RealD> beta_sup_; // upper superdiagonal: T_{k,k+1}
std::vector<int> gsign_; // G_k = ±1
CVec evals_;
std::vector<Field> evecs_;
std::vector<RealD> residuals_;
public:
bool doEvalCheck = false;
Gamma5ScalarLanczos(LinearOperatorBase<Field>& op, GridBase* grid,
Gamma5Func g5, RealD tol = 1e-8)
: Linop(op), Grid_(grid), applyGamma5(g5), Tolerance(tol), nSteps(0)
{}
CVec getEvals() { return evals_; }
std::vector<Field> getEvecs() { return evecs_; }
std::vector<RealD> getResiduals() { return residuals_; }
// v0 : starting vector
// maxSteps : maximum Lanczos steps
// Nstop : target Ritz pairs to return (0 = all)
// reorthog : full γ5-reorthogonalisation at each step
// filter : eigenvalue selection criterion
void operator()(const Field& v0, int maxSteps, int Nstop,
bool reorthog = false, RitzFilter filter = EvalImNormSmall)
{
basis.clear(); alpha_.clear(); beta_sub_.clear(); beta_sup_.clear(); gsign_.clear();
nSteps = 0;
// Normalise starting vector in γ5-norm: q_0 = v0 / √|v0†γ5 v0|
Field q(Grid_); q = v0;
Field g5q(Grid_);
applyGamma5(q, g5q);
RealD omega = std::real(toStdCmplx(innerProduct(g5q, q)));
RealD beta0 = std::sqrt(std::abs(omega));
assert(beta0 > 1e-14 && "starting vector has zero γ5-norm");
q *= (1.0 / beta0);
gsign_.push_back(omega >= 0 ? +1 : -1);
basis.push_back(q);
std::cout << GridLogMessage
<< "Gamma5ScalarLanczos: start (v0†γ5 v0)=" << omega
<< " G[0]=" << gsign_[0] << std::endl;
for (int k = 0; k < maxSteps; k++) {
// ── Apply operator ─────────────────────────────────────────────────
Field p(Grid_);
Linop.Op(basis[k], p);
// ── α_k = (q_k†γ5 D_W q_k) / G_k (real) ─────────────────────────
applyGamma5(basis[k], g5q);
RealD a = std::real(toStdCmplx(innerProduct(g5q, p))) / RealD(gsign_[k]);
alpha_.push_back(a);
// ── r = D_W q_k − α_k q_k − β_sup[k-1] q_{k-1} ───────────────────
Field r(Grid_); r = p;
r -= basis[k] * a;
if (k > 0)
r -= basis[k-1] * beta_sup_[k-1];
// ── Optional γ5-reorthogonalisation ────────────────────────────────
if (reorthog) {
for (int i = 0; i <= k; i++) {
applyGamma5(basis[i], g5q);
ComplexD c = toStdCmplx(innerProduct(g5q, r)) / RealD(gsign_[i]);
r -= basis[i] * c;
}
}
// ── β_sub[k] = √|r†γ5 r| ───────────────────────────────────────────
applyGamma5(r, g5q);
omega = std::real(toStdCmplx(innerProduct(g5q, r)));
RealD beta = std::sqrt(std::abs(omega));
std::cout << GridLogMessage
<< "Gamma5ScalarLanczos: k=" << k
<< " alpha=" << a
<< " beta=" << beta
<< " G[k]=" << gsign_[k]
<< " (r†γ5r)=" << omega << std::endl;
beta_sub_.push_back(beta);
if (beta < Tolerance) {
std::cout << GridLogMessage
<< "Gamma5ScalarLanczos: invariant subspace at step " << k << std::endl;
nSteps = k + 1;
break;
}
int G_next = (omega >= 0) ? +1 : -1;
gsign_.push_back(G_next);
// β_sup[k] = T_{k,k+1} = β_sub[k] * G_{k+1} / G_k
beta_sup_.push_back(beta * RealD(G_next) / RealD(gsign_[k]));
Field qnew(Grid_); qnew = r; qnew *= (1.0 / beta);
basis.push_back(qnew);
nSteps = k + 1;
}
if (nSteps == 0) return;
computeRitzPairs(nSteps, Nstop, filter);
}
private:
void computeRitzPairs(int m, int Nstop, RitzFilter filter)
{
// Assemble real tridiagonal T_m
RMat Tm = RMat::Zero(m, m);
for (int k = 0; k < m; k++) {
Tm(k, k) = alpha_[k];
if (k + 1 < m) {
Tm(k+1, k) = beta_sub_[k];
Tm(k, k+1) = beta_sup_[k];
}
}
std::cout << GridLogMessage
<< "Gamma5ScalarLanczos: T_m (" << m << "×" << m << "):" << std::endl;
for (int i = 0; i < m; i++) {
for (int j = 0; j < m; j++)
std::cout << " " << std::setprecision(8) << std::setw(16) << Tm(i,j);
std::cout << std::endl;
}
Eigen::EigenSolver<RMat> es(Tm);
CVec lambdas = es.eigenvalues();
CMat Y = es.eigenvectors();
// Sort by filter criterion
ComplexComparator cComp(filter);
std::vector<int> idx(m);
std::iota(idx.begin(), idx.end(), 0);
std::sort(idx.begin(), idx.end(), [&](int a, int b){
return cComp(toStdCmplx(lambdas(a)), toStdCmplx(lambdas(b)));
});
int Nout = (Nstop > 0 && Nstop < m) ? Nstop : m;
evals_.resize(Nout);
evecs_.clear(); evecs_.reserve(Nout);
residuals_.clear(); residuals_.reserve(Nout);
// Trailing β for residual estimate: β_sub[m-1] * |y_j[m-1]|
RealD beta_trail = (m - 1 < (int)beta_sub_.size()) ? beta_sub_[m-1] : 0.0;
for (int ji = 0; ji < Nout; ji++) {
int j = idx[ji];
evals_(ji) = lambdas(j);
CVec yj = Y.col(j);
// Ritz vector: u_j = Σ_k q_k y_j[k]
Field uj(Grid_); uj = Zero();
for (int k = 0; k < m && k < (int)basis.size(); k++)
uj += basis[k] * toStdCmplx(yj(k));
evecs_.push_back(uj);
RealD res = beta_trail * std::abs(toStdCmplx(yj(m-1)));
if (!std::isfinite(res)) res = std::numeric_limits<RealD>::infinity();
residuals_.push_back(res);
std::cout << GridLogMessage
<< " Ritz[" << std::setw(3) << ji << "]"
<< " lambda=" << evals_(ji)
<< " |res|=" << res;
if (doEvalCheck) {
Field w(Grid_);
Linop.Op(uj, w);
w -= uj * evals_(ji);
RealD actual = std::sqrt(norm2(w));
std::cout << " ||D_W u - λu||=" << actual;
}
std::cout << std::endl;
}
}
};
NAMESPACE_END(Grid);
#endif // GRID_GAMMA5_SCALAR_LANCZOS_H
examples/fft5d.cc
+16 -12
View File
@@ -237,36 +237,39 @@ static void trajFFT(const std::vector<LatticeComplexD>& in,
GridBase* g = in[0].Grid();
long lsites = g->lSites();
// Pack: buf[traj * lsites + lsite] (traj varies slowly, site varies fast)
std::vector<fftw_complex> ibuf((long)Ntraj * lsites);
// fftw_complex is double[2] — std::vector<double[2]> is not supported by nvcc.
// Use std::vector<double> with 2x size and reinterpret_cast.
std::vector<double> ibuf((long)Ntraj * lsites * 2, 0.0);
for (int n = 0; n < Ntraj; n++) {
std::vector<ComplexD> lc;
unvectorizeToLexOrdArray(lc, in[n]);
for (long s = 0; s < lsites; s++) {
ibuf[(long)n*lsites + s][0] = lc[s].real();
ibuf[(long)n*lsites + s][1] = lc[s].imag();
ibuf[((long)n*lsites + s)*2 + 0] = lc[s].real();
ibuf[((long)n*lsites + s)*2 + 1] = lc[s].imag();
}
}
// lsites transforms of length Ntraj, stride=lsites, dist=1
std::vector<fftw_complex> obuf((long)Ntraj * lsites);
std::vector<double> obuf((long)Ntraj * lsites * 2, 0.0);
fftw_complex* iptr = reinterpret_cast<fftw_complex*>(ibuf.data());
fftw_complex* optr = reinterpret_cast<fftw_complex*>(obuf.data());
int n1[1] = {Ntraj};
fftw_plan p = fftw_plan_many_dft(
1, n1, (int)lsites,
ibuf.data(), nullptr, (int)lsites, 1,
obuf.data(), nullptr, (int)lsites, 1,
iptr, nullptr, (int)lsites, 1,
optr, nullptr, (int)lsites, 1,
FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
// vector::assign triggers _M_fill_assign which needs a default constructor;
// Lattice has none. Use explicit push_back instead.
// Lattice has none. Use explicit emplace_back instead.
out.clear(); out.reserve(Ntraj);
for (int k = 0; k < Ntraj; k++) out.emplace_back(g);
for (int k = 0; k < Ntraj; k++) {
std::vector<ComplexD> lc(lsites);
for (long s = 0; s < lsites; s++)
lc[s] = ComplexD(obuf[(long)k*lsites+s][0], obuf[(long)k*lsites+s][1]);
lc[s] = ComplexD(obuf[((long)k*lsites+s)*2], obuf[((long)k*lsites+s)*2+1]);
vectorizeFromLexOrdArray(lc, out[k]);
}
}
@@ -311,10 +314,11 @@ static void analyzeTraj(const std::vector<LatticeComplexD>& traj,
}
// Autocorrelation via IFFT of the power spectrum
std::vector<fftw_complex> Pc(Nf);
for (int k = 0; k < Nf; k++) { Pc[k][0] = Pavg[k]; Pc[k][1] = 0.0; }
std::vector<double> Pc_buf(Nf * 2, 0.0);
for (int k = 0; k < Nf; k++) { Pc_buf[k*2] = Pavg[k]; Pc_buf[k*2+1] = 0.0; }
fftw_complex* Pc = reinterpret_cast<fftw_complex*>(Pc_buf.data());
std::vector<double> acorr(Ntraj, 0.0);
fftw_plan ip = fftw_plan_dft_c2r(1, &Ntraj, Pc.data(), acorr.data(), FFTW_ESTIMATE);
fftw_plan ip = fftw_plan_dft_c2r(1, &Ntraj, Pc, acorr.data(), FFTW_ESTIMATE);
fftw_execute(ip); fftw_destroy_plan(ip);
double c0 = acorr[0] / Ntraj;
{
tests/core/Test_plaquette_stats.cc
+28
View File
@@ -34,14 +34,17 @@ See the full license in the file "LICENSE" in the top level distribution directo
*
* Usage:
* ./Test_plaquette_stats [Grid options] [--file <nersc_config>] [--hot]
* [--threshold <val>]
*
* --file <path> Read gauge field from NERSC-format file
* --hot Use random (hot) SU(3) start (default: cold/unit start)
* --threshold <val> Print coordinates of every plaquette with Re Tr/Nc < val
*
* Grid size defaults to 8^3 x 16; override with --grid (e.g. --grid 4.4.4.8).
*/
#include <Grid/Grid.h>
#include <limits>
using namespace std;
using namespace Grid;
@@ -91,6 +94,16 @@ int main(int argc, char** argv)
if (std::string(argv[i]) == "--hot") { doHot = true; break; }
}
double threshold = std::numeric_limits<double>::quiet_NaN();
bool doThresh = false;
for (int i = 1; i < argc - 1; i++) {
if (std::string(argv[i]) == "--threshold") {
threshold = std::stod(argv[i+1]);
doThresh = true;
break;
}
}
if (!config_file.empty()) {
std::cout << GridLogMessage << "Reading gauge field from " << config_file << std::endl;
FieldMetaData header;
@@ -175,6 +188,21 @@ int main(int argc, char** argv)
<< std::setw(20) << std::setprecision(10) << local_min
<< std::setw(20) << std::setprecision(10) << avg
<< std::endl;
if (doThresh) {
for (int s = 0; s < (int)sv.size(); s++) {
RealD val = TensorRemove(sv[s]).real() / Nc;
if (val < threshold) {
Coordinate lc(Nd), gc(Nd);
Lexicographic::CoorFromIndex(lc, s, grid._ldimensions);
for (int d = 0; d < Nd; d++)
gc[d] = grid._processor_coor[d] * grid._ldimensions[d] + lc[d];
std::cout << "BELOW_THRESHOLD plane=" << planeName(mu, nu)
<< " site=(" << gc[0] << "," << gc[1] << "," << gc[2] << "," << gc[3] << ")"
<< " P=" << std::setprecision(10) << val << std::endl;
}
}
}
}
}