From af1fb41d55921d64aa6fa7c5bac7086cff537b9b Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Wed, 8 Apr 2026 22:50:37 -0400 Subject: [PATCH] Precision reduction in Block KS verify() for better readability Adding Plaquette stat testing --- Grid/algorithms/iterative/BlockKrylovSchur.h | 6 +- .../iterative/HarmonicBlockKrylovSchur.h | 6 +- tests/core/Test_plaquette_stats.cc | 211 ++++++++++++++++++ 3 files changed, 217 insertions(+), 6 deletions(-) create mode 100644 tests/core/Test_plaquette_stats.cc diff --git a/Grid/algorithms/iterative/BlockKrylovSchur.h b/Grid/algorithms/iterative/BlockKrylovSchur.h index c3835c88d..c48c8b360 100644 --- a/Grid/algorithms/iterative/BlockKrylovSchur.h +++ b/Grid/algorithms/iterative/BlockKrylovSchur.h @@ -329,7 +329,7 @@ public: 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::setprecision(4) << std::setw(14) << H(i, j); std::cout << std::endl; } @@ -337,7 +337,7 @@ public: std::cout << GridLogMessage << "B (" << Nfull << " x " << nF << "):" << std::endl; for (int i = 0; i < Nfull; i++) { for (int t = 0; t < nF; t++) - std::cout << " " << std::setw(14) << B(i, t); + std::cout << " " << std::setprecision(4) << std::setw(14) << B(i, t); std::cout << std::endl; } @@ -353,7 +353,7 @@ public: std::cout << GridLogMessage << "M = (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << M(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << M(i, j); std::cout << std::endl; } diff --git a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h index 6b7daca67..117ef6a07 100644 --- a/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h +++ b/Grid/algorithms/iterative/HarmonicBlockKrylovSchur.h @@ -326,7 +326,7 @@ public: std::cout << GridLogMessage << "H (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << H(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << H(i, j); std::cout << std::endl; } @@ -334,7 +334,7 @@ public: std::cout << GridLogMessage << "B (" << nBasis << " x " << nF << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int t = 0; t < nF; t++) - std::cout << " " << std::setw(14) << B(i, t); + std::cout << " " << std::setprecision(4) << std::setw(14) << B(i, t); std::cout << std::endl; } @@ -350,7 +350,7 @@ public: std::cout << GridLogMessage << "M = (" << nBasis << " x " << nBasis << "):" << std::endl; for (int i = 0; i < nBasis; i++) { for (int j = 0; j < nBasis; j++) - std::cout << " " << std::setw(14) << M(i, j); + std::cout << " " << std::setprecision(4) << std::setw(14) << M(i, j); std::cout << std::endl; } diff --git a/tests/core/Test_plaquette_stats.cc b/tests/core/Test_plaquette_stats.cc new file mode 100644 index 000000000..10b33535c --- /dev/null +++ b/tests/core/Test_plaquette_stats.cc @@ -0,0 +1,211 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: tests/core/Test_plaquette_stats.cc + +Copyright (C) 2015 + +Author: Chulwoo Jung + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +/** + * Test_plaquette_stats + * + * Measure every plaquette Re Tr[U_mu(x) U_nu(x+mu) U_mu†(x+nu) U_nu†(x)] / Nc + * across all sites and all (mu,nu) planes and report max, min, and average. + * + * Usage: + * ./Test_plaquette_stats [Grid options] [--file ] [--hot] + * + * --file Read gauge field from NERSC-format file + * --hot Use random (hot) SU(3) start (default: cold/unit start) + * + * Grid size defaults to 8^3 x 16; override with --grid (e.g. --grid 4.4.4.8). + */ + +#include + +using namespace std; +using namespace Grid; + +// Return the plane label string for output +static std::string planeName(int mu, int nu) +{ + const char dirs[] = "xyzt"; + std::string s; + s += dirs[mu]; + s += dirs[nu]; + return s; +} + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + // ---- Lattice geometry ---- + Coordinate latt_size = GridDefaultLatt(); + if (latt_size.size() == 0) { + latt_size = {8, 8, 8, 16}; + } + Coordinate simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + GridCartesian grid(latt_size, simd_layout, mpi_layout); + + std::cout << GridLogMessage << "Lattice: "; + for (int d = 0; d < Nd; d++) + std::cout << latt_size[d] << (d < Nd-1 ? "x" : "\n"); + + // ---- Gauge field ---- + LatticeGaugeField Umu(&grid); + + // Check for --file argument + std::string config_file = ""; + for (int i = 1; i < argc - 1; i++) { + if (std::string(argv[i]) == "--file") { + config_file = argv[i+1]; + break; + } + } + + bool doHot = false; + for (int i = 1; i < argc; i++) { + if (std::string(argv[i]) == "--hot") { doHot = true; break; } + } + + if (!config_file.empty()) { + std::cout << GridLogMessage << "Reading gauge field from " << config_file << std::endl; + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, config_file); + } else { + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&grid); + pRNG.SeedFixedIntegers(seeds); + if (doHot) { + std::cout << GridLogMessage << "Generating hot (random SU(3)) start" << std::endl; + SU::HotConfiguration(pRNG, Umu); + } else { + std::cout << GridLogMessage << "Using cold (unit) gauge start" << std::endl; + SU::ColdConfiguration(pRNG, Umu); + } + } + + // ---- Extract link matrices ---- + std::vector U(Nd, &grid); + for (int mu = 0; mu < Nd; mu++) + U[mu] = PeekIndex(Umu, mu); + + // ---- Per-plane plaquette statistics ---- + // + // For each (mu, nu) plane (mu > nu) compute + // P_munu(x) = Re Tr[plaquette] / Nc + // then report max, min, mean over all sites. + // + // WilsonLoops::traceDirPlaquette gives Tr[U_mu U_nu U_mu† U_nu†] (complex). + + double vol = grid.gSites(); + + // Accumulate site-average plaquette (sum over planes / Nplanes / Nc) + LatticeComplex plaq_all(&grid); + plaq_all = Zero(); + int Nplanes = 0; + + std::cout << GridLogMessage + << "======== Per-plane plaquette statistics ========" << std::endl; + std::cout << GridLogMessage + << std::setw(6) << "plane" + << std::setw(20) << "max" + << std::setw(20) << "min" + << std::setw(20) << "average" + << std::endl; + + for (int mu = 1; mu < Nd; mu++) { + for (int nu = 0; nu < mu; nu++) { + + // Per-site trace of plaquette in (mu,nu) plane + LatticeComplex sitePlaq(&grid); + ColourWilsonLoops::traceDirPlaquette(sitePlaq, U, mu, nu); + + plaq_all = plaq_all + sitePlaq; + Nplanes++; + + // --- global average via sum() --- + TComplex Tsum = sum(sitePlaq); + ComplexD csum = TensorRemove(Tsum); + RealD avg = csum.real() / vol / Nc; + + // --- global max and min via unvectorize + GlobalMax/GlobalMin --- + std::vector sv; + unvectorizeToLexOrdArray(sv, sitePlaq); + + RealD local_max = -1e38, local_min = 1e38; + for (auto& tc : sv) { + RealD r = TensorRemove(tc).real() / Nc; + if (r > local_max) local_max = r; + if (r < local_min) local_min = r; + } + + // Reduce across MPI ranks + grid.GlobalMax(local_max); + grid.GlobalMin(local_min); + + std::cout << GridLogMessage + << std::setw(6) << planeName(mu, nu) + << std::setw(20) << std::setprecision(10) << local_max + << std::setw(20) << std::setprecision(10) << local_min + << std::setw(20) << std::setprecision(10) << avg + << std::endl; + } + } + + // ---- Overall (averaged over all planes) statistics ---- + // plaq_all = sum of Tr[...] over all 6 (mu,nu) planes + // Normalise to Re Tr / Nc per plane + plaq_all = plaq_all * (1.0 / Nc / Nplanes); + + TComplex Tsum_all = sum(plaq_all); + RealD avg_all = TensorRemove(Tsum_all).real() / vol; + + std::vector sv_all; + unvectorizeToLexOrdArray(sv_all, plaq_all); + + RealD max_all = -1e38, min_all = 1e38; + for (auto& tc : sv_all) { + RealD r = TensorRemove(tc).real(); + if (r > max_all) max_all = r; + if (r < min_all) min_all = r; + } + grid.GlobalMax(max_all); + grid.GlobalMin(min_all); + + // Cross-check with built-in avgPlaquette + RealD avg_builtin = ColourWilsonLoops::avgPlaquette(Umu); + + std::cout << GridLogMessage + << "======== Overall plaquette statistics (all planes) ========" << std::endl; + std::cout << GridLogMessage << " max = " << max_all << std::endl; + std::cout << GridLogMessage << " min = " << min_all << std::endl; + std::cout << GridLogMessage << " average = " << avg_all << std::endl; + std::cout << GridLogMessage << " avgPlaquette (builtin check) = " << avg_builtin << std::endl; + + Grid_finalize(); + return 0; +}