From bb99123e6d4cd8e862872aaf793135389901c3cf Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Fri, 8 May 2026 18:53:30 -0400 Subject: [PATCH] Adding threshold in Plaquette testing fft5d bugfix --- .../iterative/Gamma5ScalarLanczos.h | 230 ++++++++++++++++++ examples/fft5d.cc | 28 ++- tests/core/Test_plaquette_stats.cc | 32 ++- 3 files changed, 276 insertions(+), 14 deletions(-) create mode 100644 Grid/algorithms/iterative/Gamma5ScalarLanczos.h diff --git a/Grid/algorithms/iterative/Gamma5ScalarLanczos.h b/Grid/algorithms/iterative/Gamma5ScalarLanczos.h new file mode 100644 index 000000000..b9795fed7 --- /dev/null +++ b/Grid/algorithms/iterative/Gamma5ScalarLanczos.h @@ -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 +#include +#include + +NAMESPACE_BEGIN(Grid); + +template +class Gamma5ScalarLanczos { +public: + using Gamma5Func = std::function; + +private: + typedef Eigen::MatrixXcd CMat; + typedef Eigen::VectorXcd CVec; + typedef Eigen::MatrixXd RMat; + + LinearOperatorBase& Linop; + GridBase* Grid_; + Gamma5Func applyGamma5; + RealD Tolerance; + int nSteps; + + std::vector basis; + std::vector alpha_; // diagonal of T_m + std::vector beta_sub_; // lower subdiagonal: T_{k+1,k} + std::vector beta_sup_; // upper superdiagonal: T_{k,k+1} + std::vector gsign_; // G_k = ±1 + + CVec evals_; + std::vector evecs_; + std::vector residuals_; + +public: + bool doEvalCheck = false; + + Gamma5ScalarLanczos(LinearOperatorBase& 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 getEvecs() { return evecs_; } + std::vector 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 es(Tm); + CVec lambdas = es.eigenvalues(); + CMat Y = es.eigenvectors(); + + // Sort by filter criterion + ComplexComparator cComp(filter); + std::vector 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::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 diff --git a/examples/fft5d.cc b/examples/fft5d.cc index abfa8f6fe..da1aeee1d 100644 --- a/examples/fft5d.cc +++ b/examples/fft5d.cc @@ -237,36 +237,39 @@ static void trajFFT(const std::vector& in, GridBase* g = in[0].Grid(); long lsites = g->lSites(); - // Pack: buf[traj * lsites + lsite] (traj varies slowly, site varies fast) - std::vector ibuf((long)Ntraj * lsites); + // fftw_complex is double[2] — std::vector is not supported by nvcc. + // Use std::vector with 2x size and reinterpret_cast. + std::vector ibuf((long)Ntraj * lsites * 2, 0.0); for (int n = 0; n < Ntraj; n++) { std::vector 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 obuf((long)Ntraj * lsites); + std::vector obuf((long)Ntraj * lsites * 2, 0.0); + fftw_complex* iptr = reinterpret_cast(ibuf.data()); + fftw_complex* optr = reinterpret_cast(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 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& traj, } // Autocorrelation via IFFT of the power spectrum - std::vector Pc(Nf); - for (int k = 0; k < Nf; k++) { Pc[k][0] = Pavg[k]; Pc[k][1] = 0.0; } + std::vector 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(Pc_buf.data()); std::vector 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; { diff --git a/tests/core/Test_plaquette_stats.cc b/tests/core/Test_plaquette_stats.cc index 2ba4c57be..a81cfc60c 100644 --- a/tests/core/Test_plaquette_stats.cc +++ b/tests/core/Test_plaquette_stats.cc @@ -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 ] [--hot] + * [--threshold ] * - * --file Read gauge field from NERSC-format file - * --hot Use random (hot) SU(3) start (default: cold/unit start) + * --file Read gauge field from NERSC-format file + * --hot Use random (hot) SU(3) start (default: cold/unit start) + * --threshold 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 +#include 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::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; + } + } + } } }