diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index b0817932d..90a720151 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -90,6 +90,7 @@ NAMESPACE_CHECK(multigrid); #include #include #include +#include #include #include #include diff --git a/Grid/algorithms/iterative/Gamma5BlockLanczos.h b/Grid/algorithms/iterative/Gamma5BlockLanczos.h index e4156c875..0ff243f6c 100644 --- a/Grid/algorithms/iterative/Gamma5BlockLanczos.h +++ b/Grid/algorithms/iterative/Gamma5BlockLanczos.h @@ -112,19 +112,26 @@ public: assert(nrm > 1e-14 && "first starting vector is zero"); u0 *= (1.0 / nrm); +#if 0 +//HACK + applyGamma5(u0,u1); +#else u1 = v1; ComplexD proj = innerProduct(u0, u1); u1 -= u0 * proj; nrm = std::sqrt(norm2(u1)); assert(nrm > 1e-14 && "second starting vector is linearly dependent on first"); u1 *= (1.0 / nrm); +#endif CMat2 G1 = gramMatrix(u0, u1); { Eigen::SelfAdjointEigenSolver es(G1); auto evals = es.eigenvalues(); std::cout << GridLogMessage - << "Gamma5BlockLanczos: G1 eigenvalues = " + << "Gamma5BlockLanczos: G1 = " <doVerify) verify("iter= "+std::to_string(iter)); + if(this->doVerify){ verify("iter= "+std::to_string(iter)); exit(-42);} int nRitz = (int)residuals_.size(); if (nRitz == 0) { @@ -1308,6 +1315,15 @@ 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); + { + CMat2 gamma1 = gramMatrix(r1, r2); + Eigen::ComplexEigenSolver esB(gamma1); + auto evB = esB.eigenvalues(); + std::cout << GridLogMessage + << "Gamma5BlockLanczos: gamma "< + + γ5-Block Lanczos example for the Wilson Dirac operator. + + Reads a gauge configuration from "config" (NERSC format) and Lanczos + parameters from "LanParams.xml". Runs Gamma5BlockLanczos to + compute eigenvalues of D_W directly (not H_W = γ5 D_W). + + 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 */ + +#include + +#include +#include +#include + +using namespace std; +using namespace Grid; + +namespace Grid { + +struct LanczosParameters : Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters, + RealD, mass, + Integer, Nstop, + Integer, Nk, + Integer, Np, + Integer, maxIter, + Integer, reorthog, + Integer, verify, + Integer, ReadEvec, + RealD, resid) + + LanczosParameters() + : mass(-0.5), Nstop(10), Nk(20), maxIter(100), + reorthog(1), verify(0), ReadEvec(0), resid(1e-8) + {} + + template + LanczosParameters(Reader& r) { initialize(r); } + + template + void initialize(Reader& r) { + read(r, "LanczosParameters", *this); + } +}; + +} // namespace Grid + +template +void writeField(T& in, std::string const& fname) { +#ifdef HAVE_LIME + std::cout << GridLogMessage << "Writing to: " << fname << std::endl; + Grid::emptyUserRecord record; + Grid::ScidacWriter WR(in.Grid()->IsBoss()); + WR.open(fname); + WR.writeScidacFieldRecord(in, record, 0); + WR.close(); +#endif +} + +typedef WilsonFermionD WilsonOp; +typedef typename WilsonFermionD::FermionField FermionField; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), + GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + std::vector seeds4({1, 2, 3, 4}); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Read gauge configuration + LatticeGaugeField Umu(UGrid); + FieldMetaData header; + std::string configFile("config"); + NerscIO::readConfiguration(Umu, header, configFile); + std::cout << GridLogMessage << "Loaded gauge configuration: " << configFile << std::endl; + + // Read Lanczos parameters + LanczosParameters Params; + { + XmlReader rd("LanParams.xml"); + read(rd, "LanczosParameters", Params); + } + std::cout << GridLogMessage << Params << std::endl; + + // Build Wilson Dirac operator and wrap in a non-Hermitian linear operator + std::vector boundary = {1, 1, 1, -1}; + WilsonOp::ImplParams WilsonParams(boundary); + WilsonOp Dwilson(Umu, *UGrid, *UrbGrid, Params.mass, WilsonParams); + NonHermitianLinearOperator DLinOp(Dwilson); + + // γ5 functor: for 4D Wilson fermions γ5 is Gamma(Gamma5) + Gamma G5(Gamma::Algebra::Gamma5); + auto gamma5 = [&G5](const FermionField& in, FermionField& out) { + out = G5 * in; + }; + + // Starting vectors: two independent random vectors + FermionField src(UGrid), src2(UGrid); + random(RNG4, src); + random(RNG4, src2); + std::cout << GridLogMessage << "Using two random starting vectors" << std::endl; + + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "*******************************************" << std::endl; + std::cout << GridLogMessage << " Running γ5-Block Lanczos" << std::endl; + std::cout << GridLogMessage << " mass = " << Params.mass << std::endl; + std::cout << GridLogMessage << " Nk = " << Params.Nk << std::endl; + std::cout << GridLogMessage << " maxIter = " << Params.maxIter << std::endl; + std::cout << GridLogMessage << " Nstop = " << Params.Nstop << std::endl; + std::cout << GridLogMessage << " resid = " << Params.resid << std::endl; + std::cout << GridLogMessage << " reorthog = " << Params.reorthog << std::endl; + std::cout << GridLogMessage << " verify = " << Params.verify << std::endl; + std::cout << GridLogMessage << "*******************************************" << std::endl; + std::cout << GridLogMessage << std::endl; + + Gamma5BlockLanczos G5BL(DLinOp, UGrid, gamma5, Params.resid); + G5BL.doEvalCheck = (Params.verify != 0); + G5BL.doVerify = (Params.verify != 0); +// G5BL(src, src2, Params.maxIter, Params.Nstop, Params.reorthog != 0); + G5BL.restart(src, src2, Params.maxIter, Params.Nk+Params.Np, Params.Nk, Params.Nstop, Params.reorthog != 0, EvalNormSmall); +// G5BL.implicitRestart(src, src2, Params.maxIter, Params.Nk+Params.Np, Params.Nk, Params.Nstop, Params.reorthog != 0, EvalNormSmall); + if (Params.verify ) G5BL.verify("after restart"); + + // Summary of results + Eigen::VectorXcd evals = G5BL.getEvals(); + std::vector residuals = G5BL.getResiduals(); + std::vector evecs = G5BL.getEvecs(); + + int Nout = (int)evals.size(); + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "*******************************************" << std::endl; + std::cout << GridLogMessage << " γ5-Block Lanczos: " << Nout << " Ritz pairs" << std::endl; + std::cout << GridLogMessage << "*******************************************" << std::endl; + for (int i = 0; i < Nout; i++) { + std::cout << GridLogMessage + << " [" << std::setw(3) << i << "]" + << " lambda = " << evals(i) + << " |res| = " << residuals[i] << std::endl; + } + + // Write the first Nstop eigenvectors (in SCIDAC format when LIME is available) + int Nwrite = std::min((int)Params.Nstop, Nout); + for (int i = 0; i < Nwrite; i++) { + std::string fname = "./g5bl_evec_m" + std::to_string(Params.mass) + + "_" + std::to_string(i); + writeField(evecs[i], fname); + } + + std::cout << GridLogMessage << std::endl; + std::cout << GridLogMessage << "Done" << std::endl; + + Grid_finalize(); + return 0; +} diff --git a/examples/LanParams.xml b/examples/LanParams.xml index 0697cae5a..af518497d 100644 --- a/examples/LanParams.xml +++ b/examples/LanParams.xml @@ -5,15 +5,16 @@ -0.025 1.8 48 - 10 - 10 - 10 + 800 + 800 + 100 0 1000 1 4 1 - 1e-10 + 1.5 + 1e-8 1 100 51 diff --git a/examples/fft5d.cc b/examples/fft5d.cc new file mode 100644 index 000000000..abfa8f6fe --- /dev/null +++ b/examples/fft5d.cc @@ -0,0 +1,465 @@ +/************************************************************************************* + + fft5d.cc — Fourier analysis of a time series of 4-D lattice scalar fields. + + Designed for force-norm files from FTHMC (one RealD per site, SCIDAC or binary). + Assembles a (4+1)-D structure [trajectory][t][z][y][x] using Grid lattice fields + and performs three FFT analyses: + + --fft spatial + 4-D spatial FFT using Grid's FFT class (MPI+GPU parallel via Cshift/cufft). + Outputs shell-averaged P(|k|^2) averaged over trajectories, and a + per-mode table P(kt,kz,ky,kx). + + --fft traj + 1-D trajectory-axis FFT at each lattice site using FFTW3 locally at each + MPI rank. Outputs site-averaged P(f_traj) and autocorrelation C(lag). + Results are combined with GlobalSumVector. + + --fft all + Both of the above, plus a 2-D cross spectrum P(f_traj, |k_spatial|^2): + spatial FFT per trajectory, then trajectory FFTW on the momentum-space series. + + Normalisations: + Spatial FFT (volume V): P_spatial(k) = |F(k)|^2 / V^2 + => (1/V) * sum_k P(k) = site mean-square per trajectory + Trajectory FFT (Ntraj): P_traj(f) = |F(f)|^2 / Ntraj^2 + => sum_f P(f) = site mean-square over trajectories + + Build: add to examples/Make.inc (see bottom of this file), then make. + + Usage (follows Grid conventions): + fft5d --grid Lx.Ly.Lz.Lt [--mpi Px.Py.Pz.Pt] \ + [--fft spatial|traj|all] [--format scidac|binary] \ + [--output PREFIX] file1 file2 ... + +*************************************************************************************/ + +#include +// Grid's FFT.h uses cufft on CUDA builds; for the trajectory axis we also need +// CPU-side FFTW3 (already linked via -lfftw3 in GRID_LIBS). +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace Grid; + +// ───────────────────────────────────────────────────────────────────────────── +// Global coordinate of a local site (osite, lane) on this MPI rank +// ───────────────────────────────────────────────────────────────────────────── +static void globalCoor(int osite, int lane, GridCartesian* g, Coordinate& gc) +{ + Coordinate oc, ic; + Lexicographic::CoorFromIndex(oc, osite, g->_rdimensions); + Lexicographic::CoorFromIndex(ic, lane, g->_simd_layout); + for (int d = 0; d < Nd; d++) + gc[d] = g->_processor_coor[d] * g->_ldimensions[d] + + oc[d] * g->_simd_layout[d] + ic[d]; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Shell-average a power LatticeRealD over |k|^2 bins (MPI-aware via GlobalSumVector) +// ───────────────────────────────────────────────────────────────────────────── +static void shellAverage(const LatticeRealD& power, GridCartesian* g, + double norm, std::ostream& ofs) +{ + const int Nsimd = vRealD::Nsimd(); + Coordinate fdims = g->_fdimensions; + int maxk2 = 0; + for (int d = 0; d < Nd; d++) { int h = fdims[d]/2; maxk2 += h*h; } + + std::vector psum(maxk2+1, 0.0), cnt(maxk2+1, 0.0); + { + auto pv = power.View(CpuRead); + Coordinate gc(Nd); + for (int os = 0; os < (int)g->oSites(); os++) { + for (int lane = 0; lane < Nsimd; lane++) { + globalCoor(os, lane, g, gc); + int k2 = 0; + for (int d = 0; d < Nd; d++) { + int kd = std::min(gc[d], fdims[d] - gc[d]); + k2 += kd*kd; + } + psum[k2] += (RealD)extractLane(lane, pv[os]); + cnt [k2] += 1.0; + } + } + } + g->GlobalSumVector(psum.data(), (int)psum.size()); + g->GlobalSumVector(cnt .data(), (int)cnt .size()); + + for (int k2 = 0; k2 <= maxk2; k2++) + if (cnt[k2] > 0) + ofs << k2 << " " << std::sqrt((double)k2) << " " + << (psum[k2] / cnt[k2]) * norm << "\n"; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Convert LatticeRealD → LatticeComplexD (zero imaginary part) +// ───────────────────────────────────────────────────────────────────────────── +static LatticeComplexD toComplex(const LatticeRealD& r) +{ + std::vector lr; + unvectorizeToLexOrdArray(lr, r); + std::vector lc(lr.size()); + for (size_t i = 0; i < lr.size(); i++) lc[i] = ComplexD(lr[i], 0.0); + LatticeComplexD c(r.Grid()); + vectorizeFromLexOrdArray(lc, c); + return c; +} + +// ───────────────────────────────────────────────────────────────────────────── +// File readers +// ───────────────────────────────────────────────────────────────────────────── +static LatticeComplexD readScidac(const std::string& fname, GridCartesian* g) +{ + LatticeRealD field(g); + emptyUserRecord rec; + ScidacReader RD; + RD.open(fname); + RD.readScidacFieldRecord(field, rec); + RD.close(); + return toComplex(field); +} + +static LatticeComplexD readBinary(const std::string& fname, GridCartesian* g) +{ + // Raw IEEE doubles in lex order [x][y][z][t] written serially. + // Rank 0 reads the file, broadcasts to all ranks via GlobalSumVector. + Coordinate fdims = g->_fdimensions; + long vol4 = 1; for (int d = 0; d < Nd; d++) vol4 *= fdims[d]; + + std::vector buf(vol4, 0.0); + if (g->IsBoss()) { + std::ifstream f(fname, std::ios::binary); + if (!f) throw std::runtime_error("Cannot open: " + fname); + f.read(reinterpret_cast(buf.data()), vol4 * sizeof(RealD)); + } + g->GlobalSumVector(buf.data(), (int)vol4); // broadcast from rank 0 + + LatticeRealD field(g); + Coordinate ldims = g->_ldimensions, pcoor = g->_processor_coor; + for (long ls = 0; ls < g->lSites(); ls++) { + Coordinate lc; + g->LocalIndexToLocalCoor(ls, lc); + long glex = 0, stride = 1; + for (int d = 0; d < Nd; d++) { + glex += (pcoor[d]*ldims[d] + lc[d]) * stride; + stride *= fdims[d]; + } + pokeLocalSite(buf[glex], field, lc); + } + return toComplex(field); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Load all files into a trajectory vector +// ───────────────────────────────────────────────────────────────────────────── +static std::vector +loadFiles(const std::vector& files, GridCartesian* g, + const std::string& fmt) +{ + std::vector traj; + traj.reserve(files.size()); + for (int n = 0; n < (int)files.size(); n++) { + std::cout << GridLogMessage << "[" << n << "] reading " << files[n] << "\n"; + if (fmt == "scidac") traj.push_back(readScidac(files[n], g)); + else traj.push_back(readBinary (files[n], g)); + } + return traj; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Spatial FFT analysis +// ───────────────────────────────────────────────────────────────────────────── +static void analyzeSpatial(const std::vector& traj, + GridCartesian* g, const std::string& pfx) +{ + int Ntraj = (int)traj.size(); + long vol4 = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d]; + + FFT theFFT(g); + LatticeRealD pavg(g); pavg = Zero(); + for (int n = 0; n < Ntraj; n++) { + LatticeComplexD fk(g); + theFFT.FFT_all_dim(fk, traj[n], FFT::forward); + // Evaluate product before applying real() — real() is not defined for + // unevaluated LatticeBinaryExpression. + LatticeComplexD fk_sq(g); fk_sq = conjugate(fk) * fk; + std::vector lc; unvectorizeToLexOrdArray(lc, fk_sq); + std::vector lr(lc.size()); + for (size_t i = 0; i < lc.size(); i++) lr[i] = lc[i].real(); + LatticeRealD pk(g); vectorizeFromLexOrdArray(lr, pk); + pavg += pk; + } + pavg *= (1.0 / Ntraj); + + // Shell-averaged spectrum + { + std::ofstream fs(pfx + "_spatial_shell.dat"); + fs << "# k2 |k| P_shell_avg (P = |F|^2 / V^2 / shell_count)\n"; + shellAverage(pavg, g, 1.0 / ((double)vol4 * vol4), fs); + std::cout << GridLogMessage << "Written: " << pfx << "_spatial_shell.dat\n"; + } + + // Per-mode table (rank 0 only, via peekSite) + if (g->IsBoss()) { + std::ofstream fm(pfx + "_spatial_modes.dat"); + fm << "# kt kz ky kx k2 |k| P_traj_avg\n"; + Coordinate fdims = g->_fdimensions, site(Nd); + double norm = 1.0 / ((double)vol4 * vol4); + for (site[3]=0; site[3]& in, + std::vector& out) +{ + int Ntraj = (int)in.size(); + 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); + 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(); + } + } + + // lsites transforms of length Ntraj, stride=lsites, dist=1 + std::vector obuf((long)Ntraj * lsites); + 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, + 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. + 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]); + vectorizeFromLexOrdArray(lc, out[k]); + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Trajectory-axis FFT analysis +// ───────────────────────────────────────────────────────────────────────────── +static void analyzeTraj(const std::vector& traj, + GridCartesian* g, const std::string& pfx) +{ + int Ntraj = (int)traj.size(); + int Nf = Ntraj / 2 + 1; + long vol4 = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d]; + + std::vector ftraj; + trajFFT(traj, ftraj); + + // P(k) = (1/vol4) * sum_sites |F_traj(k)|^2 / Ntraj^2 + std::vector Pavg(Ntraj, 0.0); + for (int k = 0; k < Ntraj; k++) { + LatticeComplexD tmp(g); tmp = conjugate(ftraj[k]) * ftraj[k]; + std::vector lc; unvectorizeToLexOrdArray(lc, tmp); + std::vector lr(lc.size()); + for (size_t i = 0; i < lc.size(); i++) lr[i] = lc[i].real(); + LatticeRealD pk(g); vectorizeFromLexOrdArray(lr, pk); + RealD s = 0.0; + for (long ls = 0; ls < g->lSites(); ls++) { + Coordinate lc; g->LocalIndexToLocalCoor(ls, lc); + RealD v; peekLocalSite(v, pk, lc); + s += v; + } + g->GlobalSum(s); + Pavg[k] = s / ((double)vol4 * Ntraj * Ntraj); + } + + { + std::ofstream f(pfx + "_traj_power.dat"); + f << "# k freq P_avg (sum_k P = site mean-square)\n"; + for (int k = 0; k < Nf; k++) + f << k << " " << (double)k/Ntraj << " " << Pavg[k] << "\n"; + std::cout << GridLogMessage << "Written: " << pfx << "_traj_power.dat\n"; + } + + // 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 acorr(Ntraj, 0.0); + fftw_plan ip = fftw_plan_dft_c2r(1, &Ntraj, Pc.data(), acorr.data(), FFTW_ESTIMATE); + fftw_execute(ip); fftw_destroy_plan(ip); + double c0 = acorr[0] / Ntraj; + { + std::ofstream f(pfx + "_traj_autocorr.dat"); + f << "# lag C(lag) C(lag)/C(0)\n"; + for (int lag = 0; lag < Ntraj/2; lag++) { + double c = acorr[lag] / Ntraj; + f << lag << " " << c << " " << (c0 > 0 ? c/c0 : 0.0) << "\n"; + } + std::cout << GridLogMessage << "Written: " << pfx << "_traj_autocorr.dat\n"; + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Full 5-D analysis: spatial FFT → trajectory FFT → 2-D cross spectrum +// ───────────────────────────────────────────────────────────────────────────── +static void analyzeAll5D(const std::vector& traj, + GridCartesian* g, const std::string& pfx) +{ + int Ntraj = (int)traj.size(); + int Nf = Ntraj / 2 + 1; + long vol4 = 1; for (int d = 0; d < Nd; d++) vol4 *= g->_fdimensions[d]; + Coordinate fdims = g->_fdimensions; + + // Step 1: spatial FFT for each trajectory + FFT theFFT(g); + std::vector sfft(Ntraj, LatticeComplexD(g)); + for (int n = 0; n < Ntraj; n++) + theFFT.FFT_all_dim(sfft[n], traj[n], FFT::forward); + + // Step 2: trajectory FFTW on the momentum-space series + std::vector tfft; + trajFFT(sfft, tfft); + + // Step 3: P(f_traj, |k_spatial|^2) shell-averaged + int maxk2 = 0; + for (int d = 0; d < Nd; d++) { int h = fdims[d]/2; maxk2 += h*h; } + + long nb = (long)Nf * (maxk2+1); + std::vector p2d(nb, 0.0), cnt2d(nb, 0.0); + const int Nsimd = vComplexD::Nsimd(); + double norm = 1.0 / ((double)Ntraj * Ntraj * vol4 * vol4); + + for (int kf = 0; kf < Nf; kf++) { + LatticeComplexD tmp(g); tmp = conjugate(tfft[kf]) * tfft[kf]; + std::vector lc2; unvectorizeToLexOrdArray(lc2, tmp); + std::vector lr2(lc2.size()); + for (size_t i = 0; i < lc2.size(); i++) lr2[i] = lc2[i].real(); + LatticeRealD pk(g); vectorizeFromLexOrdArray(lr2, pk); + auto pv = pk.View(CpuRead); + Coordinate gc(Nd); + for (int os = 0; os < (int)g->oSites(); os++) { + for (int lane = 0; lane < Nsimd; lane++) { + globalCoor(os, lane, g, gc); + int k2 = 0; + for (int d = 0; d < Nd; d++) { + int kd = std::min(gc[d], fdims[d] - gc[d]); + k2 += kd*kd; + } + long b = (long)kf*(maxk2+1) + k2; + p2d [b] += (RealD)extractLane(lane, pv[os]); + cnt2d[b] += 1.0; + } + } + } + g->GlobalSumVector(p2d .data(), (int)nb); + g->GlobalSumVector(cnt2d.data(), (int)nb); + + if (g->IsBoss()) { + std::ofstream f(pfx + "_5d_cross.dat"); + f << "# k_traj freq_traj k2_spatial |k_spatial| P_avg\n"; + for (int kf = 0; kf < Nf; kf++) { + double freq = (double)kf / Ntraj; + for (int k2 = 0; k2 <= maxk2; k2++) { + long b = (long)kf*(maxk2+1) + k2; + if (cnt2d[b] > 0) + f << kf << " " << freq << " " << k2 << " " + << std::sqrt((double)k2) << " " + << (p2d[b]/cnt2d[b]) * norm << "\n"; + } + f << "\n"; // blank line between kf slices for gnuplot pm3d + } + std::cout << GridLogMessage << "Written: " << pfx << "_5d_cross.dat\n"; + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Main +// ───────────────────────────────────────────────────────────────────────────── +static void usage(const char* prog) +{ + std::cerr + << "Usage: " << prog << " --grid Lx.Ly.Lz.Lt [--mpi Px.Py.Pz.Pt]\n" + << " [--fft spatial|traj|all] [--format scidac|binary]\n" + << " [--output PREFIX] file1 [file2 ...]\n"; + exit(1); +} + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); // consumes --grid, --mpi, etc. + + std::string fftMode = "all"; + std::string fmt = "scidac"; + std::string pfx = "fft5d"; + std::vector files; + + for (int i = 1; i < argc; i++) { + std::string a = argv[i]; + if (a == "--fft" && i+1 < argc) fftMode = argv[++i]; + else if (a == "--format" && i+1 < argc) fmt = argv[++i]; + else if (a == "--output" && i+1 < argc) pfx = argv[++i]; + else if (!a.empty() && a[0] != '-') files.push_back(a); + // unknown --flags skipped (may be Grid flags already consumed) + } + + if (files.empty()) usage(argv[0]); + + GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), + GridDefaultMpi()); + + Coordinate latt = GridDefaultLatt(); + std::cout << GridLogMessage << "Lattice : " + << latt[0]<<"x"<