diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h index aef53cc72..c7f61b217 100644 --- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h +++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h @@ -85,6 +85,9 @@ public: // the iteration. GCRCoefficients *Recorder = nullptr; void SetCoefficientRecorder(GCRCoefficients *r) { Recorder = r; if(r) r->mmax = mmax; }; + // Free the persistent history (e.g. when this solver is replaced by a + // replayed polynomial): re-made on the next call if ever needed again. + void ReleaseHistory(void) { q.clear(); p.clear(); qq.clear(); hist_grid = nullptr; }; PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase &_Linop,LinearFunction &Prec,int _mmax,int _nstep) : Tolerance(tol), diff --git a/Grid/algorithms/multigrid/BlockCyclic.h b/Grid/algorithms/multigrid/BlockCyclic.h index 3012bd904..c8c7879be 100644 --- a/Grid/algorithms/multigrid/BlockCyclic.h +++ b/Grid/algorithms/multigrid/BlockCyclic.h @@ -62,13 +62,19 @@ public: // For fixed P the per-rank SUMMA volume N^2 (1/Pr + 1/Pc) is minimised // at the most square grid. P=288 -> 16 x 18. /////////////////////////////////////////////////////////////////////////// + // Nearest-to-square factorisation with Pr >= Pc. The orientation matters + // for SUMMA: B panels travel along process columns in slots of + // S = ceil(Pc/Pr) panels, padded; with Pr >= Pc, S == 1 and there is no + // padding. (288 as 16x18 had S=2 -- half the B-ring bytes were zeros; + // 18x16 moves ~30% fewer bytes for the same inverse.) static void ChooseProcessGrid(int P, int &Pr, int &Pc) { GRID_ASSERT(P >= 1); - Pr = 1; - for(int r=1; (int64_t)r*r <= (int64_t)P; r++) - if ( P % r == 0 ) Pr = r; - Pc = P / Pr; + int r = 1; + for(int f=1; (int64_t)f*f <= (int64_t)P; f++) + if ( P % f == 0 ) r = f; + Pc = r; + Pr = P / r; } /////////////////////////////////////////////////////////////////////////// diff --git a/Grid/algorithms/multigrid/BlockCyclicSumma.h b/Grid/algorithms/multigrid/BlockCyclicSumma.h index 625667e58..33f3cefa0 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSumma.h +++ b/Grid/algorithms/multigrid/BlockCyclicSumma.h @@ -132,6 +132,15 @@ public: // row, ring B along the process column, then the local GEMMs. The rings // are synchronous SendToRecvFrom, so tRing is time the GPU is idle unless // a future version overlaps them with the GEMMs. + // PERSISTENT ring buffers: allocated once (grow-only) and reused by every + // Multiply, so the device addresses handed to MPI never change. Fresh + // per-call buffers rotated through the caching allocator's blocks, and a + // GPU-direct RDMA registration cache keyed on address then re-registers + // per message: measured 1.26 GB/s/rank on 13 MB ring messages (production, + // GRID_ALLOC_NCACHE_LARGE=64) against 62 s total for the same inverse when + // hipMalloc returned a stable address. + deviceVector Abuf; + deviceVector Bbuf; double tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0; uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0; void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; } @@ -199,8 +208,8 @@ public: const uint64_t slotB = (uint64_t)S*slotB1; nMultiply++; tAlloc -= usecond(); - deviceVector Abuf( slotA*Pc ? slotA*Pc : 1 ); - deviceVector Bbuf( slotB*Pr ? slotB*Pr : 1 ); + if ( Abuf.size() < std::max(slotA*Pc,1) ) Abuf.resize( std::max(slotA*Pc,1) ); + if ( Bbuf.size() < std::max(slotB*Pr,1) ) Bbuf.resize( std::max(slotB*Pr,1) ); tAlloc += usecond(); deviceVector ap(1), bp(1), cp(1); diff --git a/Grid/algorithms/multigrid/Smoothers.h b/Grid/algorithms/multigrid/Smoothers.h index 08fef288b..938a188e2 100644 --- a/Grid/algorithms/multigrid/Smoothers.h +++ b/Grid/algorithms/multigrid/Smoothers.h @@ -91,6 +91,8 @@ public: LinearOperatorBase &Linop; RealD lo, hi; int order; std::vector Coeffs; + int Verbose = 0; + std::string name = "cheb"; ChebyshevNonHermitianSmoother(RealD _lo,RealD _hi,int _order,LinearOperatorBase &Op) : Linop(Op), lo(_lo), hi(_hi), order(_order) { @@ -126,6 +128,7 @@ public: if ( Coeffs[n] != 0.0 ) axpy(out,Coeffs[n],*Tnp,out); Field *swizzle=Tnm; Tnm=Tn; Tn=Tnp; Tnp=swizzle; } + if ( Verbose ) { Linop.Op(out,y); y = y - in; std::cout << GridLogMessage << " " << name << " cheb |r|/|r0| = " << std::sqrt(norm2(y)/norm2(in)) << std::endl; } } }; @@ -205,16 +208,27 @@ public: std::vector > b; GridBase *hist_grid = nullptr; std::vector p; + int Verbose = 0; // 1: print |r_m|/|r_0| per call (one extra reduction) + std::string name = "replay"; GCRReplaySmoother(LinearOperatorBase &Op, const GCRCoefficients &c) : Linop(Op) { mmax = c.mmax; GRID_ASSERT(mmax>=1); nstep = c.Steps(); GRID_ASSERT(nstep>=1); a.resize(nstep); b.resize(nstep); + RealD amax=0.0, bmax=0.0; for(int k=0;k(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM)); + FineCheb->Verbose = PolyVerbose; FineCheb->name = "Fsmoother"; FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev"); } if ( CoarseSmootherMode == "cheb" ) { CoarseCheb.reset(new ChebyshevNonHermitianSmoother(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC)); + CoarseCheb->Verbose = PolyVerbose; CoarseCheb->name = "Csmoother"; CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev"); } if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF); @@ -1005,13 +1009,17 @@ int main (int argc, char ** argv) SmootherGCR.SetCoefficientRecorder(nullptr); recF.Report("Fsmoother"); FineReplay.reset(new GCRReplaySmoother(ShiftedPVdagM,recF)); + FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother"; FineSmootherSlot.Set(*FineReplay,"Fsmoother replay"); + SmootherGCR.ReleaseHistory(); // memory-neutral swap: the GCR's history goes as the replay's comes } if ( CoarseSmootherMode == "replay" ) { CoarseSmootherGCR.SetCoefficientRecorder(nullptr); recC.Report("Csmoother"); CoarseReplay.reset(new GCRReplaySmoother(ShiftedC,recC)); + CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother"; CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay"); + CoarseSmootherGCR.ReleaseHistory(); } }; std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << " coarse " << CoarseSmootherMode diff --git a/systems/Frontier/schur2d_ladder.job b/systems/Frontier/schur2d_ladder.job index fa28bf100..71016a364 100644 --- a/systems/Frontier/schur2d_ladder.job +++ b/systems/Frontier/schur2d_ladder.job @@ -80,7 +80,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $root/tests/debug/Test_schur2d_scale \ ############################################################################## echo "=========================================================" -echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 16x18)" +echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 18x16)" echo " THE number: invert phase vs the 1D baseline 328 s" echo "=========================================================" ############################################################################## diff --git a/systems/Frontier/slate_vs_schur2d.job b/systems/Frontier/slate_vs_schur2d.job index d04a2d330..1aa92fa55 100644 --- a/systems/Frontier/slate_vs_schur2d.job +++ b/systems/Frontier/slate_vs_schur2d.job @@ -20,7 +20,7 @@ # # S1 : 1 node, 8 ranks, N=4096 -- API/device/thread-level shakeout # S2 : 1 node, 8 ranks, N=13824 -- one-node production-shape rehearsal -# S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18) -- THE comparison +# S3 : 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16) -- THE comparison # # S3 is gated on S1 passing (rc=0 and both certificates printed): a broken # SLATE leg should not burn the 36-node allocation. @@ -100,7 +100,7 @@ S2D_N=13824 srun -N1 -n8 ./select_gpu $BIN --mpi 1.1.2.4 --grid 16.16.16.16 $OPT ############################################################################## echo "=========================================================" -echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18)" +echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 18x16)" echo "=========================================================" ############################################################################## if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -ge 2 ] diff --git a/systems/Frontier/smoother_modes.job b/systems/Frontier/smoother_modes.job index 664219233..befb1dc24 100644 --- a/systems/Frontier/smoother_modes.job +++ b/systems/Frontier/smoother_modes.job @@ -23,6 +23,7 @@ # heavy one at Nrhs=1) # M4 cheb / gcr fine Chebyshev [FineChebLo,FineChebHi] order Fso # M5 cheb / cheb +# M6 gcr / replay coarse frozen only # # Laptop 8^4 findings (hot config, Ls=4, NBASIS=8): replay/replay converges # (28 vs 23 outer); cheb on the FINE level diverges there while cheb on the @@ -96,6 +97,7 @@ export SmootherCoeffLog=0 # frozen-polynomial controls export PolyRecordIters=4 # outer steps recorded before the switch +export PolyVerbose=1 # frozen smoothers print |r_m|/|r_0| per call: separates 'bad polynomial' from 'linear V-cycle stagnates the outer' export FineChebLo=3.0 # harvested |R|<0.1 edge / PowerIteration edge x1.05 export FineChebHi=137.0 export CoarseChebLo=8.0 @@ -105,8 +107,9 @@ run_mode () { name=$1; export FineSmootherMode=$2; export CoarseSmootherMode=$3 echo "----- $name : FineSmootherMode=$FineSmootherMode CoarseSmootherMode=$CoarseSmootherMode -----" fname=log.modes.$name - srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ + srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1 + echo " exit $?"; sleep 60 # let a faulted step's GPUs be released before the next srun (M3 after M2 faulted instantly) echo " $(grep -h 'V2 3-level solve Nrhs' $fname | tr '\n' ' ')" echo " $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | tr '\n' ' ')" echo " $(grep -h 'FINAL Nrhs .: worst' $fname | tr '\n' ' ')" @@ -118,6 +121,7 @@ run_mode M2_replay_replay replay replay run_mode M3_replay_gcr replay gcr run_mode M4_cheb_gcr cheb gcr run_mode M5_cheb_cheb cheb cheb +run_mode M6_gcr_replay gcr replay # coarse frozen only: M2 showed Couter 16 -> 5 steps with it live echo "=========================================================" echo "summary" diff --git a/tests/debug/Test_blockcyclic.cc b/tests/debug/Test_blockcyclic.cc index f9e6388ad..e5b0d04ff 100644 --- a/tests/debug/Test_blockcyclic.cc +++ b/tests/debug/Test_blockcyclic.cc @@ -38,7 +38,7 @@ Author: Peter Boyle // equal mloc*nloc; LocalOffset is a bijection onto [0,mloc*nloc). // T5 : block contiguity: within any owned global block, consecutive // global rows are consecutive local rows (what SUMMA panels rely on). -// T6 : ChooseProcessGrid: exact factorisation, Pr<=Pc, most-square. +// T6 : ChooseProcessGrid: exact factorisation, Pr>=Pc, most-square. ////////////////////////////////////////////////////////////////////////////// #include @@ -227,15 +227,15 @@ int main(int argc, char **argv) int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(P,Pr,Pc); if ( Pr*Pc != P ) ok = false; - if ( Pr > Pc ) ok = false; - // most-square: no divisor r with Pr < r <= sqrt(P) - for(int r=Pr+1; (int64_t)r*r <= (int64_t)P; r++) + if ( Pr < Pc ) ok = false; // Pr >= Pc (SUMMA B panels unpadded) + // most-square: no divisor r with Pc < r <= sqrt(P) + for(int r=Pc+1; (int64_t)r*r <= (int64_t)P; r++) if ( P % r == 0 ) ok = false; } int Pr,Pc; BlockCyclicLayout::ChooseProcessGrid(288,Pr,Pc); - if ( !(Pr==16 && Pc==18) ) ok = false; - Report("T6 ChooseProcessGrid exact, Pr<=Pc, most-square (288 -> 16x18)", ok); + if ( !(Pr==18 && Pc==16) ) ok = false; // Pr >= Pc: SUMMA B panels unpadded (S = ceil(Pc/Pr) = 1) + Report("T6 ChooseProcessGrid exact, Pr>=Pc, most-square (288 -> 18x16)", ok); } //////////////////////////////////////////////////////////////////////// diff --git a/tests/debug/Test_poly_smoother.cc b/tests/debug/Test_poly_smoother.cc new file mode 100644 index 000000000..b3b7b2dae --- /dev/null +++ b/tests/debug/Test_poly_smoother.cc @@ -0,0 +1,200 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/debug/Test_poly_smoother.cc + + Copyright (C) 2026 + +Author: Peter Boyle + + 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. + + See the full license in the file "LICENSE" in the top level distribution + directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////// +// Fixed-polynomial smoothers (Smoothers.h) against the adaptive GCR they +// replace, on a shifted Wilson operator (non-Hermitian, spectrum in the right +// half plane): +// +// record : GCR(mmax=2, 8 steps) with a coefficient recorder on 16 sources +// T1 : GCRReplaySmoother on a FRESH source reduces the residual to +// within a factor 2 of the live GCR on the same source +// T2 : replay and GCR solutions agree to the coefficient spread (<10%) +// T3 : ChebyshevNonHermitianSmoother with Op=HermOp reproduces the legacy +// HermOp ChebyshevSmoother (real spectrum); on the complex-spectrum +// Wilson op its (poor) reduction is printed for information only +// T4 : the historical HermOp ChebyshevSmoother compiles with both +// constructor signatures and reduces the MdagM residual +// T5 : replay is bitwise repeatable (no reductions => no reordering) +// +// mpirun -n 1 ./Test_poly_smoother --grid 8.8.8.8 --mpi 1.1.1.1 +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +using namespace Grid; + +static int failures = 0; +static void Report(const std::string &name, bool pass, const std::string &detail="") +{ + std::cout << GridLogMessage << " " << name << (pass ? " PASS" : " ** FAIL **"); + if ( detail.size() ) std::cout << " " << detail; + std::cout << std::endl; + if ( !pass ) failures++; +} + +template +class ShiftedOp : public LinearOperatorBase { + LinearOperatorBase &_Op; RealD shift; +public: + ShiftedOp(RealD s, LinearOperatorBase &Op) : _Op(Op), shift(s) {} + void OpDiag (const Field &in, Field &out) { GRID_ASSERT(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { GRID_ASSERT(0); } + void OpDirAll(const Field &in, std::vector &out) { GRID_ASSERT(0); } + void Op (const Field &in, Field &out) { _Op.Op(in,out); out = out + shift*in; } + void AdjOp (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } + void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } +}; + +template +RealD Residual(LinearOperatorBase &Op, const Field &src, const Field &x) +{ + Field r(src.Grid()); Op.Op(x,r); r = r - src; + return std::sqrt(norm2(r)/norm2(src)); +} + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), + GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + + LatticeGaugeFieldD Umu(UGrid); + SU::HotConfiguration(RNG4, Umu); + + RealD mass = 0.5, shift = 0.1; + WilsonFermionD Dw(Umu, *UGrid, *UrbGrid, mass); + NonHermitianLinearOperator Op(Dw); + ShiftedOp SOp(shift, Op); + TrivialPrecon simple; + + const int mmax = 2, nstep = 8, ncal = 16; + + ////////////////////////////////////////////////////////////////////// + // record + ////////////////////////////////////////////////////////////////////// + PrecGeneralisedConjugateResidualNonHermitian GCR(0.0, 1, SOp, simple, mmax, nstep); + GCR.SetZeroGuess(1); GCR.Name("smoother"); + GCRCoefficients rec; + GCR.SetCoefficientRecorder(&rec); + LatticeFermionD src(UGrid), x(UGrid); + for(int c=0;c Replay(SOp, rec); + Replay(src,xr); + RealD rr = Residual(SOp,src,xr); + Report("T1 replay residual within 2x of live GCR", rr < 2.0*rg, + "GCR |r|/|r0| = "+std::to_string(rg)+" replay "+std::to_string(rr)); + + d = xr - xg; + RealD rel = std::sqrt(norm2(d)/norm2(xg)); + Report("T2 replay solution vs GCR solution", rel < 0.1, "rel "+std::to_string(rel)); + + ////////////////////////////////////////////////////////////////////// + // Chebyshev 1/x on [lo,hi], hi from a power iteration on SOp + ////////////////////////////////////////////////////////////////////// + RealD hi; + { + LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v); + RealD n = std::sqrt(norm2(v)); v = v*(1.0/n); + for(int i=0;i<60;i++){ SOp.Op(v,Av); hi = std::sqrt(norm2(Av)); v = Av*(1.0/hi); } + } + RealD lo = 0.5; + std::cout << GridLogMessage << "power iteration |lambda_max| ~ " << hi << " Chebyshev range [" << lo << "," << 1.05*hi << "]" << std::endl; + ChebyshevNonHermitianSmoother Cheb(lo, 1.05*hi, nstep, SOp); + Cheb(src,xc); + RealD rc = Residual(SOp,src,xc); + // Not gated: Wilson's spectrum has O(1) imaginary parts and a real-interval + // Chebyshev fit to 1/x degrades exponentially off the axis (Bernstein + // ellipse). Printed as the reminder of what a non-real spectrum does. + std::cout << GridLogMessage << " (info) ChebyshevNonHermitian on the COMPLEX-spectrum Wilson op: |r|/|r0| = " << rc + << " (GCR " << rg << ") -- expected poor; PVdagM smoother ops have real coefficients" << std::endl; + + ////////////////////////////////////////////////////////////////////// + // legacy HermOp Chebyshev smoother, both constructor forms, on MdagM + ////////////////////////////////////////////////////////////////////// + { + MdagMLinearOperator HermOp(Dw); + LatticeFermionD hsrc(UGrid), hx(UGrid), hr(UGrid); gaussian(RNG4,hsrc); + RealD hhi = 0.0; + { LatticeFermionD v(UGrid), Av(UGrid); gaussian(RNG4,v); RealD n=std::sqrt(norm2(v)); v=v*(1.0/n); + for(int i=0;i<40;i++){ HermOp.HermOp(v,Av); hhi=std::sqrt(norm2(Av)); v=Av*(1.0/hhi); } } + ChebyshevSmoother S4(0.5, 1.05*hhi, 12, HermOp); + ChebyshevSmoother S5(0.5, 1.05*hhi, 12, HermOp, Dw); // historical 5-arg form + S4(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc; + RealD r4 = std::sqrt(norm2(hr)/norm2(hsrc)); + S5(hsrc,hx); HermOp.HermOp(hx,hr); hr = hr - hsrc; + RealD r5 = std::sqrt(norm2(hr)/norm2(hsrc)); + Report("T4 legacy ChebyshevSmoother (4- and 5-arg) reduces MdagM residual", r4 < 0.5 && r5 == r4, + "|r|/|r0| = "+std::to_string(r4)+" / "+std::to_string(r5)); + + // T3: the Op()-based Clenshaw on a REAL-spectrum operator must reproduce + // the legacy HermOp smoother: same coefficients, same recurrence. + struct HermAsOp : public LinearOperatorBase { + LinearOperatorBase &H; + HermAsOp(LinearOperatorBase &h) : H(h) {} + void OpDiag (const LatticeFermionD &in, LatticeFermionD &out) { GRID_ASSERT(0); } + void OpDir (const LatticeFermionD &in, LatticeFermionD &out,int dir,int disp) { GRID_ASSERT(0); } + void OpDirAll(const LatticeFermionD &in, std::vector &out) { GRID_ASSERT(0); } + void Op (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); } + void AdjOp (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); } + void HermOpAndNorm(const LatticeFermionD &in, LatticeFermionD &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } + void HermOp (const LatticeFermionD &in, LatticeFermionD &out) { H.HermOp(in,out); } + } HOp(HermOp); + ChebyshevNonHermitianSmoother C3(0.5, 1.05*hhi, 12, HOp); + LatticeFermionD hx3(UGrid), dd(UGrid); + C3(hsrc,hx3); HermOp.HermOp(hx3,hr); hr = hr - hsrc; + RealD r3 = std::sqrt(norm2(hr)/norm2(hsrc)); + S4(hsrc,hx); dd = hx - hx3; + RealD reldiff = std::sqrt(norm2(dd)/norm2(hx)); + Report("T3 ChebyshevNonHermitian(Op=HermOp) == legacy ChebyshevSmoother", r3 < 0.5 && reldiff < 1.0e-12, + "|r|/|r0| = "+std::to_string(r3)+" rel diff of solutions "+std::to_string(reldiff)); + } + + ////////////////////////////////////////////////////////////////////// + // determinism + ////////////////////////////////////////////////////////////////////// + Replay(src,xr2); + d = xr - xr2; + Report("T5 replay bitwise repeatable", norm2(d)==0.0); + + std::cout << GridLogMessage << (failures ? "Test_poly_smoother: FAILURES" : "Test_poly_smoother: ALL PASS") << std::endl; + Grid_finalize(); + return failures ? 1 : 0; +}