diff --git a/CLAUDE.md b/CLAUDE.md index ebd2bffae..f478cb6ac 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -106,7 +106,7 @@ Tests and benchmarks that need optional fermion representations are guarded by ` ### GPU acceleration and the view/memory-manager discipline ### Multigrid (`Grid/algorithms/multigrid/`) -Aggregation-based algebraic multigrid for Wilson-type fermions. Key files: `GeneralCoarsenedMatrixMultiRHSV2.h` (the multi-RHS coarse operator, "V2": coarsening by Fourier probing, batched-GEMM apply on a D+1 grid with rhs innermost and unvectorised), `MultiRHSBlockProject.h` (in `deflation/`; the batched-GEMM transfer operators), `PVdagMMultiGrid.h` (three-level non-Hermitian PVdagM chain with dense bottom) and `HDCGMultiGrid.h` (two-level Hermitian HDCG chain) with their `*Params.h` (XML-serialisable parameters), `MrhsMultiGrid.h` (mrhs V-cycle, preconditioner interface, fp64/fp32 seam), `Smoothers.h`, `MultiGridIO.h`, `Aggregates.h` (near-null vector construction), `Geometry.h` (coarse stencils: max-norm-1 boxes of Manhattan radius 1, 2 or 4), `CoarsenedMatrix.h` (the original nearest-neighbour Wilson multigrid). `deprecated/` holds the V1 coarse operators and the Aggregation-based single-RHS ADEF-2, kept only for the pre-2026 drivers. `MultiGrid.h` is the top-level include. +Aggregation-based algebraic multigrid for Wilson-type fermions. Key files: `GeneralCoarsenedMatrixMultiRHS.h` (`MultiGeneralCoarsenedOperator`, the coarse operator for any number of right-hand sides: coarsening by Fourier probing, batched-GEMM apply on a D+1 grid with rhs innermost and unvectorised), `MultiRHSBlockProject.h` (in `deflation/`; the batched-GEMM transfer operators), `PVdagMMultiGrid.h` (three-level non-Hermitian PVdagM chain with dense bottom) and `HDCGMultiGrid.h` (two-level Hermitian HDCG chain) with their `*Params.h` (XML-serialisable parameters), `MrhsMultiGrid.h` (mrhs V-cycle, preconditioner interface, fp64/fp32 seam), `Smoothers.h`, `MultiGridIO.h`, `Aggregates.h` (near-null vector construction), `Geometry.h` (coarse stencils: max-norm-1 boxes of Manhattan radius 1, 2 or 4), `CoarsenedMatrix.h` (the original nearest-neighbour Wilson multigrid). `deprecated/` holds the superseded coarse operators and the Aggregation-based single-RHS ADEF-2, with class names prefixed `Deprecated` (`DeprecatedGeneralCoarsenedMatrix`, `DeprecatedMultiGeneralCoarsenedMatrix`, `DeprecatedTwoLevelADEF2`); nothing in the library uses them, only older drivers in `examples/` and `tests/debug/`. `MultiGrid.h` is the top-level include. ### GPU acceleration @@ -174,7 +174,7 @@ Every programme is wrapped in `Grid_init(&argc, &argv)` / `Grid_finalize()` (`Gr - Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`. - **Tensor indices are positional, not labelled.** The `Grid/tensors/` arithmetic recurses structurally over the `iScalar`/`iVector`/`iMatrix` nest: each level defines only the {scalar,vector,matrix}² products at its own level, with element types resolved by automatic type deduction, so every colour/spin/lorentz combination composes from ~200 lines (versus the pre-C++11 QDP++/PETE approach of machine-generating every case). An index's meaning derives entirely from its nesting depth counted from the outside; `iScalar` is the identity/broadcast case at every level. Never insert or remove a nesting level casually — the multiplication tables contract by position. - **Multigrid coarsening deepens the tensor nest by one level.** A coarse site vector is `iVector`, and `innerProduct` on it returns `iScalar` — one level deeper than the fine block scalar. So the block-inner-product scalar type gains one `iScalar` wrapper per MG level (fine: `vTComplex`; level 2: `iScalar`; see `examples/Example_pvdagm_3level.cc`). When calling `blockInnerProduct`/`blockZAXPY`/`blockOrthogonalise` on coarse fields, the coarse scalar type must match `decltype(innerProduct(siteVector(),siteVector()))` exactly; a wrong depth fails to compile (no viable `operator=` deep in the instantiation chain) rather than mis-contracting. -- **Grids are borrowed, never owned.** `conformable` is pointer identity, so every object that interoperates must hold the *same* `GridCartesian *`; a class that minted its own grid internally could never conform with anything else. Ownership is therefore not available, and lifetime is managed by scope discipline instead of reference counting: whoever creates a grid retains it beyond every object it handed a reference to. Anything *derived* from a grid inherits this — `~PaddedCell` dereferences its `unpadded_grid`, so a `PaddedCell` cannot even be **destroyed** after its parent grid, only used. Where a consumer must let go early, it offers an explicit hand-back (`MultiGeneralCoarsenedOperatorV2::ReleaseGrid()`) to be called *before* the grid is destroyed. +- **Grids are borrowed, never owned.** `conformable` is pointer identity, so every object that interoperates must hold the *same* `GridCartesian *`; a class that minted its own grid internally could never conform with anything else. Ownership is therefore not available, and lifetime is managed by scope discipline instead of reference counting: whoever creates a grid retains it beyond every object it handed a reference to. Anything *derived* from a grid inherits this — `~PaddedCell` dereferences its `unpadded_grid`, so a `PaddedCell` cannot even be **destroyed** after its parent grid, only used. Where a consumer must let go early, it offers an explicit hand-back (`MultiGeneralCoarsenedOperator::ReleaseGrid()`) to be called *before* the grid is destroyed. - The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`. - Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI. diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 575a03a30..986a46d02 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -667,6 +667,9 @@ public: (*this)(in[i], out[i]); } } + // The identity, if a derived class says so. A caller can then use its input + // where it would otherwise have materialised the output, and skip the copy. + virtual int isTrivial(void) { return 0; }; virtual ~LinearFunction(){}; }; diff --git a/Grid/algorithms/Preconditioner.h b/Grid/algorithms/Preconditioner.h index a95dad7c5..197575cd6 100644 --- a/Grid/algorithms/Preconditioner.h +++ b/Grid/algorithms/Preconditioner.h @@ -45,6 +45,7 @@ public: virtual void operator()(const Field &src, Field & psi){ psi = src; } + virtual int isTrivial(void) { return 1; }; TrivialPrecon(void){}; }; diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h index 9ef40b086..87225701a 100644 --- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h +++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h @@ -1,6 +1,6 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid + Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h @@ -28,213 +28,48 @@ Author: Peter Boyle /* END LEGAL */ #ifndef GRID_PREC_GCR_H #define GRID_PREC_GCR_H +#include -/////////////////////////////////////////////////////////////////////////////////////////////////////// -//VPGCR Abe and Zhang, 2005. -//INTERNATIONAL JOURNAL OF NUMERICAL ANALYSIS AND MODELING -//Computing and Information Volume 2, Number 2, Pages 147-161 -//NB. Likely not original reference since they are focussing on a preconditioner variant. -// but VPGCR was nicely written up in their paper -/////////////////////////////////////////////////////////////////////////////////////////////////////// NAMESPACE_BEGIN(Grid); -#define GCRLogLevel std::cout << GridLogMessage < class PrecGeneralisedConjugateResidual : public LinearFunction { -public: + // Declaration order is initialisation order: the adaptor is complete before + // the solver binds its reference to it. + HermOpAdaptor HermLinop; + PrecGeneralisedConjugateResidualNonHermitian GCR; +public: using LinearFunction::operator(); - RealD Tolerance; - Integer MaxIterations; - int verbose; - int mmax; - int nstep; - int steps; - int level; - GridStopWatch PrecTimer; - GridStopWatch MatTimer; - GridStopWatch LinalgTimer; - LinearFunction &Preconditioner; - LinearOperatorBase &Linop; + PrecGeneralisedConjugateResidual(RealD tol,Integer maxit,LinearOperatorBase &Linop, + LinearFunction &Prec,int mmax,int nstep) + : HermLinop(Linop), + GCR(tol,maxit,HermLinop,Prec,mmax,nstep) {}; - void Level(int lv) { level=lv; }; + void operator() (const Field &src, Field &psi) { GCR(src,psi); }; - PrecGeneralisedConjugateResidual(RealD tol,Integer maxit,LinearOperatorBase &_Linop,LinearFunction &Prec,int _mmax,int _nstep) : - Tolerance(tol), - MaxIterations(maxit), - Linop(_Linop), - Preconditioner(Prec), - mmax(_mmax), - nstep(_nstep) - { - level=1; - verbose=1; - }; - - void operator() (const Field &src, Field &psi){ - - psi=Zero(); - RealD cp, ssq,rsq; - ssq=norm2(src); - rsq=Tolerance*Tolerance*ssq; - - Field r(src.Grid()); - - PrecTimer.Reset(); - MatTimer.Reset(); - LinalgTimer.Reset(); - - GridStopWatch SolverTimer; - SolverTimer.Start(); - - steps=0; - for(int k=0;k q(mmax,grid); - std::vector p(mmax,grid); - std::vector qq(mmax); - - GCRLogLevel<< "PGCR nStep("<(mmax-1))?(mmax-1):(kp); // if more than mmax done, we orthog all mmax history. - for(int back=0;back=0); - - b=-real(innerProduct(q[peri_back],Az))/qq[peri_back]; - p[peri_kp]=p[peri_kp]+b*p[peri_back]; - q[peri_kp]=q[peri_kp]+b*q[peri_back]; - - } - qq[peri_kp]=norm2(q[peri_kp]); // could use axpy_norm - LinalgTimer.Stop(); - } - GRID_ASSERT(0); // never reached - return cp; - } + void Name(std::string n) { GCR.Name(n); }; + void Level(int n) { GCR.Level(n); }; + void SetZeroGuess(int z) { GCR.SetZeroGuess(z); }; + int Steps(void) const { return GCR.Steps(); }; + void LogCoefficients(int l) { GCR.LogCoefficients(l); }; + void SetCoefficientRecorder(GCRCoefficients *r){ GCR.SetCoefficientRecorder(r); }; + void ReleaseHistory(void) { GCR.ReleaseHistory(); }; }; + NAMESPACE_END(Grid); -#undef GCRLogLevel #endif diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h index 1feef6691..e781b78c7 100644 --- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h +++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h @@ -43,8 +43,51 @@ NAMESPACE_BEGIN(Grid); template class PrecGeneralisedConjugateResidualNonHermitian : public LinearFunction { -public: +public: using LinearFunction::operator(); + + PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase &_Linop,LinearFunction &Prec,int _mmax,int _nstep) : + Tolerance(tol), + MaxIterations(maxit), + mmax(_mmax), + nstep(_nstep), + Preconditioner(Prec), + Linop(_Linop) + { + Level(1); + verbose=1; + trivial_prec = Preconditioner.isTrivial(); + }; + + // The name is also the trace range's, so a profile separates the four + // instances (Fouter, Fsmoother, Couter, Csmoother) that otherwise nest + // indistinguishably as one shared range. + void Name(std::string _name) { name = _name; trace_step = name + " PGCR_step"; }; + + void Level(int n) { Name("Level " + std::to_string(n)); level = n; } + + void SetZeroGuess(int z) { ZeroGuess = z; }; + + // Steps taken by the last solve. + int Steps(void) const { return steps; }; + + // Coefficient logging: one line per step with the step length a_k and the + // orthogonalisation coefficients b_{k,j}. These are the data from which a + // FIXED polynomial smoother can be harvested: if they are stable from call + // to call, the adaptive GCR can be replaced by a stationary p(A) with the + // same applies and no reductions. Off by default; boss rank prints. + void LogCoefficients(int l) { LogCoeffs = l; }; + + // Optional recorder of the per-step coefficients (means over calls), for + // replay by GCRReplaySmoother (Smoothers.h). Records only; no effect on + // the iteration. + 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; }; + +private: RealD Tolerance; RealD SSQ; Integer MaxIterations; @@ -65,45 +108,14 @@ public: std::vector p; std::vector qq; int FirstCycle = 0; + int LogCoeffs = 0; + int trivial_prec; // Preconditioner is the identity: p == r, copies fold away + GCRCoefficients *Recorder = nullptr; LinearFunction &Preconditioner; LinearOperatorBase &Linop; - // The name is also the trace range's, so a profile separates the four - // instances (Fouter, Fsmoother, Couter, Csmoother) that otherwise nest - // indistinguishably as one shared range. - void Name(std::string _name) { name = _name; trace_step = name + " PGCR_step"; }; - - void Level(int n) { Name("Level " + std::to_string(n)); level = n; } - - void SetZeroGuess(int z) { ZeroGuess = z; }; - // Coefficient logging: one line per step with the step length a_k and the - // orthogonalisation coefficients b_{k,j}. These are the data from which a - // FIXED polynomial smoother can be harvested: if they are stable from call - // to call, the adaptive GCR can be replaced by a stationary p(A) with the - // same applies and no reductions. Off by default; boss rank prints. - int LogCoeffs = 0; - void LogCoefficients(int l) { LogCoeffs = l; }; - // Optional recorder of the per-step coefficients (means over calls), for - // replay by GCRReplaySmoother (Smoothers.h). Records only; no effect on - // 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), - MaxIterations(maxit), - Linop(_Linop), - Preconditioner(Prec), - mmax(_mmax), - nstep(_nstep) - { - Level(1); - verbose=1; - }; +public: void operator() (const Field &src, Field &psi){ @@ -151,6 +163,7 @@ public: // GRID_ASSERT(0); } +private: RealD GCRnStep(const Field &src, Field &psi,RealD rsq){ RealD cp; @@ -189,18 +202,20 @@ public: // contract (enforced here), so r0 = src exactly; skip the apply. // Restart cycles (psi!=0) always do the full computation. ////////////////////////////////// - if (ZeroGuess && FirstCycle) { - psi = Zero(); - LinalgTimer.Start(); - r = src; - LinalgTimer.Stop(); - } else { + // The residual is READ through rp and only materialised in r when it is + // first written, at the k=0 update. On a zero-guess first cycle r0 IS src, + // so no copy is made; a restart cycle computes r and rp points at it. + const Field *rp = &src; + GRID_ASSERT(nstep>=1); // psi is first written at k=0; see zero_start below + int zero_start = (ZeroGuess && FirstCycle); + if ( !zero_start ) { MatTimer.Start(); Linop.Op(psi,Az); MatTimer.Stop(); LinalgTimer.Start(); r=src-Az; LinalgTimer.Stop(); + rp = &r; } FirstCycle=0; @@ -208,9 +223,11 @@ public: // p = Prec(r) ///////////////////// - // p[0] = Prec(r), q[0] = A p[0], written straight into the history slots + // p[0] = Prec(r), q[0] = A p[0], written straight into the history slots. + // This copy is forced even for a trivial preconditioner: p[0] is history and + // must survive while the residual evolves. PrecTimer.Start(); - Preconditioner(r,p[0]); + Preconditioner(*rp,p[0]); PrecTimer.Stop(); MatTimer.Start(); @@ -221,7 +238,7 @@ public: qq[0]= norm2(q[0]); - cp =norm2(r); + cp =norm2(*rp); LinalgTimer.Stop(); GCRLogLevel<< "PGCR true residual "<< sqrt(cp/SSQ) <RecordA(k,a); - axpy(psi,a,p[peri_k],psi); + // On a zero-guess first step psi is still unwritten, and the update would + // be a*p added to zero: assign it instead, so psi is never zeroed. + if ( zero_start && (k==0) ) psi = a*p[peri_k]; + else axpy(psi,a,p[peri_k],psi); - cp = axpy_norm(r,-a,q[peri_k],r); +#ifdef GRID_GCR_RESIDUAL_RECURRENCE + // |r - a q|^2 = |r|^2 - ||^2/|q|^2 with a = /|q|^2, so the new + // residual norm follows from scalars already in hand and the reduction in + // axpy_norm is not needed. cp then DRIFTS rather than being measured; the + // true residual is recomputed at restart and at convergence, which bounds + // the exposure. The drift only matters for an instance driven to a tight + // tolerance -- the outer solver -- and there it is the same computed-vs-true + // residual gap that CG lives with; the smoother and coarse instances run at + // fixed work or loose tolerance and are indifferent to it. + // Explicit re/im: ComplexD is thrust::complex under HIP. + axpy(r,-a,q[peri_k],*rp); + cp = cp - (real(rq)*real(rq)+imag(rq)*imag(rq))/qq[peri_k]; + if ( cp < 0.0 ) cp = 0.0; // drift can undershoot; keeps the log finite +#else + cp = axpy_norm(r,-a,q[peri_k],*rp); +#endif + rp = &r; // the residual now lives in r; src is untouched from here on LinalgTimer.Stop(); if ( LogCoeffs ) { GCRLogLevel<<"coeff["< &rowStart, DenseInverseScalar *rows1d, int64_t myrows, BlockCyclicMatrix &A) - { Redistribute(+1, grid, rowStart, rows1d, myrows, A); } + { + Redistribute(+1, grid, rowStart, rows1d, myrows, A); + } static void CyclicToRows(GridBase *grid, const std::vector &rowStart, BlockCyclicMatrix &A, DenseInverseScalar *rows1d, int64_t myrows) - { Redistribute(-1, grid, rowStart, rows1d, myrows, A); } + { + Redistribute(-1, grid, rowStart, rows1d, myrows, A); + } }; NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h index 718f4370a..c436d6b9e 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h +++ b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h @@ -215,7 +215,11 @@ public: // SendToRecvFrom with root (symmetric byte count: the reverse direction // carries a same-size dummy -- a leaf-local cost, accepted for simplicity). /////////////////////////////////////////////////////////////////////////// - static int64_t FirstBlock(int64_t b0, int p, int Pg){ int64_t r = ((b0 % Pg) <= p) ? b0 - (b0 % Pg) + p : b0 - (b0 % Pg) + Pg + p; return r; } + static int64_t FirstBlock(int64_t b0, int p, int Pg) + { + int64_t r = ((b0 % Pg) <= p) ? b0 - (b0 % Pg) + p : b0 - (b0 % Pg) + Pg + p; + return r; + } void BigLeaf(BlockCyclicMatrix &A, int64_t b0, int64_t b1) { GRID_TRACE("SchurBigLeaf"); diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index 88c3697d9..4cdb74644 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -38,7 +38,7 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////////////////////////// // DenseCoarseMatrix: a coarsened operator treated as a DENSE matrix -- explicit, -// row-distributed A^{-1} of a GeneralCoarsenedMatrix. +// row-distributed A^{-1} of a coarsened operator. // // - Stencil -> dense DIRECT IMPORT. The coarse operator IS the dense matrix // unrolled: Dense[(s,a),(s+shift_p,b)] += A[p][s]_{a,b}. Rows of my sites are @@ -176,6 +176,7 @@ public: { double t0 = usecond(); { + CheckMrhsSlices(Op); // the batched apply keeps its slices apart ImportDense(Op); // slab <- my rows of A (LOCAL, no comms) ImportCertificate(Op); // dense apply == Op.M, before inversion InvertDense(Op); // slab <- my rows of A^{-1} @@ -274,32 +275,67 @@ public: GRID_ASSERT(mgrid->_ndimension == nd+1); int nr = mgrid->_fdimensions[0]; + // Every slice carries the same vector; slice 0 is the answer. The + // slices are not a check here -- see CheckMrhsSlices, which does that + // once on a field chosen for the purpose. + Field min(mgrid), mout(mgrid); + for(int r=0;r + void CheckMrhsSlices(CoarseOp &Op) + { + if ( Op.Grid() == grid ) return; // single rhs: nothing to mix + + GridBase *mgrid = Op.Grid(); + GRID_ASSERT(mgrid->_ndimension == nd+1); + int nr = mgrid->_fdimensions[0]; + if ( nr < 2 ) return; + + Field in(grid); + GridParallelRNG rng(grid); rng.SeedFixedIntegers(std::vector({7,8,9,10})); + random(rng,in); + Field min(mgrid), mout(mgrid); for(int r=0;r= otol ) { - std::cout << GridLogMessage << "DenseCoarseMatrix: oracle rhs "< - - 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 */ -#pragma once - - -NAMESPACE_BEGIN(Grid); - - -// Fine Object == (per site) type of fine field -// nbasis == number of deflation vectors -template -class MultiGeneralCoarsenedOperatorV2 : public SparseMatrixBase > > { -public: - typedef typename CComplex::scalar_object SComplex; - // The BLAS scalar of this level follows the coefficient precision: - // ComplexF for an fp32 coarse space, ComplexD for fp64. - typedef typename GridTypeMapper::scalar_type CoarseBLASScalar; - typedef MultiGeneralCoarsenedOperatorV2 MultiGeneralCoarseOp; - - typedef iVector siteVector; - typedef iMatrix siteMatrix; - typedef iVector calcVector; - typedef iMatrix calcMatrix; - typedef Lattice > CoarseComplexField; - typedef Lattice CoarseVector; - typedef Lattice > CoarseMatrix; - typedef iMatrix Cobj; - typedef iVector Cvec; - typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field - typedef Lattice FineField; - typedef CoarseVector Field; - - // Block operations on the fine vectors carry the fine layout, which need - // not be the coarse one - typedef decltype(innerProduct(Fobj(),Fobj())) FineInner; - typedef Lattice FineComplexField; - typedef Lattice BlockComplexField; - - //////////////////// - // Data members - // - // Nrhs independent: the D dimensional coarse grid, the geometry, the padded - // cell that supplies the stencil grid, the stencil, and the matrix elements. - // - // Nrhs dependent: the D+1 grid, its padded cell, and the BLAS B/C buffers - // with their pointer tables. Owned by SetNRHS(). - //////////////////// - GridCartesian * _CoarseGrid; // D dimensional - NonLocalStencilGeometry geom; - NonLocalStencilGeometry geom_srhs; - PaddedCell CellD; // D dimensional, supplies stencil grid - GeneralLocalStencil Stencil; // D dimensional - - int _Nrhs; - GridCartesian * _CoarseGridMulti; // D+1 dimensional, SetNRHS - PaddedCell * CellMulti; // D+1 dimensional, SetNRHS - - deviceVector BLAS_B; - deviceVector BLAS_C; - std::vector > BLAS_A; - - std::vector > BLAS_AP; - std::vector > BLAS_BP; - deviceVector BLAS_CP; - - /////////////////////////////////////////////////////////////////////////// - // Stencil legs carried in the BATCH dimension, LegGroup at a time. - // - // One call per stencil point batches over the local coarse volume alone, - // which at the production point is 1024 sites. That is too few for the - // single-precision kernel: it then runs at half the bandwidth the double - // one reaches and takes the same time, so an fp32 coarse space buys - // nothing in the operator. The same shape at four times the batch runs - // 1.7x faster in fp32 than fp64, so the fix is more batch, not fewer - // bytes. Grouping G legs multiplies the batch by G and divides the launch - // count by G, at the price of G partial outputs summed at the end. - // - // G must divide npoint; 3 divides both the 33 point (next-to-nearest) and - // 81 point (full box) stencils. G=1 is the ungrouped form. The only - // numerical difference is the ORDER of the sum over stencil points. - /////////////////////////////////////////////////////////////////////////// - // The group need not divide npoint: the last call carries the remainder, - // with its own smaller tables. The batch a call presents is LegGroup times - // the local coarse volume, so the group that saturates the kernel depends - // on the decomposition, and the default is chosen from it rather than - // fixed. TargetBatch is where the measured single-precision curve flattens - // (530 GB/s at 1024, 695 at 9216, 728 at 33792 for 60x12x60 on MI250X). - // - // 9216 is nine times the production local coarse volume, so the default - // lands on nine legs per call: 81 = 9x9 for the full-box stencil and - // 33 = 9+9+9+6 for the next-to-nearest one, which is the grouping the - // earlier HDCG coarse operator used. - static const int TargetBatch = 9216; - - int LegGroup = 1; - std::vector > BLAS_APg; // per group - std::vector > BLAS_BPg; - deviceVector BLAS_CPg; // LegGroup*sites - deviceVector BLAS_CPr; // remainder*sites - deviceVector BLAS_Cg; // LegGroup partials - - int LegGroups(void) const { return geom.npoint/LegGroup; } // full - int LegRemainder(void) const { return geom.npoint%LegGroup; } - - /////////////////////////////////////////////////////////////////////////// - // Matrix pointer table for the grouped call. Group g holds legs - // g*LegGroup .. and the batch index runs site fastest within a leg. The - // last group is short when LegGroup does not divide npoint. - /////////////////////////////////////////////////////////////////////////// - void BuildGroupedA(void) - { - int32_t sites = _CoarseGrid->lSites(); - int ng = LegGroups() + (LegRemainder()?1:0); - BLAS_APg.resize(ng); - for(int g=0;g=1 ); - GRID_ASSERT( G<=geom.npoint ); // every partial must get an initialising call - LegGroup = G; - BuildGroupedA(); - if ( _CoarseGridMulti ) SetGrid(_CoarseGridMulti); // rebuild the rest - } - - /////////////////////// - // Interface - /////////////////////// - GridBase * Grid(void) { CheckGridSet(); return _CoarseGridMulti; }; - GridCartesian * CoarseGrid(void) { CheckGridSet(); return _CoarseGridMulti; }; - GridCartesian * CoarseGridD(void) { return _CoarseGrid; }; // lower dimensional grid - int Nrhs(void) { CheckGridSet(); return _Nrhs; }; - - void CheckGridSet(void) - { - if ( _CoarseGridMulti == nullptr ) { - std::cout << GridLogError - << "MultiGeneralCoarsenedOperatorV2: the multiRHS grid has not been set." - << std::endl; - std::cout << GridLogError - << " Call SetGrid(CoarseGridMulti) with the D+1 dimensional grid your" - << std::endl; - std::cout << GridLogError - << " coarse vectors live on, before Grid(), Nrhs() or M()." - << std::endl; - GRID_ASSERT(_CoarseGridMulti != nullptr); - } - } - - ////////////////////////////////////////////////////////////////////////// - // Bilingual accessors, matching GeneralCoarsenedMatrix. Note Grid() is the - // D+1 multiRHS grid here, so a consumer wanting the space the elements live - // on must ask for CoarseGridD(). - ////////////////////////////////////////////////////////////////////////// - // Resident device memory this operator holds (matrix elements and the - // Nrhs-dependent BLAS buffers). Not evictable. - uint64_t DeviceBytes(void) - { - uint64_t b = (uint64_t)(BLAS_B.capacity()+BLAS_C.capacity()+BLAS_Cg.capacity())*sizeof(calcVector); - for(int p=0;p<(int)BLAS_A.size();p++) b += (uint64_t)BLAS_A[p].capacity()*sizeof(calcMatrix); - return b; - } - - NonLocalStencilGeometry & Geometry(void) { return geom_srhs; }; - void ExtractMatrix(int p,CoarseMatrix &A) { BLAStoGrid(A,BLAS_A[p]); }; - - // I/O on the operator matrices, via the BLAS layout array. The parameter is - // a vector over the geometry points; the body indexes A[p]. - void SetMatrix (int p,std::vector & A) - { - GRID_ASSERT(A.size()==geom_srhs.npoint); - GridtoBLAS(A[p],BLAS_A[p]); - } - void GetMatrix (int p,std::vector & A) - { - GRID_ASSERT(A.size()==geom_srhs.npoint); - BLAStoGrid(A[p],BLAS_A[p]); - } - - /////////////////////////////////////////////////////////////////////////// - // Constructor takes the D dimensional coarse grid. Everything built here - // is independent of Nrhs, in particular the matrix elements, which must - // survive a change of Nrhs untouched. - /////////////////////////////////////////////////////////////////////////// - MultiGeneralCoarsenedOperatorV2(NonLocalStencilGeometry &_geom,GridCartesian *CoarseGrid) : - _CoarseGrid(CoarseGrid), - geom_srhs(_geom), - geom(CoarseGrid,_geom.hops,_geom.skip), - CellD(geom.Depth(),CoarseGrid), - Stencil(CellD.grids.back(),geom.shifts), // D dimensional padded cell stencil - _Nrhs(-1), - _CoarseGridMulti(nullptr), - CellMulti(nullptr) - { - int32_t unpadded_sites = _CoarseGrid->lSites(); - - ///////////////////////////////////////////////// - // Matrix elements and their pointer table - ///////////////////////////////////////////////// - BLAS_A.resize(geom.npoint); - BLAS_AP.resize(geom.npoint); - for(int p=0;p geom.npoint ) LegGroup = geom.npoint; - BuildGroupedA(); - std::cout << GridLogMessage << "MultiGeneralCoarsenedOperatorV2: stencil " - << geom.npoint << " points, local coarse volume " << unpadded_sites - << ", legs per GEMM " << LegGroup - << " (batch " << LegGroup*unpadded_sites << ", " - << LegGroups() << " full call(s)" - << (LegRemainder() ? " + a remainder of "+std::to_string(LegRemainder()) : "") - << ")" << std::endl; - } - - virtual ~MultiGeneralCoarsenedOperatorV2() - { - ReleaseGrid(); - } - - /////////////////////////////////////////////////////////////////////////// - // Free everything SetGrid allocated. The D+1 grid is borrowed from the - // caller and is never deleted here. Safe to call repeatedly and before - // the destructor. - /////////////////////////////////////////////////////////////////////////// - void ReleaseGrid(void) - { - if ( CellMulti != nullptr ) { delete CellMulti; CellMulti = nullptr; } - - _CoarseGridMulti = nullptr; // borrowed, not owned - _Nrhs = -1; - - BLAS_B.resize(0); - BLAS_C.resize(0); - BLAS_Cg.resize(0); - BLAS_CPg.resize(0); - BLAS_CPr.resize(0); - for(int g=0;g<(int)BLAS_BPg.size();g++){ BLAS_BPg[g].resize(0); } - BLAS_BPg.resize(0); - for(int p=0;p_ndimension; - - GRID_ASSERT(CoarseGridMulti->_ndimension == nd+1); - GRID_ASSERT(CoarseGridMulti->_processors[0] == 1); // rhs is not distributed - for(int d=0;d_fdimensions[d+1] == _CoarseGrid->_fdimensions[d]); - GRID_ASSERT(CoarseGridMulti->_processors [d+1] == _CoarseGrid->_processors [d]); - GRID_ASSERT(CoarseGridMulti->_simd_layout[d+1] == _CoarseGrid->_simd_layout[d]); - } - - _CoarseGridMulti = CoarseGridMulti; - _Nrhs = CoarseGridMulti->_fdimensions[0]; - GRID_ASSERT(_Nrhs>=1); - - int nrhs = _Nrhs; - - CellMulti = new PaddedCell(geom.Depth(),_CoarseGridMulti); - - int32_t padded_sites = CellD.grids.back()->lSites(); // D dimensional - int32_t unpadded_sites = _CoarseGrid->lSites(); // D dimensional - - // The neighbour offset multiplication by nrhs is exact only if the D+1 - // padded volume is nrhs copies of the D dimensional one. Check it. - GRID_ASSERT(CellMulti->grids.back()->lSites() == nrhs*padded_sites); - GRID_ASSERT(_CoarseGridMulti->lSites() == nrhs*unpadded_sites); - - ///////////////////////////////////////////////// - // Device data vector storage - ///////////////////////////////////////////////// - BLAS_B.resize(nrhs *padded_sites); // includes ghost zone - BLAS_C.resize(nrhs *unpadded_sites); // no ghost zone - BLAS_BP.resize(geom.npoint); - for(int p=0;p lSite, D dim - nbr = nbr*nrhs; // D -> D+1, rhs innermost - GRID_ASSERT(nbr void GridtoBLAS(const Lattice &from,deviceVector &to) - { - typedef typename vobj::scalar_object sobj; - typedef typename vobj::scalar_type scalar_type; - typedef typename vobj::vector_type vector_type; - - GridBase *Fg = from.Grid(); - GRID_ASSERT(!Fg->_isCheckerBoarded); - int nd = Fg->_ndimension; - - to.resize(Fg->lSites()); - - Coordinate LocalLatt = Fg->LocalDimensions(); - size_t nsite = 1; - for(int i=0;i_ostride; - Coordinate f_istride = Fg->_istride; - Coordinate f_rdimensions = Fg->_rdimensions; - - autoView(from_v,from,AcceleratorRead); - auto to_v = &to[0]; - - const int words=sizeof(vobj)/sizeof(vector_type); - accelerator_for(idx,nsite,1,{ - - Coordinate from_coor, base; - Lexicographic::CoorFromIndex(base,idx,LocalLatt); - for(int i=0;i void BLAStoGrid(Lattice &grid,deviceVector &in) - { - typedef typename vobj::scalar_object sobj; - typedef typename vobj::scalar_type scalar_type; - typedef typename vobj::vector_type vector_type; - - GridBase *Tg = grid.Grid(); - GRID_ASSERT(!Tg->_isCheckerBoarded); - int nd = Tg->_ndimension; - - GRID_ASSERT(in.size()==Tg->lSites()); - - Coordinate LocalLatt = Tg->LocalDimensions(); - size_t nsite = 1; - for(int i=0;i_ostride; - Coordinate t_istride = Tg->_istride; - Coordinate t_rdimensions = Tg->_rdimensions; - - autoView(to_v,grid,AcceleratorWrite); - auto from_v = &in[0]; - - const int words=sizeof(vobj)/sizeof(vector_type); - accelerator_for(idx,nsite,1,{ - - Coordinate to_coor, base; - Lexicographic::CoorFromIndex(base,idx,LocalLatt); - for(int i=0;i - // = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} - // = M_{kl} A_ji^{b.b+l} - // - // Where q_k = delta_k . (2*M_PI/global_nb[mu]) - // Then A{ji}^{b,b+l} = M^{-1}_{lm} ComputeProj_{m,b,i,j} - /////////////////////////////////////////////////////////////////////////// - /////////////////////////////////////////////////////////////////////////// - // Probe momenta. The stencil DISPLACEMENTS are the +-1 box; the momenta - // used to separate them are free, and are chosen here to span the - // Brillouin zone: k_mu = K_mu * shift_mu with K_mu ~ L_mu/3, so a phase - // is ~2pi/3 per unit displacement rather than 2pi/L_mu. - // - // This is what conditions the extraction. With K_mu = 1 every entry of - // the phase matrix tends to 1 as the coarse lattice grows, the matrix - // tends to rank one, and its inverse amplifies any error in the measured - // projections: at 24.24.16.32 the condition number is 2.7e5, so fp32 - // projections give coarse matrix elements wrong by ~1e-3. Spread momenta - // bring it to ~20, independent of the lattice size. - // - // K_mu is stepped away from L_mu/2, where +k and -k alias into the same - // momentum and the matrix is singular. - /////////////////////////////////////////////////////////////////////////// - void CoarsenMomenta(GridBase *CoarseGrid,Coordinate &K) - { - Coordinate clatt = CoarseGrid->GlobalDimensions(); - int Nd = CoarseGrid->Nd(); - K.resize(Nd); - for(int mu=0;mu=3) ? (int)std::lround(L/3.0) : 1; - if ( k < 1 ) k = 1; - if ( (L>2) && ((2*k)%L == 0) ) k = k-1; // +k and -k must differ - if ( k < 1 ) k = 1; - K[mu] = k; - } - } - - void CoarsenFourierMatrix(GridBase *CoarseGrid,Eigen::MatrixXcd &invMkl) - { - const int npoint = geom_srhs.npoint; - Coordinate clatt = CoarseGrid->GlobalDimensions(); - int Nd = CoarseGrid->Nd(); - Coordinate K; CoarsenMomenta(CoarseGrid,K); - - Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint); - ComplexD ci(0.0,1.0); - for(int k=0;k svd(Mkl); - RealD cond = svd.singularValues()(0)/svd.singularValues()(npoint-1); - std::cout << GridLogMessage << "CoarsenOperator: probe momenta "<_ndimension; - latt.resize(nd); simd.resize(nd); mpi.resize(nd); - for(int d=0;d_fdimensions[d]; - simd[d] = grid->_simd_layout[d]; - mpi [d] = CoarseGrid->_processors[d]; - } - } - - // D+1 coarse grid holding the batch, rhs innermost and unvectorised - void CoarsenBatchGridLayout(GridBase *CoarseGrid,int batch, - Coordinate &latt,Coordinate &simd,Coordinate &mpi) - { - latt.resize(1,batch); simd.resize(1,1); mpi.resize(1,1); - latt[0]=batch; simd[0]=1; mpi[0]=1; - for(int d=0;d_ndimension;d++){ - latt.push_back(CoarseGrid->_fdimensions[d]); - simd.push_back(CoarseGrid->_simd_layout[d]); - mpi .push_back(CoarseGrid->_processors[d]); - } - } - - /////////////////////////////////////////////////////////////////////////// - // The Fourier inverse needs the phase in the coarse layout and the basis - // phasing needs it in the fine layout; each is built from its own - // coordinates rather than transferred. - /////////////////////////////////////////////////////////////////////////// - void CoarsenPhases(GridBase *grid,GridBase *CoarseGrid,GridCartesian *BlockGrid, - std::vector &pha, - std::vector &phaF) - { - const int npoint = geom_srhs.npoint; - Coordinate clatt = CoarseGrid->GlobalDimensions(); - int Nd = CoarseGrid->Nd(); - Coordinate K; CoarsenMomenta(CoarseGrid,K); - ComplexD ci(0.0,1.0); - - // The fine-side scratch carries the FINE level's precision, which need - // not be the coarse one (fp64 fine, fp32 coarse at L1). - typedef typename GridTypeMapper::scalar_type FineScalar; - FineComplexField one(grid); one=FineScalar(1.0); - FineComplexField zz(grid); zz = Zero(); - BlockComplexField pha_blk (BlockGrid); - BlockComplexField blk_coor(BlockGrid); - - for(int p=0;p &pha, - CoarseComplexField &phaB, - CoarseVector &TmpProj, - std::vector &_A, - GridBase *CoarseGrid) - { - typedef typename CComplex::scalar_type SComplex; - const int npoint = geom_srhs.npoint; - - for(int b=0;boSites(); - for(int k=0;k > &linop, - GridCartesian *FineGridMulti, - std::vector &Subspace, - GridBase *CoarseGrid) - { - MultiRHSBlockProject > Projector; - CoarsenOperator(linop,FineGridMulti,Subspace,CoarseGrid,Projector); - } - template - void CoarsenOperator(LinearOperatorBase > &linop, - GridCartesian *FineGridMulti, - std::vector &Subspace, - GridBase *CoarseGrid, - Projector_t &Projector) - { - RealD tproj=0.0, tmat=0.0, tphase=0.0, tphaseBZ=0.0, tslice=0.0, tinv=0.0; - - std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs (multiRHS fine operator)"<< std::endl; - - GRID_ASSERT(Subspace.size()==nbasis); - GridBase *grid = Subspace[0].Grid(); - - GRID_ASSERT(FineGridMulti->_ndimension == grid->_ndimension+1); - GRID_ASSERT(FineGridMulti->_processors[0] == 1); - for(int d=0;d_ndimension;d++){ - GRID_ASSERT(FineGridMulti->_fdimensions[d+1] == grid->_fdimensions[d]); - GRID_ASSERT(FineGridMulti->_processors [d+1] == grid->_processors [d]); - } - int batch = FineGridMulti->_fdimensions[0]; - - Coordinate blatt,bsimd,bmpi; - CoarsenBlockGridLayout(grid,CoarseGrid,blatt,bsimd,bmpi); - GridCartesian BlockGrid(blatt,bsimd,bmpi); - - BlockComplexField InnerProd(&BlockGrid); - blockOrthogonalise(InnerProd,Subspace); - - // the caller owns it; import the basis we have just orthonormalised - Projector.Allocate(nbasis,grid,CoarseGrid); - Projector.ImportBasis(Subspace); - - const int npoint = geom_srhs.npoint; - - Eigen::MatrixXcd invMkl; - CoarsenFourierMatrix(CoarseGrid,invMkl); - - FineField phaV(grid); - std::vector phaF(npoint,grid); - std::vector pha (npoint,CoarseGrid); - - tphase=-usecond(); - CoarsenPhases(grid,CoarseGrid,&BlockGrid,pha,phaF); - tphase+=usecond(); - - std::vector _A; - _A.resize(npoint,CoarseGrid); - for(int k=0;k > &linop, - std::vector &Subspace, - GridBase *CoarseGrid, - int batch) - { - MultiRHSBlockProject > Projector; - CoarsenOperator(linop,Subspace,CoarseGrid,batch,Projector); - } - template - void CoarsenOperator(LinearOperatorBase > &linop, - std::vector &Subspace, - GridBase *CoarseGrid, - int batch, - Projector_t &Projector) - { - RealD tproj=0.0, tmat=0.0, tphase=0.0, tphaseBZ=0.0, tslice=0.0, tinv=0.0; - - std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs (single RHS fine operator)"<< std::endl; - - GRID_ASSERT(Subspace.size()==nbasis); - GRID_ASSERT(batch>=1); - GridBase *grid = Subspace[0].Grid(); - - Coordinate blatt,bsimd,bmpi; - CoarsenBlockGridLayout(grid,CoarseGrid,blatt,bsimd,bmpi); - GridCartesian BlockGrid(blatt,bsimd,bmpi); - - BlockComplexField InnerProd(&BlockGrid); - blockOrthogonalise(InnerProd,Subspace); - - // the caller owns it; import the basis we have just orthonormalised - Projector.Allocate(nbasis,grid,CoarseGrid); - Projector.ImportBasis(Subspace); - - const int npoint = geom_srhs.npoint; - - Eigen::MatrixXcd invMkl; - CoarsenFourierMatrix(CoarseGrid,invMkl); - - FineField phaV(grid); - std::vector phaF(npoint,grid); - std::vector pha (npoint,CoarseGrid); - - tphase=-usecond(); - CoarsenPhases(grid,CoarseGrid,&BlockGrid,pha,phaF); - tphase+=usecond(); - - std::vector _A; - _A.resize(npoint,CoarseGrid); - for(int k=0;k MphaV(batch,grid); - - for(int i0=0;i0M(in,out); - } - void M (const CoarseVector &in, CoarseVector &out) - { - // std::cout << GridLogMessage << "New Mrhs coarse"<ExchangePeriodic(in); }(); //padded input - t_exch+=usecond(); - - int npoint = geom.npoint; - typedef calcMatrix* Aview; - typedef LatticeView Vview; - - const int Nsimd = CComplex::Nsimd(); - - int64_t nrhs =pin.Grid()->GlobalDimensions()[0]; - GRID_ASSERT(nrhs>=1); - - RealD flops,bytes; - int64_t osites=in.Grid()->oSites(); // unpadded - int64_t unpadded_vol = CoarseGrid()->lSites()/nrhs; - - flops = 1.0* npoint * nbasis * nbasis * 8.0 * osites * CComplex::Nsimd(); - bytes = 1.0*osites*sizeof(siteMatrix)*npoint/pin.Grid()->GlobalDimensions()[0] - + 2.0*osites*sizeof(siteVector)*npoint; - - - t_GtoB=-usecond(); - { GRID_TRACE("CoarseV2GridToBLAS"); - GridtoBLAS(pin,BLAS_B); - } - t_GtoB+=usecond(); - - GridBLAS BLAS; - - t_mult=-usecond(); - { GRID_TRACE("CoarseV2StencilGEMM"); - // The scalar type selects the GEMM: Cgemm for an fp32 coarse space, - // Zgemm for fp64. - if ( LegGroup == 1 ) { - for(int p=0;p 1 ) { GRID_TRACE("CoarseV2LegSum"); - // Sum the LegGroup partial results into the single output buffer. - // Flat over scalars: a calcVector is nbasis of them and the partials - // are contiguous, one whole C volume after another. - int64_t nscalar = (int64_t)BLAS_C.size()*nbasis; // one whole C volume - const int G = LegGroup; - CoarseBLASScalar *dst = (CoarseBLASScalar *)&BLAS_C[0]; - CoarseBLASScalar *src = (CoarseBLASScalar *)&BLAS_Cg[0]; - accelerator_for(i,nscalar,1,{ - CoarseBLASScalar sum = src[i]; - for(int l=1;l &out){assert(0);}; -}; - -NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/HDCGMultiGrid.h b/Grid/algorithms/multigrid/HDCGMultiGrid.h index 065d4bd9f..89c49a477 100644 --- a/Grid/algorithms/multigrid/HDCGMultiGrid.h +++ b/Grid/algorithms/multigrid/HDCGMultiGrid.h @@ -34,7 +34,7 @@ NAMESPACE_BEGIN(Grid); // HDCGCoarsening everything that is a function of the gauge // field: the raw near-null basis (and its // refinement), the transfer operators, the -// Galerkin coarse operator (V2, Hermitian), +// Galerkin coarse operator (Hermitian), // the coarse eigenvectors for deflation // HDCGSolver the solve chain on a borrowed coarsening: // smoother, deflated coarse solve, ADEF-2 @@ -51,7 +51,7 @@ template class HDCGCoarsening { public: typedef Lattice FineField; - typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; + typedef MultiGeneralCoarsenedOperator CoarseOperator; typedef typename CoarseOperator::CoarseVector CoarseVector; typedef typename GridTypeMapper::SinglePrecision FobjF; typedef Lattice FineFieldF; diff --git a/Grid/algorithms/multigrid/MrhsMultiGrid.h b/Grid/algorithms/multigrid/MrhsMultiGrid.h index cf7fabb6a..d678511f8 100644 --- a/Grid/algorithms/multigrid/MrhsMultiGrid.h +++ b/Grid/algorithms/multigrid/MrhsMultiGrid.h @@ -50,23 +50,74 @@ public: LinearOperatorBase &Linop; MrhsLinearFunction &Preconditioner; std::function OnStep; // called with the outer step count after every step - void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } - void Name(std::string n){ name = n; trace_op = name+" MrhsPGCR::vOp"; trace_orthog = name+" MrhsPGCR orthog"; } - void SetZeroGuess(int z){ ZeroGuess=z; } + + void Level(int lv){ + name = "Level " + std::to_string(lv); + level=lv; + } + + void Name(std::string n){ + name = n; + trace_op = name+" MrhsPGCR::vOp"; + trace_orthog = name+" MrhsPGCR orthog"; + } + + void SetZeroGuess(int z) + { + ZeroGuess=z; + } + MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase &_Linop,MrhsLinearFunction &Prec,int _mmax,int _nstep) - : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; } - static RealD vnorm2(std::vector &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; } - static ComplexD vinnerProduct(std::vector &x,std::vector &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; } - static void vaxpy(std::vector &z,ComplexD a,std::vector &x,std::vector &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); } - void vOp(std::vector &in,std::vector &out){ GRID_TRACE(trace_op.c_str()); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); } + : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep) + { + level=1; + } + + static RealD vnorm2(std::vector &x) + { + RealD s=0; + for(auto &f:x) s+=norm2(f); + return s; + } + + static ComplexD vinnerProduct(std::vector &x,std::vector &y){ + ComplexD s(0); + for(int r=0;r<(int)x.size();r++) { + s+=innerProduct(x[r],y[r]); + } + return s; + } + + static void vaxpy(std::vector &z,ComplexD a,std::vector &x,std::vector &y) + { + for(int r=0;r<(int)z.size();r++) { + axpy(z[r],a,x[r],y[r]); + } + } + + void vOp(std::vector &in,std::vector &out) + { + GRID_TRACE(trace_op.c_str()); + for(int r=0;r<(int)in.size();r++) { + Linop.Op(in[r],out[r]); + } + } void operator()(std::vector &src,std::vector &psi){ - RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid(); - ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq; + RealD cp,ssq,rsq; + int nrhs=src.size(); + GridBase *grid=src[0].Grid(); + ssq=vnorm2(src); + rsq=Tolerance*Tolerance*ssq; + std::vector r(nrhs,grid); - GridStopWatch T; T.Start(); steps=0; FirstCycle=1; + GridStopWatch T; + T.Start(); + steps=0; + FirstCycle=1; for(int k=0;k &src,std::vector &psi,RealD rsq){ - RealD cp; ComplexD a,rq; int nrhs=src.size(); GridBase *grid=src[0].Grid(); + RealD GCRnStep(std::vector &src,std::vector &psi,RealD rsq) + { + RealD cp; + ComplexD a,rq; + int nrhs=src.size(); + GridBase *grid=src[0].Grid(); + std::vector r(nrhs,grid),Az(nrhs,grid); // Az: restart residual scratch only std::vector< std::vector > q(mmax,std::vector(nrhs,grid)); std::vector< std::vector > p(mmax,std::vector(nrhs,grid)); std::vector qq(mmax); - if (ZeroGuess && FirstCycle) { for(int rr=0;rr(mmax-1))?(mmax-1):(kp); + { GRID_TRACE(trace_orthog.c_str()); // Classical Gram-Schmidt: all coefficients against the UN-updated new q @@ -110,15 +194,30 @@ public: std::vector bcoef(northog,ComplexD(0.0)), part; for(int rr=0;rr qwin(northog); - for(int back=0;back=0); qwin[back]=&q[peri_back][rr]; } + for(int back=0;back=0); + qwin[back]=&q[peri_back][rr]; + } + rankInnerProductMulti(part,qwin,q[peri_kp][rr]); + for(int back=0;backGlobalSumVector(&bcoef[0],northog); - for(int back=0;back qwin(northog), pwin(northog); - for(int back=0;back::operator(); virtual void operator()(const CoarseField &in, CoarseField &out) { + GRID_TRACE("MGCoarseVcycle"); CoarseField vec1(in.Grid()); CoarseField vec2(in.Grid()); - out = in; - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - - // restrict, through the mixed blockProject: D+1 coarse in, D+1 cc out CoarseCoarseField CCsrc(_CoarseCoarseMrhs); CoarseCoarseField CCsol(_CoarseCoarseMrhs); - _Projector.blockProject(vec1,CCsrc); - - _CoarseCoarseSolve(CCsrc,CCsol); - - _Projector.blockPromote(vec1,CCsol); - add(out,out,vec1); - - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - _CoarseSmoother(vec1,vec2); - add(out,out,vec2); + // The cycle starts from x0 = in, one unit-step Richardson iteration ahead of + // a zero start; measured, it pays for its apply by cutting Couter steps 22%. + // x0 is never materialised in out: the operator is applied to in directly and + // the copy folds into the add that lands the coarse-coarse correction. + { + GRID_TRACE("MGCoarseResidual"); + _CoarseOp.Op(in,vec1); sub(vec1,in,vec1); + } + // restrict, through the mixed blockProject: D+1 coarse in, D+1 cc out + { + GRID_TRACE("MGCCProject"); + _Projector.blockProject(vec1,CCsrc); + } + { + GRID_TRACE("MGCCSolve"); + _CoarseCoarseSolve(CCsrc,CCsol); + } + { + GRID_TRACE("MGCCPromote"); + _Projector.blockPromote(vec1,CCsol); + add(out,in,vec1); + } + { + GRID_TRACE("MGCoarseResidual2"); + _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); + } + { + GRID_TRACE("MGCoarseSmooth"); + _CoarseSmoother(vec1,vec2); + add(out,out,vec2); + } } }; @@ -211,37 +328,47 @@ public: GridBase *_CoarseGrid, *_CoarseGridMrhs; std::function SetSloppy = [](int){}; int SloppyComms = 0; // value passed to SetSloppy on entry + MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, Projector_t &Projector, LinearFunction &CoarseSolve, GridBase *CoarseGrid, GridBase *CoarseGridMrhs) : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve), _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){} - virtual void operator()(std::vector &in, std::vector &out){ + + virtual void operator()(std::vector &in, std::vector &out) + { GRID_TRACE("MGVcycle"); SetSloppy(SloppyComms); int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); - for(int r=0;r D+1 coarse, via the mixed blockProject CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); - { GRID_TRACE("MGProject"); + { + GRID_TRACE("MGProject"); _Projector.blockProject(vec1,CsrcMrhs); } CsolMrhs=Zero(); - { GRID_TRACE("MGCoarseSolve"); + { + GRID_TRACE("MGCoarseSolve"); _CoarseSolve(CsrcMrhs,CsolMrhs); } - { GRID_TRACE("MGPromote"); + { + GRID_TRACE("MGPromote"); _Projector.blockPromote(vec1,CsolMrhs); - for(int r=0;r _in_f, _out_f; + MrhsMixedPrecPreconditioner(MrhsPreconditioner &Inner, GridBase *gridD, GridBase *gridF, int nrhs) : _Inner(Inner), _ws_d2f(gridF,gridD), _ws_f2d(gridD,gridF), _gridF(gridF) {} + void Scratch(int nrhs){ if ( (int)_in_f.size() == nrhs ) return; _in_f.clear(); _in_f.reserve(nrhs); _out_f.clear(); _out_f.reserve(nrhs); for(int r=0;r &in, std::vector &out){ + + virtual void operator()(std::vector &in, std::vector &out) + { GRID_TRACE("MGPrecisionSeam"); int nrhs=in.size(); Scratch(nrhs); for(int r=0;r &x, std::vector &src){ + virtual void Vstart(std::vector &x, std::vector &src) + { GRID_TRACE("MGPrecisionSeamVstart"); int nrhs=src.size(); Scratch(nrhs); for(int r=0;r #include #include #include -#include -// DEPRECATED: the V1 coarse operators and the Aggregation-based single-RHS -// ADEF-2. Nothing in the library uses them; the PVdagM and HDCG chains are -// on V2 (PVdagMMultiGrid.h, HDCGMultiGrid.h). Kept so the pre-2026 drivers -// in tests/debug and examples still build; removing this block is the -// deletion gate for them. +#include +// DEPRECATED: the superseded coarse operators and the Aggregation-based +// single-RHS ADEF-2. Nothing in the library uses them; the PVdagM and HDCG +// chains are on MultiGeneralCoarsenedOperator (PVdagMMultiGrid.h, +// HDCGMultiGrid.h). Kept so the pre-2026 drivers in tests/debug and +// examples still build; removing this block is the deletion gate for them. #include #include #include diff --git a/Grid/algorithms/multigrid/PVdagMMultiGrid.h b/Grid/algorithms/multigrid/PVdagMMultiGrid.h index 1e646044d..4d5fb8771 100644 --- a/Grid/algorithms/multigrid/PVdagMMultiGrid.h +++ b/Grid/algorithms/multigrid/PVdagMMultiGrid.h @@ -212,14 +212,14 @@ template class PVdagMMultiGridCoarsening { public: typedef Lattice FineField; - typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; + typedef MultiGeneralCoarsenedOperator CoarseOperator; typedef typename CoarseOperator::CoarseVector CoarseVector; typedef typename CoarseVector::vector_object CoarseSiteObj; // Every coarse level carries the same site type iVector: // the coefficient scalar does not deepen with the level (the operator keeps // its own deeper scratch type for the block inner products). Levels are // told apart by their grids, not their C++ types. - typedef MultiGeneralCoarsenedOperatorV2 CoarseCoarseOperator; + typedef MultiGeneralCoarsenedOperator CoarseCoarseOperator; typedef typename CoarseCoarseOperator::CoarseVector CoarseCoarseVector; typedef DenseCoarseMatrix DenseBottom; // The fp32 fine field, derived from the fp64 one: the fp32 fine level of @@ -359,7 +359,7 @@ public: rawPsi.reserve(nbasis); for(int k=0;k LinOpCoarse(CoarseOpPV); std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L2 CoarsenOperator, batch " diff --git a/Grid/algorithms/multigrid/PVdagMOperators.h b/Grid/algorithms/multigrid/PVdagMOperators.h index dcdde3d0b..4c90734ed 100644 --- a/Grid/algorithms/multigrid/PVdagMOperators.h +++ b/Grid/algorithms/multigrid/PVdagMOperators.h @@ -36,17 +36,54 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////////// template class PVdagMLinearOperator : public LinearOperatorBase { - Matrix &_Mat; Matrix &_PV; + Matrix &_Mat; + Matrix &_PV; public: - PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {}; - void SloppyComms(int sloppy) { _Mat.SloppyComms(sloppy); _PV.SloppyComms(sloppy); } - 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){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); } - void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); } - void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); } - void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): + _Mat(Mat),_PV(PV) + {}; + void SloppyComms(int sloppy) + { + _Mat.SloppyComms(sloppy); + _PV.SloppyComms(sloppy); + } + 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) + { + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out) + { + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2) + { + HermOp(in,out); + ComplexD d=innerProduct(in,out); + n1=real(d); + n2=norm2(out); + } + void HermOp(const Field &in, Field &out) + { + Field tmp(in.Grid()); + Op(in,tmp); + AdjOp(tmp,out); + } }; ////////////////////////////////////////////////////////////////////// @@ -55,17 +92,46 @@ public: ////////////////////////////////////////////////////////////////////// template class ShiftedPVdagMLinearOperator : public LinearOperatorBase { - Matrix &_Mat; Matrix &_PV; + Matrix &_Mat; + Matrix &_PV; public: RealD shift; ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; - 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){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; } - void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,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); } + 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) + { + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + out = out + shift*in; + } + void AdjOp (const Field &in, Field &out) + { + Field tmp(in.Grid()); + _PV.M(in,tmp); + _Mat.Mdag(tmp,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); + } }; ////////////////////////////////////////////////////////////////////// @@ -73,16 +139,41 @@ public: ////////////////////////////////////////////////////////////////////// template class ShiftedLinearOperator : public LinearOperatorBase { - LinearOperatorBase &_Op; RealD shift; + LinearOperatorBase &_Op; + RealD shift; public: ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : _Op(Op), shift(_shift) {} - 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); } + 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); + } }; NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h index 2b988b64c..f03e14ed9 100644 --- a/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h +++ b/Grid/algorithms/multigrid/deprecated/GeneralCoarsenedMatrix.h @@ -38,10 +38,10 @@ NAMESPACE_BEGIN(Grid); // Fine Object == (per site) type of fine field // nbasis == number of deflation vectors template -class GeneralCoarsenedMatrix : public SparseMatrixBase > > { +class DeprecatedGeneralCoarsenedMatrix : public SparseMatrixBase > > { public: - typedef GeneralCoarsenedMatrix GeneralCoarseOp; + typedef DeprecatedGeneralCoarsenedMatrix GeneralCoarseOp; typedef iVector siteVector; typedef iMatrix siteMatrix; typedef Lattice > CoarseComplexField; @@ -103,7 +103,7 @@ public: void ProjectNearestNeighbour(RealD shift, GeneralCoarseOp &CopyMe) { int nfound=0; - std::cout << GridLogMessage <<"GeneralCoarsenedMatrix::ProjectNearestNeighbour "<< CopyMe._A[0].Grid()< -class MultiGeneralCoarsenedMatrix : public SparseMatrixBase > > { +class DeprecatedMultiGeneralCoarsenedMatrix : public SparseMatrixBase > > { public: typedef typename CComplex::scalar_object SComplex; - typedef GeneralCoarsenedMatrix GeneralCoarseOp; - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarseOp; + typedef DeprecatedGeneralCoarsenedMatrix GeneralCoarseOp; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarseOp; typedef iVector siteVector; typedef iMatrix siteMatrix; @@ -78,13 +78,33 @@ public: GridCartesian * CoarseGrid(void) { return _CoarseGridMulti; }; // this is all the linalg routines need to know ////////////////////////////////////////////////////////////////////////// - // Bilingual accessors, matching GeneralCoarsenedMatrix. Grid() here is the + // Accessors shared with DeprecatedGeneralCoarsenedMatrix. Grid() here is the // D+1 multiRHS grid and this class never holds the D dimensional one, so // ExtractMatrix writes into whatever grid the caller's lattice is on. ////////////////////////////////////////////////////////////////////////// NonLocalStencilGeometry & Geometry(void) { return geom_srhs; }; void ExtractMatrix(int p,CoarseMatrix &A) { BLAStoGrid(A,BLAS_A[p]); }; + // One stencil point in or out, matching MultiGeneralCoarsenedOperator, so a + // caller can read or write a point without knowing which of the two layouts + // (per-point here, site-major there) holds it. + void MatrixPointIn(int p,const deviceVector &src) + { + int64_t sites = BLAS_A[p].size(); + GRID_ASSERT((int64_t)src.size() == sites); + calcMatrix *dst = &BLAS_A[p][0]; + const calcMatrix *sp = &src[0]; + accelerator_for(ss,sites,1,{ dst[ss] = sp[ss]; }); + } + void MatrixPointOut(int p,deviceVector &dst) + { + int64_t sites = BLAS_A[p].size(); + dst.resize(sites); + const calcMatrix *sp = &BLAS_A[p][0]; + calcMatrix *dp = &dst[0]; + accelerator_for(ss,sites,1,{ dp[ss] = sp[ss]; }); + } + // I/O on the operator matrices, via the BLAS layout array. The parameter is // a vector over the geometry points; the body indexes A[p]. void SetMatrix (int p,std::vector & A) @@ -121,7 +141,7 @@ public: } */ - MultiGeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *CoarseGridMulti) : + DeprecatedMultiGeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *CoarseGridMulti) : _CoarseGridMulti(CoarseGridMulti), geom_srhs(_geom), geom(_CoarseGridMulti,_geom.hops,_geom.skip+1), diff --git a/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h b/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h index 52ff3e851..850daf318 100644 --- a/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h +++ b/Grid/algorithms/multigrid/deprecated/TwoLevelADEF2.h @@ -27,15 +27,15 @@ Author: Peter Boyle /* END LEGAL */ #pragma once -// DEPRECATED with the V1 coarse operators: the single-RHS ADEF-2 on an +// DEPRECATED with the other classes in deprecated/: the single-RHS ADEF-2 on an // Aggregation (ProjectToSubspace/PromoteFromSubspace) and a D-dimensional // coarse operator. The mrhs solver in AdefMrhs.h covers one right-hand -// side through its LinearFunction interface, on the V2 coarse operator. +// side through its LinearFunction interface, on MultiGeneralCoarsenedOperator. NAMESPACE_BEGIN(Grid); template -class TwoLevelADEF2 : public TwoLevelCG +class DeprecatedTwoLevelADEF2 : public TwoLevelCG { public: /////////////////////////////////////////////////////////////////////////////////// @@ -50,7 +50,7 @@ class TwoLevelADEF2 : public TwoLevelCG /////////////////////////////////////////////////////////////////////////////////// // more most opertor functions - TwoLevelADEF2(RealD tol, + DeprecatedTwoLevelADEF2(RealD tol, Integer maxit, LinearOperatorBase &FineLinop, LinearFunction &Smoother, diff --git a/Grid/communicator/Communicator_base.h b/Grid/communicator/Communicator_base.h index e98e26ec5..fdf35f2e6 100644 --- a/Grid/communicator/Communicator_base.h +++ b/Grid/communicator/Communicator_base.h @@ -133,6 +133,7 @@ public: template void GlobalSumP2P(obj &o) { + GRID_TRACE("GlobalSumP2P"); std::vector column; obj accum = o; int source,dest; diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h index 4fd86c674..7d534546f 100644 --- a/Grid/lattice/Lattice_reduction.h +++ b/Grid/lattice/Lattice_reduction.h @@ -898,7 +898,8 @@ template void axpyMultiChunk(Lattice &z,const ComplexD *b, const std::vector*> &x,int m, int do_norm, - decltype(innerProduct(vobj(),vobj())) *inner_tmp_v) + decltype(innerProduct(vobj(),vobj())) *inner_tmp_v, + const Lattice *from=nullptr) { typedef decltype(z.View(AcceleratorRead)) View; GRID_ASSERT(m>=1 && m<=B); @@ -917,13 +918,29 @@ void axpyMultiChunk(Lattice &z,const ComplexD *b, } for(int j=m;j &y = *from; // named: autoView is a macro, *p.View() misparses + autoView(y_v,y,AcceleratorRead); + autoView(z_v,z,AcceleratorWriteDiscard); + accelerator_for(ss,sites,nsimd,{ + auto acc = coalescedRead(y_v[ss]); + for(int j=0;j &z,const ComplexD *b, // (last) pass. Chunks of 16 for windows longer than 16. template RealD axpyMultiNormImpl(Lattice &z,const std::vector &b, - const std::vector*> &x,int do_norm) + const std::vector*> &x,int do_norm, + const Lattice *from=nullptr) { typedef decltype(innerProduct(vobj(),vobj())) inner_t; int m = x.size(); @@ -943,6 +961,7 @@ RealD axpyMultiNormImpl(Lattice &z,const std::vector &b, inner_t *inner_tmp_v = &inner_tmp[0]; if ( m==0 ) { + if ( from ) z = *from; if ( do_norm ) return norm2(z); return 0.0; } @@ -951,10 +970,12 @@ RealD axpyMultiNormImpl(Lattice &z,const std::vector &b, int last = (j0+mm>=m); std::vector*> sub(x.begin()+j0,x.begin()+j0+mm); int dn = do_norm && last; - if ( mm<=2 ) axpyMultiChunk<2> (z,&b[j0],sub,mm,dn,inner_tmp_v); - else if ( mm<=4 ) axpyMultiChunk<4> (z,&b[j0],sub,mm,dn,inner_tmp_v); - else if ( mm<=8 ) axpyMultiChunk<8> (z,&b[j0],sub,mm,dn,inner_tmp_v); - else axpyMultiChunk<16>(z,&b[j0],sub,mm,dn,inner_tmp_v); + // Only the first chunk is seeded from `from`; later chunks accumulate into z. + const Lattice *f = (j0==0) ? from : nullptr; + if ( mm<=2 ) axpyMultiChunk<2> (z,&b[j0],sub,mm,dn,inner_tmp_v,f); + else if ( mm<=4 ) axpyMultiChunk<4> (z,&b[j0],sub,mm,dn,inner_tmp_v,f); + else if ( mm<=8 ) axpyMultiChunk<8> (z,&b[j0],sub,mm,dn,inner_tmp_v,f); + else axpyMultiChunk<16>(z,&b[j0],sub,mm,dn,inner_tmp_v,f); } RealD nrm = 0.0; @@ -969,6 +990,15 @@ void axpyMulti(Lattice &z,const std::vector &b,const std::vector { axpyMultiNormImpl(z,b,x,0); } +// z = y + sum_j b[j] x[j]. Same single pass as axpyMulti, but the accumulator +// starts from y, so a caller whose z would otherwise have to be primed with a +// copy of y does not make that copy. +template +void axpyMultiFrom(Lattice &z,const Lattice &y, + const std::vector &b,const std::vector*> &x) +{ + axpyMultiNormImpl(z,b,x,0,&y); +} template RealD axpyMultiNorm(Lattice &z,const std::vector &b,const std::vector*> &x) { diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index 3de398e88..f6a6f2eb9 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -800,6 +800,39 @@ void localCopyRegion(const Lattice &From,Lattice & To,Coordinate Fro autoView(from_v,From,AcceleratorRead); autoView(to_v,To,AcceleratorWrite); + if constexpr ( vobj::Nsimd() == 1 ) { + + // Unvectorised: there is one lane, so the lane indices are identically zero + // and the internal word index is free to carry the thread. Threading over + // (site,word) with the word fastest makes consecutive threads read and write + // consecutive elements; one thread per site instead leaves them `words` + // apart, which is what made this uncoalesced. The site count and the word + // count are folded into one range so the large product lands on the only + // block dimension that can hold it. + accelerator_for(ss,nsite*words,1,{ + + int w = ss % words; + uint64_t idx = ss / words; + + Coordinate from_coor, to_coor, base; + Lexicographic::CoorFromIndex(base,idx,RegionSize); + for(int i=0;i &From,Lattice & To,Coordinate Fro const vector_type* from = (const vector_type *)&from_v[from_oidx]; vector_type* to = (vector_type *)&to_v[to_oidx]; - + scalar_type stmp; for(int w=0;w diff --git a/Grid/lattice/PaddedCell.h b/Grid/lattice/PaddedCell.h index 67d78d0ca..ef271dbf4 100644 --- a/Grid/lattice/PaddedCell.h +++ b/Grid/lattice/PaddedCell.h @@ -94,7 +94,7 @@ template inline void ScatterSlice(const deviceVector &buf, // FIXME -- can put internal indices into thread loop auto buf_p = & buf[0]; autoView(lat_v, lat, AcceleratorWrite); - accelerator_for(ss, face_ovol/simd[dim],Nsimd,{ + accelerator_forNB(ss, face_ovol/simd[dim],Nsimd,{ // scalar layout won't coalesce #ifdef GRID_SIMT @@ -180,7 +180,7 @@ template inline void GatherSlice(deviceVector &buf, //for cross platform //For CPU perhaps just run a loop over Nsimd auto buf_p = & buf[0]; - accelerator_for(ss, face_ovol/simd[dim],Nsimd,{ + accelerator_forNB(ss, face_ovol/simd[dim],Nsimd,{ // scalar layout won't coalesce #ifdef GRID_SIMT @@ -303,10 +303,17 @@ public: template inline Lattice Exchange(const Lattice &in, const CshiftImplBase &cshift = CshiftImplDefault()) const { - GridBase *old_grid = in.Grid(); - int dims = old_grid->Nd(); - Lattice tmp = in; - for(int d=0;d_processors; + int dims = in.Grid()->Nd(); + // An undecomposed dimension is not padded, so Expand on it would only copy + // its input onto the same grid: those are skipped. The first decomposed + // dimension expands from `in` itself, so the chain needs no seed copy. + int first=-1; + for(int d=0;d 1 ) { first=d; break; } + if ( first < 0 ) return in; // nothing decomposed: the padded grid IS the input grid + Lattice tmp = Expand(first,in,cshift); + for(int d=first+1;d inline Lattice ExchangePeriodic(const Lattice &in) const { - GridBase *old_grid = in.Grid(); - int dims = old_grid->Nd(); - Lattice tmp = in; - for(int d=0;d_processors; + int dims = in.Grid()->Nd(); + // As Exchange: undecomposed dimensions are pure copies and are skipped, and + // the first decomposed dimension expands from `in` so there is no seed copy. + // On the D+1 mrhs coarse grid (nrhs,1,x,y,z,t) this removes three full + // coarse-field copies per operator application. + int first=-1; + for(int d=0;d 1 ) { first=d; break; } + if ( first < 0 ) return in; // nothing decomposed: the padded grid IS the input grid + Lattice tmp = ExpandPeriodic(first,in); + for(int d=first+1;d for time from the plaquette form -//(Lüscher: https://arxiv.org/pdf/1006.4518 eq. 3.1) +//(Luscher: https://arxiv.org/pdf/1006.4518 eq. 3.1) //E(t) = 2 * sum_p Retr{ 1 - Vt(p) } = // = 2 * sum_p ( Nc - Retr Vt(p) ) = // = 2 * Nc * sum_p ( 1 - Retr Vt(p)/Nc ) diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h index 86b71e16c..7e7421510 100644 --- a/Grid/simd/Grid_vector_types.h +++ b/Grid/simd/Grid_vector_types.h @@ -1002,10 +1002,10 @@ accelerator_inline void precisionChange(vRealD *out,const vRealF *in,int n Optimization::PrecisionChange::StoD(in[m].v,out[n].v,out[n+1].v); // Bug in gcc 10.0.1 and gcc 10.1 using fixed-size SVE ACLE data types CAS-159553-Y1K4C6 // function call results in compile-time error: - // In function ‘void Grid::precisionChange(Grid::vRealD*, Grid::vRealF*, int)’: + // In function 'void Grid::precisionChange(Grid::vRealD*, Grid::vRealF*, int)': // .../Grid_vector_types.h:961:56: error: - // cannot bind non-const lvalue reference of type ‘vecd&’ {aka ‘svfloat64_t&’} - // to an rvalue of type ‘vecd’ {aka ‘svfloat64_t’} + // cannot bind non-const lvalue reference of type 'vecd&' {aka 'svfloat64_t&'} + // to an rvalue of type 'vecd' {aka 'svfloat64_t'} // 961 | Optimization::PrecisionChange::StoD(in[m].v,out[n].v,out[n+1].v); // | ~~~~~~~^ } diff --git a/MPI_benchmark/halo_mpi.cc b/MPI_benchmark/halo_mpi.cc index 433e3fcf5..b50fa5704 100644 --- a/MPI_benchmark/halo_mpi.cc +++ b/MPI_benchmark/halo_mpi.cc @@ -104,6 +104,119 @@ inline double usecond(void) { gettimeofday(&tv,NULL); return 1.0e6*tv.tv_sec + 1.0*tv.tv_usec; } +/************************************************************** + * Point to point cost, per Cartesian direction. + * + * PaddedCell::Face_exchange issues one MPI_Sendrecv per dimension and waits, + * so the quantity that bounds a halo round is the sendrecv time for that + * direction, not a one way ping-pong. Sweeping the size gives the latency + * intercept and the bandwidth slope separately. + * + * Each direction is labelled on- or off-node, which is what decides whether a + * halo round rides the intra-node fabric or the network, and therefore which + * dimension ordering costs least. Times are reduced across ranks: the max is + * the one that bounds a collective halo round. + ************************************************************** + */ +void PingPong(std::vector cart_geom,bool use_device,int ncall) +{ + int Nd=cart_geom.size(); + std::vector periodic(Nd,1); + std::vector coor(Nd); + int rank; + + MPI_Comm communicator; + MPI_Cart_create(WorldComm,Nd,&cart_geom[0],&periodic[0],0,&communicator); + MPI_Comm_rank(communicator,&rank); + MPI_Cart_coords(communicator,rank,Nd,&coor[0]); + + // An identifier every rank on a node agrees on, so a neighbour can be + // classified by comparing it. + int node_id = WorldRank; + MPI_Bcast(&node_id,1,MPI_INT,0,WorldShmComm); + + size_t max_bytes = 2*1024*1024; + void *xmit, *recv; + if ( use_device ) { + xmit = acceleratorAllocDevice(max_bytes); + recv = acceleratorAllocDevice(max_bytes); + } else { + xmit = malloc(max_bytes); + recv = malloc(max_bytes); + } + + if ( !WorldRank ) { + printf("= dim dir bytes us(max) us(min) MB/s off-node ranks\n"); + fflush(stdout); + } + + for(int d=0;d 10 ? ncall/20 : 10); + + // Untimed warm-up: the first exchange on a fresh buffer pays connection + // setup and registration, and would otherwise land in the first row. + MPI_Sendrecv(xmit,bytes,MPI_CHAR,to,rank, + recv,bytes,MPI_CHAR,from,from, + communicator,MPI_STATUS_IGNORE); + + MPI_Barrier(communicator); + double t0=usecond(); + for(int i=0;i0?"+":"-", bytes, us_max, us_min, + (double)bytes/us_max, offnode_count); + fflush(stdout); + } + } + } + } + + if ( use_device ) { + acceleratorFreeDevice(xmit); + acceleratorFreeDevice(recv); + } else { + free(xmit); + free(recv); + } + MPI_Comm_free(&communicator); +} + /************************************************************** * Main benchmark routine ************************************************************** @@ -303,6 +416,20 @@ int main(int argc, char **argv) } + if( !WorldRank ) { + printf("=========================================================\n"); + printf("= Point to point cost per direction, HOST memory \n"); + printf("=========================================================\n");fflush(stdout); + } + PingPong(mpi,false,1000); + + if( !WorldRank ) { + printf("=========================================================\n"); + printf("= Point to point cost per direction, DEVICE memory \n"); + printf("=========================================================\n");fflush(stdout); + } + PingPong(mpi,true,1000); + if( !WorldRank ) { printf("=========================================================\n"); printf("= Benchmarking HOST memory MPI performance \n"); diff --git a/examples/Example_mdagm.cc b/examples/Example_mdagm.cc index 3b31b5c24..4b2c89304 100644 --- a/examples/Example_mdagm.cc +++ b/examples/Example_mdagm.cc @@ -87,7 +87,7 @@ int main (int argc, char ** argv) HermOpAdaptor HermFineOp(MdagMOp); // ── Coarse geometry ──────────────────────────────────────────────────── - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; @@ -171,7 +171,7 @@ int main (int argc, char ** argv) LatticeFermionD src(FGrid); random(RNG5, src); LatticeFermionD result(FGrid); result = Zero(); - TwoLevelADEF2 + DeprecatedTwoLevelADEF2 HDCG(1.0e-8, 1000, HermFineOp, Smoother, diff --git a/examples/Example_mdagm_cg.cc b/examples/Example_mdagm_cg.cc index 1a67a9ffb..0fe213405 100644 --- a/examples/Example_mdagm_cg.cc +++ b/examples/Example_mdagm_cg.cc @@ -31,7 +31,7 @@ Author: Peter Boyle // // Operator hierarchy: // Fine: M†M, acted on via HermOpAdaptor so Op = HermOp = M†M -// Coarse: 33-point GeneralCoarsenedMatrix (NextToNearest, 2-hop M†M) +// Coarse: 33-point DeprecatedGeneralCoarsenedMatrix (NextToNearest, 2-hop M†M) // // Setup: // 1. Chebyshev filter (T_600 x T_2500) on M†M to build near-null subspace {ψᵢ} @@ -198,7 +198,7 @@ int main(int argc, char **argv) /////////////////////////////////////////////////////////// // Coarse operator: 33-point stencil for 2-hop M†M /////////////////////////////////////////////////////////// - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNearestStencilGeometry5D geom(Coarse5d); diff --git a/examples/Example_pvdagm.cc b/examples/Example_pvdagm.cc index 24b2709f4..149b45f4c 100644 --- a/examples/Example_pvdagm.cc +++ b/examples/Example_pvdagm.cc @@ -364,7 +364,7 @@ void runMG( ) { // typedef Aggregation Subspace; - // typedef GeneralCoarsenedMatrix LittleDiracOperator; + // typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; // typedef LittleDiracOperator::CoarseVector CoarseVector; ParseEnvironment(); @@ -620,7 +620,7 @@ int main (int argc, char ** argv) // assert(nbasis <= Nevecs); // need to have enough evecs - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNearestStencilGeometry5D geom(Coarse5d); diff --git a/examples/Example_pvdagm_3level.cc b/examples/Example_pvdagm_3level.cc index f7c4387b6..27d028be3 100644 --- a/examples/Example_pvdagm_3level.cc +++ b/examples/Example_pvdagm_3level.cc @@ -190,8 +190,8 @@ public: } }; -// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. -// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src. template class LuscherGuesser : public LinearFunction { const std::vector ψ @@ -342,7 +342,7 @@ void runMG( TrivialPrecon simple_fine; ////////////////////////////////////////////////////////////////////// - // Level 0→1: coarsen PVdagM, build LinOpCoarse + // Level 0->1: coarsen PVdagM, build LinOpCoarse ////////////////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); @@ -366,7 +366,7 @@ void runMG( ////////////////////////////////////////////////////////////////////// // psi_coarse: coarse projections of pre-GS fine null vectors. // These are the Level 1 near-null vectors, promoted from Level 0. - // Used as the aggregation basis for Level 1→2 coarsening. + // Used as the aggregation basis for Level 1->2 coarsening. ////////////////////////////////////////////////////////////////////// std::vector psi_coarse(nbasis, Coarse5d); for (int k = 0; k < nbasis; k++) @@ -396,22 +396,22 @@ void runMG( RealD normC = C.norm(); RealD normCmCdag = (C - C.adjoint()).norm(); std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; - std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C^dag||/||C|| = " << normCmCdag/normC << std::endl; std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; } ////////////////////////////////////////////////////////////////////// - // Level 1→2: set up aggregation using psi_coarse as subspace. + // Level 1->2: set up aggregation using psi_coarse as subspace. // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI - // geometry 3×6×4×4 where z-local at Level 1 is 6). + // geometry 3x6x4x4 where z-local at Level 1 is 6). // psi_coarse are assigned directly; CoarsenOperator performs // block-GS orthogonalisation before building LinOpCoarseCoarse. ////////////////////////////////////////////////////////////////////// // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex - // for the L1→L2 level must be iScalar, not vTComplex. + // for the L1->L2 level must be iScalar, not vTComplex. typedef typename CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; typedef MGPreconditioner L1to2MG; @@ -429,12 +429,12 @@ void runMG( TrivialPrecon simpleCC; ////////////////////////////////////////////////////////////////////// - // Lüscher deflation guesser for L3PGCR. + // Luscher deflation guesser for L3PGCR. // Step 1: project psi_coarse[k] (promoted fine null vectors) to - // CoarseCoarseVector space — these cover the zero-momentum + // CoarseCoarseVector space -- these cover the zero-momentum // component of the near-null space of LinOpCC. // Step 2: breed Nextra additional null vectors directly on LinOpCC - // using GCR with random sources — these pick up near-null + // using GCR with random sources -- these pick up near-null // modes at all spatial frequencies not spanned by step 1. // Step 3: build C_{st} = over the // full augmented basis and invert directly via Eigen LU. @@ -476,7 +476,7 @@ void runMG( RealD normCcc = Ccc.norm(); RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; - std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; } Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); @@ -489,7 +489,7 @@ void runMG( L3PGCR.Name("CCouter"); ////////////////////////////////////////////////////////////////////// - // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Coarse-level GCR smoother for Level 1->2 V-cycle. // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. ////////////////////////////////////////////////////////////////////// @@ -505,7 +505,7 @@ void runMG( CoarseSmootherGCR.Name("Csmoother"); ////////////////////////////////////////////////////////////////////// - // Level 1→2 V-cycle preconditioner. + // Level 1->2 V-cycle preconditioner. ////////////////////////////////////////////////////////////////////// L1to2MG L1to2Precon(AggregatesL2, LinOpCoarse, @@ -513,7 +513,7 @@ void runMG( CoarseSmootherGCR, // post-smoother: 12 GCR steps LinOpCC, L3PGCR, - CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag ////////////////////////////////////////////////////////////////////// // Standalone Level 1 two-level solve test. @@ -548,7 +548,7 @@ void runMG( f_src = one; // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. - // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle). TwoLevelMG ThreeLevelPrecon(AggregatesPD, PVdagM, simple_fine, @@ -592,7 +592,7 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + // Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1) Coordinate clatt = lat_size; for (int d = 0; d < 4; d++) clatt[d] /= 2; std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; @@ -600,13 +600,13 @@ int main (int argc, char ** argv) GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); - // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, Ls=1). // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. - // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) -> must use 3. // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). - // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ - // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + // t-block=2 gives t2-local=6, 6 mod 4 != 0, fails Grid SIMD assertion. + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). Coordinate clatt2 = clatt; clatt2[0] /= 2; clatt2[1] /= 2; @@ -635,7 +635,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; typedef MGPreconditioner TwoLevelMG; diff --git a/examples/Example_pvdagm_3level_dense.cc b/examples/Example_pvdagm_3level_dense.cc index 84db5730f..c85ec32a7 100644 --- a/examples/Example_pvdagm_3level_dense.cc +++ b/examples/Example_pvdagm_3level_dense.cc @@ -212,8 +212,8 @@ public: } }; -// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. -// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src. template class LuscherGuesser : public LinearFunction { const std::vector ψ @@ -704,7 +704,7 @@ void runMG( TrivialPrecon simple_fine; ////////////////////////////////////////////////////////////////////// - // Level 0→1: coarsen PVdagM, build LinOpCoarse + // Level 0->1: coarsen PVdagM, build LinOpCoarse ////////////////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); @@ -736,7 +736,7 @@ void runMG( // Optional sigma-ordering of psi_coarse (SVD_REORDER set): replace the crude // first-NB_CC slice with the NB_CC most-null directions of span(psi_coarse) // under LinOpCoarse. Nullness measure = singular values (eig of the - // Gram-whitened Psi†A†APsi), NOT the numerical range Q†AQ which + // Gram-whitened Psi^dag A^dag APsi), NOT the numerical range Q^dag AQ which // non-normality contaminates. The printed sigma spectrum shows where the // truncation cliff sits. Unset => raw first-30 (crude GS-ordered slice). ////////////////////////////////////////////////////////////////////// @@ -787,7 +787,7 @@ void runMG( } ////////////////////////////////////////////////////////////////////// - // Level 1→2: SUPERCOARSE aggregation using psi_coarse as subspace. + // Level 1->2: SUPERCOARSE aggregation using psi_coarse as subspace. // Maximal block {8,4,3,6}: CC = [3,6,8,8] = the dense-invertible floor. // // UNBLOCKING TRUNCATION NB_CC = 30 (first-30 slice of psi_coarse): @@ -811,7 +811,7 @@ void runMG( typedef typename CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; typedef MGPreconditioner L1to2MG; @@ -830,7 +830,7 @@ void runMG( ////////////////////////////////////////////////////////////////////// // CC solve: DENSE (exact, non-iterative) by default; DENSE_CC=0 gives - // the previous iterative L3PGCR + Lüscher-guesser path for A/B. + // the previous iterative L3PGCR + Luscher-guesser path for A/B. ////////////////////////////////////////////////////////////////////// int use_dense = 1; if (getenv("DENSE_CC")) use_dense = atoi(getenv("DENSE_CC")); @@ -859,7 +859,7 @@ void runMG( ccGuess = &simpleCC; // exact solve ignores/overwrites any guess } else { //////////////////////////////////////////////////////////////////// - // Lüscher deflation guesser (arXiv:0706.2298 A.3) for the iterative CC + // Luscher deflation guesser (arXiv:0706.2298 A.3) for the iterative CC // solve, as in the earlier supercoarse configuration. //////////////////////////////////////////////////////////////////// psi_cc.resize(nbasis, CoarseCoarse5d); @@ -900,7 +900,7 @@ void runMG( } ////////////////////////////////////////////////////////////////////// - // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Coarse-level GCR smoother for Level 1->2 V-cycle. ////////////////////////////////////////////////////////////////////// RealD coarse_smoother_shift = 0.1; int coarse_smoother_nstep = 2; @@ -914,7 +914,7 @@ void runMG( CoarseSmootherGCR.SetZeroGuess(1); // post-smoother slot: caller zeroes vec2 (NOT L2MGsolver: it takes the Luscher guess) ////////////////////////////////////////////////////////////////////// - // Level 1→2 V-cycle preconditioner. + // Level 1->2 V-cycle preconditioner. ////////////////////////////////////////////////////////////////////// L1to2MG L1to2Precon(AggregatesL2, LinOpCoarse, @@ -998,7 +998,7 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + // Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1) Coordinate clatt = lat_size; // Coordinate Block1({2,2,2,2}); // Coordinate Block2({8,4,3,6}); @@ -1047,7 +1047,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; typedef MGPreconditioner TwoLevelMG; diff --git a/examples/Example_pvdagm_3level_madj.cc b/examples/Example_pvdagm_3level_madj.cc index 53ad40b9e..e0c4d34d8 100644 --- a/examples/Example_pvdagm_3level_madj.cc +++ b/examples/Example_pvdagm_3level_madj.cc @@ -323,7 +323,7 @@ void runMG( TrivialPrecon simple_fine; ////////////////////////////////////////////////////////////////////// - // Level 0→1: coarsen PVdagM, build LinOpCoarse + // Level 0->1: coarsen PVdagM, build LinOpCoarse ////////////////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); @@ -348,7 +348,7 @@ void runMG( ////////////////////////////////////////////////////////////////////// // psi_coarse: coarse projections of pre-GS fine null vectors. // These are the Level 1 near-null vectors, promoted from Level 0. - // Used as the aggregation basis for Level 1→2 coarsening. + // Used as the aggregation basis for Level 1->2 coarsening. ////////////////////////////////////////////////////////////////////// std::vector psi_coarse(nbasis, Coarse5d); for (int k = 0; k < nbasis; k++) @@ -378,22 +378,22 @@ void runMG( RealD normC = C.norm(); RealD normCmCdag = (C - C.adjoint()).norm(); std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; - std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C^dag||/||C|| = " << normCmCdag/normC << std::endl; std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; } ////////////////////////////////////////////////////////////////////// - // Level 1→2: set up aggregation using psi_coarse as subspace. + // Level 1->2: set up aggregation using psi_coarse as subspace. // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI - // geometry 3×6×4×4 where z-local at Level 1 is 6). + // geometry 3x6x4x4 where z-local at Level 1 is 6). // psi_coarse are assigned directly; CoarsenOperator performs // block-GS orthogonalisation before building LinOpCoarseCoarse. ////////////////////////////////////////////////////////////////////// // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex - // for the L1→L2 level must be iScalar, not vTComplex. + // for the L1->L2 level must be iScalar, not vTComplex. typedef typename CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; typedef MGPreconditioner L1to2MG; @@ -412,14 +412,14 @@ void runMG( // Level 2 solver: plain GCR, no further coarsening ////////////////////////////////////////////////////////////////////// TrivialPrecon simpleCC; - // L3PGCR is an inner solver inside the L1→2 V-cycle; does not need to converge + // L3PGCR is an inner solver inside the L1->2 V-cycle; does not need to converge // to fine-grid precision. Loose tolerance (3e-2) and large restart (64) to allow // the Krylov space to span enough of the near-null spectrum of LinOpCC per cycle. PrecGeneralisedConjugateResidualNonHermitian L3PGCR(1.0e-4,5,LinOpCC,simpleCC,64,64); L3PGCR.Level(3); ////////////////////////////////////////////////////////////////////// - // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Coarse-level GCR smoother for Level 1->2 V-cycle. // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. ////////////////////////////////////////////////////////////////////// @@ -433,7 +433,7 @@ void runMG( CoarseSmootherGCR.Level(2); ////////////////////////////////////////////////////////////////////// - // Level 1→2 V-cycle preconditioner. + // Level 1->2 V-cycle preconditioner. ////////////////////////////////////////////////////////////////////// L1to2MG L1to2Precon(AggregatesL2, LinOpCoarse, @@ -474,7 +474,7 @@ void runMG( f_src = one; // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. - // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle). TwoLevelMG ThreeLevelPrecon(AggregatesPD, PVdagM, simple_fine, @@ -518,7 +518,7 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + // Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1) Coordinate clatt = lat_size; for (int d = 0; d < 4; d++) clatt[d] /= 2; std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; @@ -526,13 +526,13 @@ int main (int argc, char ** argv) GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); - // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, Ls=1). // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. - // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) -> must use 3. // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). - // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ - // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + // t-block=2 gives t2-local=6, 6 mod 4 != 0, fails Grid SIMD assertion. + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). Coordinate clatt2 = clatt; clatt2[0] /= 2; clatt2[1] /= 2; @@ -566,7 +566,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; typedef MGPreconditioner TwoLevelMG; diff --git a/examples/Example_pvdagm_4level.cc b/examples/Example_pvdagm_4level.cc index 2d8d67964..7fe8ed97c 100644 --- a/examples/Example_pvdagm_4level.cc +++ b/examples/Example_pvdagm_4level.cc @@ -190,8 +190,8 @@ public: } }; -// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. -// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src. template class LuscherGuesser : public LinearFunction { const std::vector ψ @@ -343,7 +343,7 @@ void runMG( TrivialPrecon simple_fine; ////////////////////////////////////////////////////////////////////// - // Level 0→1: coarsen PVdagM, build LinOpCoarse + // Level 0->1: coarsen PVdagM, build LinOpCoarse ////////////////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); @@ -367,7 +367,7 @@ void runMG( ////////////////////////////////////////////////////////////////////// // psi_coarse: coarse projections of pre-GS fine null vectors. // These are the Level 1 near-null vectors, promoted from Level 0. - // Used as the aggregation basis for Level 1→2 coarsening. + // Used as the aggregation basis for Level 1->2 coarsening. ////////////////////////////////////////////////////////////////////// std::vector psi_coarse(nbasis, Coarse5d); for (int k = 0; k < nbasis; k++) @@ -397,22 +397,22 @@ void runMG( RealD normC = C.norm(); RealD normCmCdag = (C - C.adjoint()).norm(); std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; - std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C^dag||/||C|| = " << normCmCdag/normC << std::endl; std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; } ////////////////////////////////////////////////////////////////////// - // Level 1→2: set up aggregation using psi_coarse as subspace. + // Level 1->2: set up aggregation using psi_coarse as subspace. // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI - // geometry 3×6×4×4 where z-local at Level 1 is 6). + // geometry 3x6x4x4 where z-local at Level 1 is 6). // psi_coarse are assigned directly; CoarsenOperator performs // block-GS orthogonalisation before building LinOpCoarseCoarse. ////////////////////////////////////////////////////////////////////// // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex - // for the L1→L2 level must be iScalar, not vTComplex. + // for the L1->L2 level must be iScalar, not vTComplex. typedef typename CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; typedef MGPreconditioner L1to2MG; @@ -430,12 +430,12 @@ void runMG( TrivialPrecon simpleCC; ////////////////////////////////////////////////////////////////////// - // Lüscher deflation guesser for L3PGCR. + // Luscher deflation guesser for L3PGCR. // Step 1: project psi_coarse[k] (promoted fine null vectors) to - // CoarseCoarseVector space — these cover the zero-momentum + // CoarseCoarseVector space -- these cover the zero-momentum // component of the near-null space of LinOpCC. // Step 2: breed Nextra additional null vectors directly on LinOpCC - // using GCR with random sources — these pick up near-null + // using GCR with random sources -- these pick up near-null // modes at all spatial frequencies not spanned by step 1. // Step 3: build C_{st} = over the // full augmented basis and invert directly via Eigen LU. @@ -477,13 +477,13 @@ void runMG( RealD normCcc = Ccc.norm(); RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; - std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; } Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); ////////////////////////////////////////////////////////////////////// - // Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation + // Level 2->3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation // to build the Level 4 (coarse-coarse-coarse) operator. // psi_cc[0..nbasis-1] are the coarse-coarse near-null vectors, projected // from the RAW psi_coarse (themselves projected from the RAW fine null @@ -492,11 +492,11 @@ void runMG( // so assign COPIES of psi_cc and keep psi_cc itself raw. // // Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns - // iScalar, so CComplex for the L2→L3 level is iScalar>. + // iScalar, so CComplex for the L2->L3 level is iScalar>. ////////////////////////////////////////////////////////////////////// typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj; typedef iScalar vTTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL3; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL3; typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector; typedef Aggregation SubspaceL3; typedef MGPreconditioner L2to3MG; @@ -531,7 +531,7 @@ void runMG( L4PGCR.Name("CCCouter"); ////////////////////////////////////////////////////////////////////// - // Level 2→3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + L4 bottom. + // Level 2->3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + L4 bottom. // The shift slides the coarse-coarse field of values off the origin so a // 2-step smoother has something to bite on a non-normal operator (IRS idea). ////////////////////////////////////////////////////////////////////// @@ -556,7 +556,7 @@ void runMG( simpleCCC); // trivial guesser at the bottom ////////////////////////////////////////////////////////////////////// - // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2→L3 V-cycle. + // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2->L3 V-cycle. // Replaces the plain L3PGCR of the 3-level build -- the coarse-coarse level // is now smoothed shallowly and recursed rather than solved deeply. ////////////////////////////////////////////////////////////////////// @@ -565,7 +565,7 @@ void runMG( L3MGsolver.Name("CCouter"); ////////////////////////////////////////////////////////////////////// - // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Coarse-level GCR smoother for Level 1->2 V-cycle. // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. ////////////////////////////////////////////////////////////////////// @@ -581,15 +581,15 @@ void runMG( CoarseSmootherGCR.Name("Csmoother"); ////////////////////////////////////////////////////////////////////// - // Level 1→2 V-cycle preconditioner. + // Level 1->2 V-cycle preconditioner. ////////////////////////////////////////////////////////////////////// L1to2MG L1to2Precon(AggregatesL2, LinOpCoarse, simpleC, // no pre-smoother (matches fine-grid setup) CoarseSmootherGCR, // post-smoother: depth-2 shifted GCR LinOpCC, - L3MGsolver, // coarse-coarse solve is now the L2→L3 V-cycle - CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + L3MGsolver, // coarse-coarse solve is now the L2->L3 V-cycle + CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag ////////////////////////////////////////////////////////////////////// // Standalone Level 1 two-level solve test. @@ -624,7 +624,7 @@ void runMG( f_src = one; // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. - // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle). TwoLevelMG ThreeLevelPrecon(AggregatesPD, PVdagM, simple_fine, @@ -668,7 +668,7 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + // Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1) Coordinate clatt = lat_size; for (int d = 0; d < 4; d++) clatt[d] /= 2; std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; @@ -676,13 +676,13 @@ int main (int argc, char ** argv) GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); - // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, Ls=1). // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. - // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) -> must use 3. // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). - // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ - // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + // t-block=2 gives t2-local=6, 6 mod 4 != 0, fails Grid SIMD assertion. + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). Coordinate clatt2 = clatt; clatt2[0] /= 2; clatt2[1] /= 2; @@ -728,7 +728,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; typedef MGPreconditioner TwoLevelMG; diff --git a/examples/Example_pvdagm_5level.cc b/examples/Example_pvdagm_5level.cc index 115421cf1..fc801c696 100644 --- a/examples/Example_pvdagm_5level.cc +++ b/examples/Example_pvdagm_5level.cc @@ -190,8 +190,8 @@ public: } }; -// Lüscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. -// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi† src. +// Luscher deflated guesser (arXiv:0706.2298 Sec A.3) for a non-Hermitian solve. +// C_{st} = ; guess = sum_s c_s psi[s] where c = C^{-1} psi^dag src. template class LuscherGuesser : public LinearFunction { const std::vector ψ @@ -344,7 +344,7 @@ void runMG( TrivialPrecon simple_fine; ////////////////////////////////////////////////////////////////////// - // Level 0→1: coarsen PVdagM, build LinOpCoarse + // Level 0->1: coarsen PVdagM, build LinOpCoarse ////////////////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOpPV(geom, FGrid, Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM, AggregatesPD); @@ -368,7 +368,7 @@ void runMG( ////////////////////////////////////////////////////////////////////// // psi_coarse: coarse projections of pre-GS fine null vectors. // These are the Level 1 near-null vectors, promoted from Level 0. - // Used as the aggregation basis for Level 1→2 coarsening. + // Used as the aggregation basis for Level 1->2 coarsening. ////////////////////////////////////////////////////////////////////// std::vector psi_coarse(nbasis, Coarse5d); for (int k = 0; k < nbasis; k++) @@ -398,22 +398,22 @@ void runMG( RealD normC = C.norm(); RealD normCmCdag = (C - C.adjoint()).norm(); std::cout << GridLogMessage << "Coarse null matrix ||C|| = " << normC << std::endl; - std::cout << GridLogMessage << "Coarse null matrix ||C - C†||/||C|| = " << normCmCdag/normC << std::endl; + std::cout << GridLogMessage << "Coarse null matrix ||C - C^dag||/||C|| = " << normCmCdag/normC << std::endl; std::cout << GridLogMessage << "Galerkin check ||C||/||W|| = " << normC/normW << std::endl; } ////////////////////////////////////////////////////////////////////// - // Level 1→2: set up aggregation using psi_coarse as subspace. + // Level 1->2: set up aggregation using psi_coarse as subspace. // Block factor 2,2,3,2 (removes odd local sublattice in z given MPI - // geometry 3×6×4×4 where z-local at Level 1 is 6). + // geometry 3x6x4x4 where z-local at Level 1 is 6). // psi_coarse are assigned directly; CoarsenOperator performs // block-GS orthogonalisation before building LinOpCoarseCoarse. ////////////////////////////////////////////////////////////////////// // innerProduct(CoarseSiteObj, CoarseSiteObj) returns iScalar, so CComplex - // for the L1→L2 level must be iScalar, not vTComplex. + // for the L1->L2 level must be iScalar, not vTComplex. typedef typename CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; typedef typename LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; typedef MGPreconditioner L1to2MG; @@ -431,12 +431,12 @@ void runMG( TrivialPrecon simpleCC; ////////////////////////////////////////////////////////////////////// - // Lüscher deflation guesser for L3PGCR. + // Luscher deflation guesser for L3PGCR. // Step 1: project psi_coarse[k] (promoted fine null vectors) to - // CoarseCoarseVector space — these cover the zero-momentum + // CoarseCoarseVector space -- these cover the zero-momentum // component of the near-null space of LinOpCC. // Step 2: breed Nextra additional null vectors directly on LinOpCC - // using GCR with random sources — these pick up near-null + // using GCR with random sources -- these pick up near-null // modes at all spatial frequencies not spanned by step 1. // Step 3: build C_{st} = over the // full augmented basis and invert directly via Eigen LU. @@ -478,13 +478,13 @@ void runMG( RealD normCcc = Ccc.norm(); RealD normCccmCdag = (Ccc - Ccc.adjoint()).norm(); std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc|| = " << normCcc << std::endl; - std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc†||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; + std::cout << GridLogMessage << "Coarse-coarse deflation matrix ||Ccc-Ccc^dag||/||Ccc|| = " << normCccmCdag/normCcc << std::endl; } Eigen::MatrixXcd Ccc_inv = Ccc.inverse(); LuscherGuesser CCDeflGuesser(psi_cc, Ccc_inv); ////////////////////////////////////////////////////////////////////// - // Level 2→3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation + // Level 2->3: coarsen LinOpCC using the RAW promoted psi_cc as aggregation // to build the Level 4 (coarse-coarse-coarse) operator. // psi_cc[0..nbasis-1] are the coarse-coarse near-null vectors, projected // from the RAW psi_coarse (themselves projected from the RAW fine null @@ -493,11 +493,11 @@ void runMG( // so assign COPIES of psi_cc and keep psi_cc itself raw. // // Tensor depth deepens once more: innerProduct(CoarseCoarseSiteObj,...) returns - // iScalar, so CComplex for the L2→L3 level is iScalar>. + // iScalar, so CComplex for the L2->L3 level is iScalar>. ////////////////////////////////////////////////////////////////////// typedef typename CoarseCoarseVector::vector_object CoarseCoarseSiteObj; typedef iScalar vTTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL3; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL3; typedef typename LittleDiracOperatorL3::CoarseVector CoarseCoarseCoarseVector; typedef Aggregation SubspaceL3; typedef MGPreconditioner L2to3MG; @@ -514,7 +514,7 @@ void runMG( TrivialPrecon simpleCCC; ////////////////////////////////////////////////////////////////////// - // Level 3→4: coarsen LinOpCCC to build the Level 5 operator, using a + // Level 3->4: coarsen LinOpCCC to build the Level 5 operator, using a // TRUNCATED basis of only the first NB5 (< nbasis) raw promoted null vectors. // psi_ccc[k] = raw psi_cc projected through the (block-GS'd) L3 aggregation // -- the pre-block-GS chain continued one level deeper. We keep only the @@ -526,7 +526,7 @@ void runMG( // NB: a positive result is conservative (sigma-ordering can only help); a // negative one is inconclusive until the sigma-ordered NB5 is tried. // - // Tensor depth deepens once more: CComplex for the L3→L4 level is + // Tensor depth deepens once more: CComplex for the L3->L4 level is // iScalar. NB5 (the coarse dimension) is independent of the // depth -- it just makes the coarsest site vector NB5-dimensional. ////////////////////////////////////////////////////////////////////// @@ -542,10 +542,10 @@ void runMG( // Optional sigma-ordering of psi_ccc (SVD_REORDER set): replace the crude // first-NB5 slice with the NB5 genuinely-most-null directions of span(psi_ccc) // under LinOpCCC. For a NON-NORMAL operator the nullness measure is the - // singular value of A restricted to the span -- eig of Q†A†AQ -- NOT the - // numerical range Q†AQ (which non-normality contaminates). Robust route: + // singular value of A restricted to the span -- eig of Q^dag A^dag AQ -- NOT the + // numerical range Q^dag AQ (which non-normality contaminates). Robust route: // whiten by the Gram (drop near-dependent directions), Hermitian-eig the - // whitened A†A, rotate. The printed singular spectrum IS the SVD study: where + // whitened A^dag A, rotate. The printed singular spectrum IS the SVD study: where // it falls off tells you the natural NB5, and the same numbers illuminate why // the earlier singular-subspace deflation re-entered. Safe here because we // ORDER vectors that then feed a Galerkin projection, not REMOVE a subspace. @@ -603,7 +603,7 @@ void runMG( typedef typename CoarseCoarseCoarseVector::vector_object CoarseCoarseCoarseSiteObj; typedef iScalar vTTTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL4; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL4; typedef typename LittleDiracOperatorL4::CoarseVector CoarseCoarseCoarseCoarseVector; typedef Aggregation SubspaceL4; typedef MGPreconditioner L3to4MG; @@ -635,7 +635,7 @@ void runMG( L5PGCR.Name("CCCCouter"); ////////////////////////////////////////////////////////////////////// - // Level 3→4 V-cycle: depth-2 SHIFTED smoother on LinOpCCC + Level 5 bottom. + // Level 3->4 V-cycle: depth-2 SHIFTED smoother on LinOpCCC + Level 5 bottom. // Level 4 is no longer the bottom -- it is smoothed shallowly and recursed to // Level 5, mirroring how Level 3 recurses to Level 4. ////////////////////////////////////////////////////////////////////// @@ -660,14 +660,14 @@ void runMG( simpleCCCC); // trivial guesser at the bottom ////////////////////////////////////////////////////////////////////// - // Level 4 (coarse-coarse-coarse) solve: GCR preconditioned by the L3→L4 V-cycle. + // Level 4 (coarse-coarse-coarse) solve: GCR preconditioned by the L3->L4 V-cycle. ////////////////////////////////////////////////////////////////////// PrecGeneralisedConjugateResidualNonHermitian L4MGsolver(1.0e-1,200,LinOpCCC,L3to4Precon,16,16); L4MGsolver.Level(4); L4MGsolver.Name("CCCouter"); ////////////////////////////////////////////////////////////////////// - // Level 2→3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + Level 4 solve. + // Level 2->3 V-cycle: depth-2 SHIFTED smoother on LinOpCC + Level 4 solve. // The shift slides the coarse-coarse field of values off the origin so a // 2-step smoother has something to bite on a non-normal operator (IRS idea). ////////////////////////////////////////////////////////////////////// @@ -688,11 +688,11 @@ void runMG( simpleCC, // no pre-smoother CoarseCoarseSmootherGCR, // post-smoother: depth-2 shifted GCR LinOpCCC, - L4MGsolver, // coarse solve is now the L3→L4 V-cycle + L4MGsolver, // coarse solve is now the L3->L4 V-cycle simpleCCC); // trivial guesser ////////////////////////////////////////////////////////////////////// - // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2→L3 V-cycle. + // Level 3 (coarse-coarse) solve: GCR preconditioned by the L2->L3 V-cycle. // Replaces the plain L3PGCR of the 3-level build -- the coarse-coarse level // is now smoothed shallowly and recursed rather than solved deeply. ////////////////////////////////////////////////////////////////////// @@ -701,7 +701,7 @@ void runMG( L3MGsolver.Name("CCouter"); ////////////////////////////////////////////////////////////////////// - // Coarse-level GCR smoother for Level 1→2 V-cycle. + // Coarse-level GCR smoother for Level 1->2 V-cycle. // Mirrors fine-grid SmootherGCR: shifted operator + fixed step count. // coarse_smoother_shift and coarse_smoother_nstep are the tuning knobs. ////////////////////////////////////////////////////////////////////// @@ -717,15 +717,15 @@ void runMG( CoarseSmootherGCR.Name("Csmoother"); ////////////////////////////////////////////////////////////////////// - // Level 1→2 V-cycle preconditioner. + // Level 1->2 V-cycle preconditioner. ////////////////////////////////////////////////////////////////////// L1to2MG L1to2Precon(AggregatesL2, LinOpCoarse, simpleC, // no pre-smoother (matches fine-grid setup) CoarseSmootherGCR, // post-smoother: depth-2 shifted GCR LinOpCC, - L3MGsolver, // coarse-coarse solve is now the L2→L3 V-cycle - CCDeflGuesser); // Lüscher guesser: psi_cc C^{-1} psi_cc† + L3MGsolver, // coarse-coarse solve is now the L2->L3 V-cycle + CCDeflGuesser); // Luscher guesser: psi_cc C^{-1} psi_cc^dag ////////////////////////////////////////////////////////////////////// // Standalone Level 1 two-level solve test. @@ -760,7 +760,7 @@ void runMG( f_src = one; // Pre-smoother: none (TrivialPrecon); post-smoother: shifted PGCR. - // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1→2 V-cycle). + // Coarse solver: L2MGsolver (PGCR preconditioned by Level 1->2 V-cycle). TwoLevelMG ThreeLevelPrecon(AggregatesPD, PVdagM, simple_fine, @@ -804,7 +804,7 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Level 1 coarse grid: block 2^4 from fine (48×48×48×96 → 24×24×24×48, Ls=1) + // Level 1 coarse grid: block 2^4 from fine (48x48x48x96 -> 24x24x24x48, Ls=1) Coordinate clatt = lat_size; for (int d = 0; d < 4; d++) clatt[d] /= 2; std::cout << GridLogMessage << "Level 1 coarse lattice: " << clatt << std::endl; @@ -812,13 +812,13 @@ int main (int argc, char ** argv) GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); - // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24×24×24×48 → 12×12×8×16, Ls=1). + // Level 2 coarse-coarse grid: block 2,2,3,3 from Level 1 (24x24x24x48 -> 12x12x8x16, Ls=1). // MPI geometry 3.6.4.4 (288 ranks): fine local {16,8,12,24}. // Level 1 local {8,4,6,12}; Level 2 local {4,2,2,4}. - // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) → must use 3. + // z blocked by 3: z-Level1-local=6; 6/3=2 (even), 6/2=3 (odd) -> must use 3. // t blocked by 3: t-Level1-local=12; 12/3=4 divisible by Nsimd=4 (gen-simd-width=64). - // t-block=2 gives t2-local=6, 6 mod 4 ≠ 0, fails Grid SIMD assertion. ✓ - // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). ✓ + // t-block=2 gives t2-local=6, 6 mod 4 != 0, fails Grid SIMD assertion. + // With {4,2,2,4}: Nsimd=4 goes into x or t (both =4). Coordinate clatt2 = clatt; clatt2[0] /= 2; clatt2[1] /= 2; @@ -880,7 +880,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; typedef MGPreconditioner TwoLevelMG; diff --git a/examples/Example_pvdagm_census.cc b/examples/Example_pvdagm_census.cc index e04b81ed7..d5c3cd46c 100644 --- a/examples/Example_pvdagm_census.cc +++ b/examples/Example_pvdagm_census.cc @@ -49,7 +49,7 @@ Author: Peter Boyle // sigma_min << min|lambda| : non-normal near origin // lambda_min(H) < 0 : half-plane condition violated // -// Requires the dagger code path in GeneralCoarsenedMatrix: +// Requires the dagger code path in DeprecatedGeneralCoarsenedMatrix: // _Adag allocated, PopulateAdag active, _Adag exchanged, hermitian=0. // // Env vars: @@ -308,7 +308,7 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; diff --git a/examples/Example_pvdagm_coherence.cc b/examples/Example_pvdagm_coherence.cc index e76815b17..67e912c2a 100644 --- a/examples/Example_pvdagm_coherence.cc +++ b/examples/Example_pvdagm_coherence.cc @@ -275,7 +275,7 @@ int main (int argc, char ** argv) MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c); typedef PVdagMLinearOperator PVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; // Grid index contraction is positional (colour/spin/lorentz order is meaningful), // so each MG projection adds one index to the tensor nest rather than reusing a slot: diff --git a/examples/Example_pvdagm_mrhs.cc b/examples/Example_pvdagm_mrhs.cc index c0d14fbcb..351698439 100644 --- a/examples/Example_pvdagm_mrhs.cc +++ b/examples/Example_pvdagm_mrhs.cc @@ -38,7 +38,7 @@ Author: Peter Boyle // precon: V-cycle -- per-RHS fine post-smoother (16-step shifted GCR), // batched restriction (MultiRHSBlockProject / GEMM), // ONE coarse PGCR on the 6D mrhs coarse operator -// (MultiGeneralCoarsenedMatrix, GEMM mults -- the ~10x win), +// (DeprecatedMultiGeneralCoarsenedMatrix, GEMM mults -- the ~10x win), // batched prolongation. // // The coarse operator is coarsened once with the standard single-RHS @@ -260,8 +260,8 @@ int main (int argc, char ** argv) typedef PVdagMLinearOperator PVdagM_t; typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; diff --git a/examples/Example_pvdagm_mrhs_3level.cc b/examples/Example_pvdagm_mrhs_3level.cc index 04704a259..7bd0e8aa5 100644 --- a/examples/Example_pvdagm_mrhs_3level.cc +++ b/examples/Example_pvdagm_mrhs_3level.cc @@ -23,7 +23,7 @@ Author: Peter Boyle // with L3_DEFL=0 (NO deflation), applied to the enlarged block-diagonal mRHS system: // the coarse and coarse-coarse levels run a SINGLE Krylov (one GCR polynomial, inner // products summed over rhs) on the packed 6D mrhs fields, so both coarse levels batch -// through GEMM (MultiGeneralCoarsenedMatrix) -- the valence throughput win at BOTH levels. +// through GEMM (DeprecatedMultiGeneralCoarsenedMatrix) -- the valence throughput win at BOTH levels. // // Level structure (each coarse level is a single-field PGCR on a packed 6D mrhs field): // L1 (fine) : std::vector, MrhsPGCRNonHermitian on PVdagM, @@ -227,16 +227,16 @@ int main (int argc, char ** argv) typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; // Level 1 tensor types - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; // Level 2 tensor types (coarsening deepens the nest by one iScalar -- see CLAUDE.md) typedef CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; diff --git a/examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc index b85a5f17b..1b9415187 100644 --- a/examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc +++ b/examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc @@ -233,16 +233,16 @@ int main (int argc, char ** argv) typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; // Level 1 tensor types - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; // Level 2 tensor types (coarsening deepens the nest by one iScalar) typedef CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; diff --git a/examples/Example_pvdagm_mrhs_3level_dense.cc b/examples/Example_pvdagm_mrhs_3level_dense.cc index 320e63644..2eb7db0af 100644 --- a/examples/Example_pvdagm_mrhs_3level_dense.cc +++ b/examples/Example_pvdagm_mrhs_3level_dense.cc @@ -237,7 +237,7 @@ public: #endif // NB: the apply GEMM (Y = slab^dag X) is a tiny-output/huge-K shape that // under-fills the GPU (~13ms). The fix is a software split-K via - // GridBLAS.gemmBatched (see MultiRHSBlockCGLinalg.h / 2409.03904 Fig 11) — + // GridBLAS.gemmBatched (see MultiRHSBlockCGLinalg.h / 2409.03904 Fig 11) -- // NOT a raw rocblas strided-batched batch, which hung on Frontier and was // removed. TODO: reimplement through GridBLAS when the ~1.4 s/RHS is wanted. @@ -891,16 +891,16 @@ int main (int argc, char ** argv) typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; // Level 1 tensor types - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; typedef Aggregation Subspace; // Level 2 tensor types (coarsening deepens the nest by one iScalar) typedef CoarseVector::vector_object CoarseSiteObj; typedef iScalar vTTComplex; - typedef GeneralCoarsenedMatrix LittleDiracOperatorL2; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorL2; + typedef DeprecatedMultiGeneralCoarsenedMatrix MrhsLittleDiracOperatorL2; typedef LittleDiracOperatorL2::CoarseVector CoarseCoarseVector; typedef Aggregation SubspaceL2; diff --git a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc deleted file mode 100644 index 49e706cc4..000000000 --- a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc +++ /dev/null @@ -1,984 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.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. - - 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 */ - -// -// PVdagM three level multigrid on the V2 coarse operator. -// -// STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings. -// The dense bottom and the solves are not here yet. -// -// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc: -// -// * The coarse space is UNVECTORISED (sComplexD). The fine space stays -// vectorised. MultiRHSBlockProject carries the mixed layout. -// -// * One operator, not two. V1 needed GeneralCoarsenedMatrix to coarsen and -// MultiGeneralCoarsenedMatrix to apply, bridged by CopyMatrix. V2 does -// both, and single versus multiRHS is SetGrid on the same object with the -// matrix elements built once. -// -// * Nrhs is unconstrained. V1 required nrhs % vComplex::Nsimd() == 0 because -// its multiRHS grid carried the SIMD in the rhs direction. -// -// * CoarsenOperator takes the subspace vectors, not an Aggregation. It block -// orthonormalises them IN PLACE -- the vectors are far too large to copy -// defensively -- so rawNull is taken first and the RAW vectors are what -// define the L2 null space. Do not insert an Orthogonalise() anywhere: -// projecting a block-orthonormal vector onto its own block-orthonormalised -// aggregation gives e_k, and the near null content is silently gone. The -// || - I||_F guard below is what catches that. -// -// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH -// HOT_START CONFIG SUBSPACE_FILE V1_CHECK MRHS_COARSEN -// - -#include -#include // std::sort, for the runtime-environment dump in ParseEnvironment -#include -#include -#include -#include -#include - -#include - -using namespace std; -using namespace Grid; - -// Compile time so it can be cut down for laptop runs: -DNBASIS=8 -#ifndef NBASIS -#define NBASIS 60 -#endif - -RealD mass = 0.00078; -int Nrhs = 12; -int Ls = 24; -int CoarsenBatch = 9; -std::vector lat_size({48,48,48,96}); - -// Solver tuning. PRINCIPLE (PB, 2026-08-24): the defaults ARE the current -// optimum, so an unset environment reproduces the best banked result; they -// are updated as and when a better point is found, and every change is -// dated here. Environment variables of the same names override for sweeps. -// -// Current optimum: 2026-08-24, slurm-5335492 F4, 48^3x96 Ls=24 on 288 GCDs, -// 17.2 s/RHS at Nrhs=4, 32.2 s at Nrhs=1 (exact-halo FINAL ~1e-8 pending -// the exact-outer rerun). Smoother mmax == order (full GCR history); -// PB's mmax=1 trial gave 72 vs ~60 outer iterations and was slower. -RealD FineSmootherShift = 0.1; -int FineSmootherOrder = 6; -int FineSmootherMmax = 6; -RealD CoarseSmootherShift = 0.1; -int PowerIterations = 0; // >0: power-iterate the smoother operators before the solves (spectral edge) -// Smoother implementation per level (Smoothers.h): -// gcr : the adaptive PGCR (default, as always) -// replay : run the PGCR with a coefficient recorder for the first -// PolyRecordIters outer steps, then switch to GCRReplaySmoother -// (same polynomial, no inner products) -- 1402.2585 p.13 revisited -// cheb : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order -// = the GCR step count of that level -std::string FineSmootherMode = "gcr"; -std::string CoarseSmootherMode = "gcr"; -int PolyRecordIters = 4; -int PolyRecordStart = 0; // outer step at which recording begins (0: from the first step) -std::string PolyRecordSelect = "last"; // which recorded call to replay: last|first|mean (mean is the bad one) -int PolyRefresh = 0; // >0: every PolyRefresh outer steps, one adaptive step re-records the polynomial (HDCG: every 10) -int PolyVerbose = 0; // 1: fixed-polynomial smoothers print |r_m|/|r_0| per call -RealD FineChebLo = 3.0, FineChebHi = 137.0; // from the harvested polynomial and PowerIterations edge -RealD CoarseChebLo = 8.0, CoarseChebHi = 45.0; -int CoarseSmootherNstep = 2; -int CoarseSmootherMmax = 2; -RealD CoarseSolverTol = 0.05; -int CoarseSolverOrder = 200; -int CoarseSolverMmax = 16; -RealD OuterTol = 1.0e-8; -int OuterMmax = 4; -int OuterNstep = 8; - -// "It's legal to get the same answer faster, not to get a less correct -// answer." (PB, 2026-08-24) -// -// Halo-precision POLICY: reduced-precision (fp32 wire) -// halos belong in the PRECONDITIONER -- the smoother, the V-cycle's own -// residuals, and the coarsening -- and NEVER in the outer Krylov. The -// outer operator's applications define what "converged" means; making them -// sloppy turns the stopping criterion into a statement about the wrong -// operator (measured: solver stops at computed 9.8e-9 while the true -// residual is 3.4e-8). Exactness costs one exact fine matvec per outer -// iteration, a few percent of the solve. -// -// Stencil::SloppyComms is a free runtime setter (Stencil.h:303), so the -// policy is implemented by SCOPED toggling: SetFineSloppy(1) on entering -// the preconditioner / coarsening, SetFineSloppy(0) on leaving. The -// operators default to EXACT. FineSloppyComms therefore now means -// "sloppy inside the preconditioner"; =0 makes everything exact. -int FineSloppyComms = 1; -// SmootherCoeffLog=1 : print the GCR step lengths a_k and orthogonalisation -// coefficients b_kj of BOTH smoothers every call -- the harvest for a fixed -// polynomial smoother (stable coefficients => stationary p(A), no reductions). -int SmootherCoeffLog = 0; -std::function SetFineSloppy = [](int){}; - -void ParseEnvironment(void) -{ - if(getenv("MASS")) mass = atof(getenv("MASS")); - if(getenv("NRHS")) Nrhs = atoi(getenv("NRHS")); - if(getenv("LS")) Ls = atoi(getenv("LS")); - if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH")); - if(getenv("FineSmootherShift")) FineSmootherShift = atof(getenv("FineSmootherShift")); - if(getenv("FineSmootherOrder")) FineSmootherOrder = atoi(getenv("FineSmootherOrder")); - if(getenv("FineSmootherMmax")) FineSmootherMmax = atoi(getenv("FineSmootherMmax")); - if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift")); - if(getenv("PowerIterations")) PowerIterations = atoi(getenv("PowerIterations")); - if(getenv("FineSmootherMode")) FineSmootherMode = getenv("FineSmootherMode"); - if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode"); - if(getenv("PolyRecordIters")) PolyRecordIters = atoi(getenv("PolyRecordIters")); - if(getenv("PolyRecordStart")) PolyRecordStart = atoi(getenv("PolyRecordStart")); - if(getenv("PolyRecordSelect")) PolyRecordSelect = getenv("PolyRecordSelect"); - if(getenv("PolyRefresh")) PolyRefresh = atoi(getenv("PolyRefresh")); - if(getenv("PolyVerbose")) PolyVerbose = atoi(getenv("PolyVerbose")); - if(getenv("FineChebLo")) FineChebLo = atof(getenv("FineChebLo")); - if(getenv("FineChebHi")) FineChebHi = atof(getenv("FineChebHi")); - if(getenv("CoarseChebLo")) CoarseChebLo = atof(getenv("CoarseChebLo")); - if(getenv("CoarseChebHi")) CoarseChebHi = atof(getenv("CoarseChebHi")); - if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep")); - if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax")); - if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol")); - if(getenv("CoarseSolverOrder")) CoarseSolverOrder = atoi(getenv("CoarseSolverOrder")); - if(getenv("CoarseSolverMmax")) CoarseSolverMmax = atoi(getenv("CoarseSolverMmax")); - if(getenv("OuterTol")) OuterTol = atof(getenv("OuterTol")); - if(getenv("OuterMmax")) OuterMmax = atoi(getenv("OuterMmax")); - if(getenv("FineSloppyComms")) FineSloppyComms = atoi(getenv("FineSloppyComms")); - if(getenv("SmootherCoeffLog")) SmootherCoeffLog = atoi(getenv("SmootherCoeffLog")); - if(getenv("OuterNstep")) OuterNstep = atoi(getenv("OuterNstep")); - if(getenv("LATT")){ - Coordinate l; - GridCmdOptionIntVector(std::string(getenv("LATT")),l); - GRID_ASSERT(l.size()==4); - for(int d=0;d<4;d++) lat_size[d]=l[d]; - } - - std::cout << GridLogMessage << "PARAM: LATT " - << lat_size[0]<<"."< hits; - for(char **e = environ; e && *e; e++){ - std::string s(*e); - for(auto p : prefixes){ - if( s.compare(0,strlen(p),p)==0 ){ - if(s.size()>200) s = s.substr(0,197)+"..."; // paths can be enormous - hits.push_back(s); - break; - } - } - } - std::sort(hits.begin(),hits.end()); - std::cout << GridLogMessage << "PARAM: ---- runtime environment: "< -void saveSubspace(std::vector &subspace, std::string const fname){ -#ifdef HAVE_LIME - Grid::emptyUserRecord record; - Grid::ScidacWriter SW(subspace[0].Grid()->IsBoss()); - SW.open(fname); - for (int k = 0; k < (int)subspace.size(); k++) SW.writeScidacFieldRecord(subspace[k], record); - SW.close(); -#endif -} -template -void loadSubspace(std::vector &subspace, std::string const fname){ -#ifdef HAVE_LIME - Grid::emptyUserRecord record; - Grid::ScidacReader SR; - SR.open(fname); - for (int k = 0; k < (int)subspace.size(); k++) SR.readScidacFieldRecord(subspace[k], record); - SR.close(); -#endif -} - -////////////////////////////////////////////////////////////////////// -// A = PV^dag M (non-Hermitian) -////////////////////////////////////////////////////////////////////// -template -class PVdagMLinearOperator : public LinearOperatorBase { - Matrix &_Mat; Matrix &_PV; -public: - PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV) {}; - void OpDiag (const Field &in, Field &out) { assert(0); } - void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } - void OpDirAll (const Field &in, std::vector &out){ assert(0); }; - void Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); } - void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(in,tmp); _Mat.Mdag(tmp,out); } - void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ HermOp(in,out); ComplexD d=innerProduct(in,out); n1=real(d); n2=norm2(out); } - void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } -}; - -////////////////////////////////////////////////////////////////////// -// || - I||_F over a set of coarse vectors. Small means the raw near null -// content survived the projection; see GramGuard for where a leak lands. -////////////////////////////////////////////////////////////////////// -template -RealD GramDefect(std::vector &v) -{ - RealD s2=0.0; - for(int i=0;i<(int)v.size();i++){ - for(int j=0;j<(int)v.size();j++){ - ComplexD sij=TensorRemove(innerProduct(v[i],v[j])); - ComplexD d=sij-(i==j?ComplexD(1.0):ComplexD(0.0)); - s2+=real(d)*real(d)+imag(d)*imag(d); - } - } - return std::sqrt(s2); -} - -// On a leak every image collapses to the block unit e_k, the Gram becomes -// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2 -// of a content preserving projection. Trip well below that so a mis-set -// threshold costs a log line rather than the run. -template -void GramGuard(const std::string &name,std::vector &v,GridBase *grid) -{ - RealD defect = GramDefect(v); - RealD N = (RealD)grid->gSites(); - RealD leak = (N-1.0)*std::sqrt((RealD)v.size()); - RealD trip = std::sqrt(N); - std::cout << GridLogMessage << "GUARD: ||<"< - I||_F = " << defect - << " (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl; - GRID_ASSERT( defect < trip ); -} - - -////////////////////////////////////////////////////////////////////// -// Shifted variants for the smoothers -////////////////////////////////////////////////////////////////////// -template -class ShiftedPVdagMLinearOperator : public LinearOperatorBase { - Matrix &_Mat; Matrix &_PV; -public: - RealD shift; - ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){}; - void OpDiag (const Field &in, Field &out) { assert(0); } - void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } - void OpDirAll (const Field &in, std::vector &out){ assert(0); }; - void Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; } - void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; } - void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); } - void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } -}; - -template -class ShiftedLinearOperator : public LinearOperatorBase { - LinearOperatorBase &_Op; RealD shift; -public: - ShiftedLinearOperator(RealD _shift, LinearOperatorBase &Op) : _Op(Op), shift(_shift) {} - void OpDiag (const Field &in, Field &out) { assert(0); } - void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } - void OpDirAll (const Field &in, std::vector &out) { 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){ assert(0); } - void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } -}; - -////////////////////////////////////////////////////////////////////// -// Power iteration on a (non-Hermitian) operator: the spectral edge the -// smoother polynomial must not exceed. Reports -// step 0 : |A v|/|v| on a RANDOM unit v -- a one-sample lower bound on -// sigma_max(A). If this and the converged value agree, the -// operator is near-normal and the spectral picture (R_m(lambda) -// on the spectrum) is trustworthy; if not, the field of values -// sets the safe interval and the spectrum understates it. -// step k : |A v_k|/|v_k| -> |lambda_max| as v_k -> the dominant -// eigenvector; the complex Rayleigh quotient gives its -// phase (real => on the axis). A non-converging oscillation -// means a complex-conjugate pair of equal modulus at the top. -// Uses Op(), not HermOp(): this is the operator the smoother sees. -////////////////////////////////////////////////////////////////////// -template -void PowerIteration(const std::string &name, LinearOperatorBase &Op, GridBase *grid, int iters) -{ - GRID_TRACE("PowerIteration"); - GridParallelRNG RNG(grid); RNG.SeedFixedIntegers(std::vector({7,11,13,17})); - Field v(grid), Av(grid); - gaussian(RNG,v); - RealD nv = std::sqrt(norm2(v)); v = v*(1.0/nv); - RealD ratio=0.0, ratio0=0.0; ComplexD rq(0.0); - for(int i=0;iL2 blocking; banked optimum 2026-08-24 (env BLOCK overrides) - if ( getenv("BLOCK") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK")),Block); GRID_ASSERT(Block.size()==4); } - for(int d=0;d<4;d++){ GRID_ASSERT(lat_size[d]%Block[d]==0); clatt[d]=lat_size[d]/Block[d]; } - std::cout << GridLogMessage << "Block " << Block << " coarse lattice " << clatt << std::endl; - - ////////////////////////////////////////////////////////////////////// - // The coarse space is unvectorised. The 5D coarse grid is built here - // rather than through SpaceTimeGrid so the SIMD layout is ours. - ////////////////////////////////////////////////////////////////////// - Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); - Coordinate c5simd({1,1,1,1,1}); - Coordinate c5mpi ({1,mpi[0],mpi[1],mpi[2],mpi[3]}); - GridCartesian *Coarse5d = new GridCartesian(c5latt,c5simd,c5mpi); - - // 6D coarse multiRHS grid: rhs is dim 0, undistributed and unvectorised. - // No divisibility constraint on nrhs, unlike V1. - Coordinate cmlatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]}); - Coordinate cmsimd({1,1,1,1,1,1}); - Coordinate cmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); - GridCartesian *CoarseMrhs = new GridCartesian(cmlatt,cmsimd,cmmpi); - - // 6D coarse grid at the coarsening batch, used only while CoarsenOperator - // runs. The matrix elements survive the change back to nrhs. - Coordinate cblatt({batch,1,clatt[0],clatt[1],clatt[2],clatt[3]}); - GridCartesian *CoarseBatch = new GridCartesian(cblatt,cmsimd,cmmpi); - - // 6D fine grid carrying the coarsening batch: fine SIMD layout preserved - Coordinate fmlatt({batch,Ls,lat_size[0],lat_size[1],lat_size[2],lat_size[3]}); - Coordinate fmsimd({1,1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]}); - Coordinate fmmpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]}); - GridCartesian *FineMrhs = new GridCartesian(fmlatt,fmsimd,fmmpi); - - std::cout << GridLogMessage << "Nsimd fine " << FGrid->Nsimd() - << " coarse " << Coarse5d->Nsimd() << std::endl; - - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers({1,2,3,4}); - GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers({5,6,7,8}); - - ////////////////////////////////////////////////////////////////////// - // Gauge field - ////////////////////////////////////////////////////////////////////// - LatticeGaugeField Umu(UGrid); - if ( getenv("HOT_START") ) { - std::cout << GridLogMessage << "Hot start gauge field" << std::endl; - SU::HotConfiguration(RNG4,Umu); - } else { - std::string file("/ccs/home/poare/ckpoint_lat.1000"); - if ( getenv("CONFIG") ) file = std::string(getenv("CONFIG")); - std::cout << GridLogMessage << "Reading gauge field " << file << std::endl; - FieldMetaData header; - NerscIO::readConfiguration(Umu,header,file); - } - - MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); - MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c); - - // PVdagM and ShiftedPVdagM are thin wrappers over these same two objects, - // so this one callback controls every fine halo in the program. Default - // EXACT; the preconditioner and the coarsening turn sloppiness on for - // their own scope only (policy note at FineSloppyComms). - SetFineSloppy = [&Ddwf,&Dpv](int sloppy){ - Ddwf.SloppyComms(sloppy); - Dpv .SloppyComms(sloppy); - }; - SetFineSloppy(0); - std::cout << GridLogMessage << "Fine halo policy: preconditioner+coarsening " - << (FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact") - << ", outer Krylov EXACT" << std::endl; - - typedef PVdagMLinearOperator PVdagM_t; - typedef ShiftedPVdagMLinearOperator ShiftedPVdagM_t; - PVdagM_t PVdagM(Ddwf,Dpv); - - ////////////////////////////////////////////////////////////////////// - // Level 1 types: unvectorised coarse scalar - ////////////////////////////////////////////////////////////////////// - typedef sTComplexD CComplexS; - typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; - typedef CoarseOperator::CoarseVector CoarseVector; - typedef Aggregation Subspace; - - NextToNearestStencilGeometry5D geom(Coarse5d); - - ////////////////////////////////////////////////////////////////////// - // Subspace: load RAW (no Orthogonalise!), or generate. - // - // The Aggregation is scaffolding for CreateSubspaceGCR only. That runs - // entirely on the fine grid and ends in GlobalOrthonormalise, which is a - // whole-lattice Gram-Schmidt, so the coarse grid it holds is never - // dereferenced and may be the unvectorised one. - ////////////////////////////////////////////////////////////////////// - std::string subspace_file = "subspace_nb" + std::to_string(nbasis) + ".scidac"; - if ( getenv("SUBSPACE_FILE") ) subspace_file = std::string(getenv("SUBSPACE_FILE")); - uint64_t file_exists=0; - if ( UGrid->IsBoss() ){ std::ifstream f(subspace_file); file_exists=f.good()?1:0; } - UGrid->GlobalSum(file_exists); - - const int cb=0; - Subspace AggregatesGCR(Coarse5d,FGrid,cb); - if ( file_exists ){ - std::cout << GridLogMessage << "*** Loading subspace from disk (kept RAW) ***" << std::endl; - loadSubspace(AggregatesGCR.subspace, subspace_file); - } else { - std::cout << GridLogMessage << "*** GCR subspace generation ***" << std::endl; - AggregatesGCR.CreateSubspaceGCR(RNG5,PVdagM,nbasis); - saveSubspace(AggregatesGCR.subspace, subspace_file); - } - - // RAW copy BEFORE CoarsenOperator block-orthonormalises in place. - std::vector rawNull(nbasis,FGrid); - for(int k=0;k MrhsPVdagM(PVdagM,FGrid,batch); - CoarseOpPV.CoarsenOperator(MrhsPVdagM,FineMrhs,AggregatesGCR.subspace,Coarse5d); - } else { - // PVdagM is single RHS: apply it directly, batch on the coarse side. - CoarseOpPV.CoarsenOperator(PVdagM,AggregatesGCR.subspace,Coarse5d,batch); - } - SetFineSloppy(0); - - // (moved here 2026-08-28: it needs AggregatesGCR.subspace, which is freed right after - // the projector imports it below -- the 60 fine vectors are 20 GB of host memory per - // rank and 60 LRU-eligible device fields competing with the solver's working set) - ////////////////////////////////////////////////////////////////////// - // Optional cross check of the coarse matrix elements against the V1 - // path, which needs a vectorised coarse space. Block Gram-Schmidt is - // idempotent, so V1 may re-orthonormalise the same vectors in place - // without a second copy of the subspace. - ////////////////////////////////////////////////////////////////////// - if ( getenv("V1_CHECK") ) { - - typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef MultiGeneralCoarsenedMatrix MrhsLittleDiracOperator; - typedef Aggregation SubspaceV; - - Coordinate v5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); - Coordinate v5simd({1,fsimd[0],fsimd[1],fsimd[2],fsimd[3]}); - GridCartesian *Coarse5dV = new GridCartesian(v5latt,v5simd,c5mpi); - - int nrhs_v1 = vComplex::Nsimd(); - Coordinate vmlatt({nrhs_v1,1,clatt[0],clatt[1],clatt[2],clatt[3]}); - Coordinate vmsimd({vComplex::Nsimd(),1,1,1,1,1}); - GridCartesian *CoarseMrhsV = new GridCartesian(vmlatt,vmsimd,cmmpi); - - NextToNearestStencilGeometry5D geomV(Coarse5dV); - - SubspaceV AggV(Coarse5dV,FGrid,cb); - for(int k=0;k h1(sites),h2(sites); - acceleratorCopyFromDevice(&mrhsV1.BLAS_A[p][0], &h1[0],sites*sizeof(calcMatrix)); - acceleratorCopyFromDevice(&CoarseOpPV.BLAS_A[p][0],&h2[0],sites*sizeof(calcMatrix)); - ComplexD *w1=(ComplexD *)&h1[0]; - ComplexD *w2=(ComplexD *)&h2[0]; - int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); - for(int64_t i=0;i L3. - // - // The fine operator here is V2 at L1, which is natively multiRHS, so the - // multiRHS driver applies with no promotion adapter: its D+1 grid IS the - // batch grid the L1 operator is currently set to. - ////////////////////////////////////////////////////////////////////// - Coordinate cclatt = clatt; - Coordinate Block2({4,4,2,4}); // L2->L3 blocking; banked optimum 2026-08-24 (env BLOCK2 overrides) - if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); } - for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; } - std::cout << GridLogMessage << "Block2 " << Block2 << " coarse-coarse lattice " << cclatt << std::endl; - - Coordinate cc5latt({1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); - GridCartesian *CoarseCoarse5d = new GridCartesian(cc5latt,c5simd,c5mpi); - - Coordinate ccmlatt({nrhs,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); - GridCartesian *CoarseCoarseMrhs = new GridCartesian(ccmlatt,cmsimd,cmmpi); - - Coordinate ccblatt({batch,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); - GridCartesian *CoarseCoarseBatch = new GridCartesian(ccblatt,cmsimd,cmmpi); - - // Coarsening deepens the tensor nest by one iScalar - typedef CoarseVector::vector_object CoarseSiteObj; - typedef iScalar CComplexS2; - typedef MultiGeneralCoarsenedOperatorV2 CoarseCoarseOperator; - typedef CoarseCoarseOperator::CoarseVector CoarseCoarseVector; - - NextToNearestStencilGeometry5D geom2(CoarseCoarse5d); - - // RAW copy of the coarse null vectors, for the same reason as rawNull: - // the L2 CoarsenOperator block-orthonormalises its subspace in place, and - // the L3 basis must be defined by the vectors that still carry content. - std::vector rawPsi(nbasis,Coarse5d); - for(int k=0;k LinOpCoarse(CoarseOpPV); - - std::cout << GridLogMessage << "*** L2 CoarsenOperator, batch "< MrhsProjectorL2; - MrhsProjectorL2.Allocate(nbasis,Coarse5d,CoarseCoarse5d); - MrhsProjectorL2.ImportBasis(psi_coarse); // block orthonormal basis - - { - std::vector psi_cc(nbasis,CoarseCoarse5d); - MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in - - GramGuard("psi_cc",psi_cc,CoarseCoarse5d); - } - rawPsi.clear(); rawPsi.shrink_to_fit(); - - ////////////////////////////////////////////////////////////////////// - // Both coarse operators apply on their solve grids - ////////////////////////////////////////////////////////////////////// - { - GridParallelRNG cRNG(Coarse5d); cRNG.SeedFixedIntegers({3,4,5,6}); - CoarseVector cin(CoarseMrhs), cout_(CoarseMrhs); - random(cRNG,cin); - CoarseOpPV.M(cin,cout_); - std::cout << GridLogMessage << "L1 apply |in|^2 = " << norm2(cin) - << " |M in|^2 = " << norm2(cout_) << std::endl; - GRID_ASSERT( norm2(cout_) > 0.0 ); - - GridParallelRNG ccRNG(CoarseCoarse5d); ccRNG.SeedFixedIntegers({7,8,9,10}); - CoarseCoarseVector ccin(CoarseCoarseMrhs), ccout(CoarseCoarseMrhs); - random(ccRNG,ccin); - CoarseOpL2.M(ccin,ccout); - std::cout << GridLogMessage << "L2 apply |in|^2 = " << norm2(ccin) - << " |M in|^2 = " << norm2(ccout) << std::endl; - GRID_ASSERT( norm2(ccout) > 0.0 ); - } - - ////////////////////////////////////////////////////////////////////// - // STAGE THREE (part one): the dense bottom on L2. - // - // DenseCoarseMatrix is bilingual: it takes the elements through - // Geometry()/ExtractMatrix(), so the V2 operator serves directly. It does - // detect that a multiRHS op cannot apply on the D dimensional grid and - // skips its own certificate and VERIFY, so the equivalent check is done - // here instead, driving the L2 operator at Nrhs 1 through a slice. - ////////////////////////////////////////////////////////////////////// - typedef DenseCoarseMatrix DenseCC_t; - std::unique_ptr DenseCC; - - if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) { - - std::cout << GridLogMessage << "*** L3 dense bottom: import from the V2 L2 operator ***" << std::endl; - DenseCC.reset(new DenseCC_t(CoarseCoarse5d)); - DenseCC->Import(CoarseOpL2); - - //////////////////////////////////////////////////////////////////// - // ||A Ainv x - x|| / ||x||, the check Import could not run itself - //////////////////////////////////////////////////////////////////// - Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); - GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi); - - CoarseOpL2.SetGrid(CoarseCoarseOne); - - CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d); - GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({11,12,13,14}); - random(dRNG,x); - - (*DenseCC)(x,y); // y = Ainv x - - CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne); - InsertSliceFast(y,y1,0,0); - CoarseOpL2.M(y1,z1); // z = A y - ExtractSliceFast(z,z1,0,0); - - z = z - x; - RealD rel = std::sqrt(norm2(z)/norm2(x)); - std::cout << GridLogMessage << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl; - GRID_ASSERT( rel < 1.0e-2 ); - - CoarseOpL2.SetGrid(CoarseCoarseMrhs); - delete CoarseCoarseOne; - } - - ////////////////////////////////////////////////////////////////////// - // STAGE THREE (part two): the solves. - // - // Both operators are driven from the SAME objects at whatever Nrhs is - // asked for -- the matrix elements were built once and survive SetGrid -- - // so single RHS and multiRHS are the same code path with a different grid. - ////////////////////////////////////////////////////////////////////// - GRID_ASSERT(DenseCC != nullptr); // the PGCR bottom is not ported yet - - typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; - - ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); - TrivialPrecon simple_fine; - TrivialPrecon simpleC; - - auto RunSolve = [&](int nr) - { - std::cout << GridLogMessage << "**********************************************" << std::endl; - std::cout << GridLogMessage << " V2 THREE-level solve, Nrhs = " << nr << std::endl; - // Device-memory budget BEFORE the solve (2026-08-28: NRHS=6 died in hipMalloc at the - // first fine-smoother history allocation -- reported as an asynchronous "memory - // access fault" unless AMD_SERIALIZE_KERNEL/COPY made it a clean OOM). The outer - // mRHS GCR holds src, sol, r, Az and OuterMmax x (p,q) fine fields PER RHS; the fine - // smoother adds FineSmootherMmax x (p,q) once. The MemoryManager LRU cap - // (--device-mem) must be BELOW what is physically left after the non-LRU allocations - // (comms buffers, dense slab, stencil buffers), or the device fills before anything - // is evicted. Print the estimate, the LRU state, and the device's own free count. - { - uint64_t fieldBytes = (uint64_t)FGrid->lSites()*sizeof(typename LatticeFermionD::scalar_object); - double outerGB = (double)nr*(4 + 2*OuterMmax)*fieldBytes/1.0e9; - double smthGB = (double)(2*FineSmootherMmax + 4)*fieldBytes/1.0e9; - std::cout << GridLogMessage << "Device budget: fine field " << fieldBytes/1.0e6 << " MB; outer GCR history " - << nr << " x (4 + 2 x " << OuterMmax << ") fields = " << outerGB << " GB; fine smoother history + temps ~ " - << smthGB << " GB; MemoryManager device LRU " << MemoryManager::DeviceCacheBytes()/1.0e9 << " GB now, cap " - << MemoryManager::DeviceMaxBytes/1.0e9 << " GB" << std::endl; - // Empty the device LRU: every setup-era Lattice copy (coarse null vectors, - // coarsening temporaries) goes back to host, so the solve's working set - // starts from a clean device and the cap applies to it alone. - MemoryManager::EvictAll(); - // ...and release the allocation caches' held blocks (setup-era deviceVector - // scratch that is "free" to the caller but not to hipMalloc). - MemoryManager::DropCache(); - MemoryManager::PrintBytes(); - acceleratorMem(); - } - std::cout << GridLogMessage << "**********************************************" << std::endl; - - Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]}); - Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]}); - GridCartesian *CMrhs = new GridCartesian(cml, cmsimd,cmmpi); - GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi); - - CoarseOpPV.SetGrid(CMrhs); - CoarseOpL2.SetGrid(CCMrhs); - - NonHermitianLinearOperator LinOpC (CoarseOpPV); - NonHermitianLinearOperator LinOpCC(CoarseOpL2); - - MrhsDenseCCSolve ccSolve(*DenseCC,nr); - - ShiftedLinearOperator ShiftedC(CoarseSmootherShift, LinOpC); - if ( PowerIterations > 0 ) { - // Spectral edges the smoother polynomials must respect (see - // scripts/gcr_polynomial.py: |R_m|>1 beyond the edge = amplification). - PowerIteration ("CoarseSmootherOp(shift="+std::to_string(CoarseSmootherShift)+")", ShiftedC, CMrhs, PowerIterations); - PowerIteration ("CoarseOp(unshifted)", LinOpC, CMrhs, PowerIterations); - PowerIteration("FineSmootherOp(shift="+std::to_string(FineSmootherShift)+")", ShiftedPVdagM, Ddwf.FermionGrid(), PowerIterations); - } - PrecGeneralisedConjugateResidualNonHermitian - CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep); - CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1); - - SwitchableSmoother CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR"); - MrhsCoarseThreeLevelPrec - L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve, - Coarse5d, CoarseCoarse5d, CCMrhs, nr); - - PrecGeneralisedConjugateResidualNonHermitian - L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, CoarseSolverMmax, 16); - L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1); - - FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherMmax,FineSmootherOrder); - SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1); - SmootherGCR.LogCoefficients(SmootherCoeffLog); - CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog); - - SwitchableSmoother FineSmootherSlot(SmootherGCR,"Fsmoother GCR"); - MrhsTwoLevelMG > - ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs); - ThreeLevelPrecon.SetSloppy = SetFineSloppy; - ThreeLevelPrecon.SloppyComms = FineSloppyComms; - - MrhsPGCRNonHermitian - L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep); - L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1); - - ////////////////////////////////////////////////////////////////////// - // Smoother modes. Objects live for the duration of this RunSolve. - ////////////////////////////////////////////////////////////////////// - std::unique_ptr > FineCheb; - std::unique_ptr > CoarseCheb; - std::unique_ptr > FineReplay; - std::unique_ptr > CoarseReplay; - GCRCoefficients recF, recC; - { - GCRCoefficients::Select sel = GCRCoefficients::Last; - if ( PolyRecordSelect=="first" ) sel = GCRCoefficients::First; - if ( PolyRecordSelect=="mean" ) sel = GCRCoefficients::Mean; - recF.select = sel; recC.select = sel; - } - if ( FineSmootherMode == "cheb" ) { - FineCheb.reset(new ChebyshevNonHermitianSmoother(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"); - } - // Recording window [PolyRecordStart, PolyRecordStart+PolyRecordIters). - // M3 (2026-08-26): the GCR polynomial changes fast over the first outer - // steps (per-call |r|/|r0| 0.0034 -> 0.017 over steps 1-4) and the MEAN - // of those is a poor smoother (replay 0.033-0.049 per call); record a - // settled window instead. - if ( PolyRecordStart == 0 ) { - if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF); - if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC); - } - // Record -> replay, with optional periodic re-recording ("re-record, not - // fade away": HDCG refreshed its polynomial every 10 steps, tracking the - // evolving spectral content of the residual). Schedule on outer steps: - // [PolyRecordStart, +PolyRecordIters) adaptive GCR, recording - // then replay of the selected recorded call; - // if PolyRefresh>0: every PolyRefresh steps, ONE adaptive recording - // step, then replay of that call. - auto BuildReplays = [&](void){ - if ( FineSmootherMode == "replay" ) { - SmootherGCR.SetCoefficientRecorder(nullptr); - recF.Flush(); recF.Report("Fsmoother"); - FineReplay.reset(new GCRReplaySmoother(ShiftedPVdagM,recF)); - FineReplay->Verbose = PolyVerbose; FineReplay->name = "Fsmoother"; - FineSmootherSlot.Set(*FineReplay,"Fsmoother replay"); - SmootherGCR.ReleaseHistory(); - } - if ( CoarseSmootherMode == "replay" ) { - CoarseSmootherGCR.SetCoefficientRecorder(nullptr); - recC.Flush(); recC.Report("Csmoother"); - CoarseReplay.reset(new GCRReplaySmoother(ShiftedC,recC)); - CoarseReplay->Verbose = PolyVerbose; CoarseReplay->name = "Csmoother"; - CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay"); - CoarseSmootherGCR.ReleaseHistory(); - } - }; - auto StartRecording = [&](int step){ - if ( FineSmootherMode == "replay" ) { { auto sel=recF.select; recF = GCRCoefficients(); recF.select=sel; } SmootherGCR.SetCoefficientRecorder(&recF); FineSmootherSlot.Set(SmootherGCR,"Fsmoother GCR (recording)"); } - if ( CoarseSmootherMode == "replay" ) { { auto sel=recC.select; recC = GCRCoefficients(); recC.select=sel; } CoarseSmootherGCR.SetCoefficientRecorder(&recC); CoarseSmootherSlot.Set(CoarseSmootherGCR,"Csmoother GCR (recording)"); } - std::cout << GridLogMessage << "Smoother coefficient recording starts at outer step " << step << std::endl; - }; - int switchStep = PolyRecordStart + PolyRecordIters; - L1PGCR.OnStep = [&](int step){ - if ( FineSmootherMode != "replay" && CoarseSmootherMode != "replay" ) return; - if ( step == PolyRecordStart && PolyRecordStart > 0 ) StartRecording(step); - if ( step == switchStep ) { BuildReplays(); return; } - if ( PolyRefresh > 0 && step > switchStep ) { - int since = step - switchStep; - if ( since % PolyRefresh == 0 ) { StartRecording(step); } // one adaptive, recorded step - if ( since % PolyRefresh == 1 ) { BuildReplays(); } // then replay it - } - }; - std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << " coarse " << CoarseSmootherMode - << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? " (record outer steps "+std::to_string(PolyRecordStart)+".."+std::to_string(PolyRecordStart+PolyRecordIters)+")" : "") - << std::endl; - - std::vector src(nr,FGrid), sol(nr,FGrid); - for(int r=0;r $fname 2>&1 echo " exit $?"; sleep 10 - grep -h "V2 3-level solve\|Memory access fault\|GRID_ASSERT\|Converged on iteration" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | tail -4 | cut -c1-140 + grep -h "3-level solve\|Memory access fault\|GRID_ASSERT\|Converged on iteration" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | tail -4 | cut -c1-140 grep -n -B3 "Memory access fault" $fname | tail -5 | cut -c1-160 } @@ -124,5 +124,5 @@ run_cell F3_nrhs6_nodense 16000 NRHS=6 DENSE_CC=0 run_cell F4_nrhs6_bigmem 48000 NRHS=6 # no MemoryManager eviction echo "=========================================================" -for f in log.fault.F*; do echo "$f: $(grep -c 'Memory access fault' $f) faults; $(grep -h 'V2 3-level solve' $f | tail -1 | sed 's/.*V2/V2/' | cut -c1-60)"; done +for f in log.fault.F*; do echo "$f: $(grep -c 'Memory access fault' $f) faults; $(grep -h '3-level solve' $f | tail -1 | sed 's/.*3-level/3-level/' | cut -c1-60)"; done echo "=========================================================" diff --git a/systems/Frontier/nrhs_fault_36.job b/systems/Frontier/nrhs_fault_36.job index 23328c2f2..3083bc78c 100644 --- a/systems/Frontier/nrhs_fault_36.job +++ b/systems/Frontier/nrhs_fault_36.job @@ -67,7 +67,7 @@ export MPICH_SMP_SINGLE_COPY_MODE=CMA export MPICH_OFI_NIC_POLICY=GPU module load libfabric -BIN=$root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix +BIN=$root/examples/Example_pvdagm_3level_DenseCoarseMatrix OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000" OVERLAP="--comms-overlap" # set empty for a no-overlap job vol=48.48.48.96 @@ -112,7 +112,7 @@ run_cell () { tail -15 $(dirname $f)/Grid.stdout.$rank | cut -c1-160 else echo " no 'Memory access fault' in any Grid.stderr.*" - grep -h "V2 3-level solve\|Fouter MrhsPGCR: Converged" $GRID_STDOUT_ROOT/0/Grid.stdout.0 | tail -3 | cut -c1-120 + grep -h "3-level solve\|Fouter MrhsPGCR: Converged" $GRID_STDOUT_ROOT/0/Grid.stdout.0 | tail -3 | cut -c1-120 fi } diff --git a/systems/Frontier/pvdagm_multigrid.job b/systems/Frontier/pvdagm_multigrid.job index 017a7809c..314d18a78 100644 --- a/systems/Frontier/pvdagm_multigrid.job +++ b/systems/Frontier/pvdagm_multigrid.job @@ -70,7 +70,7 @@ export GRID_ALLOC_NCACHE_LARGE=0 # allocator large-ring depth (Grid alloca module load libfabric BIN=$root/examples/Example_pvdagm_multigrid -OPTS="--accelerator-threads 8 --shm 4096 --shm-mpi 0 --device-mem 32000 --comms-overlap" +OPTS="--accelerator-threads 8 --shm 4096 --shm-mpi 0 --device-mem 40000 --comms-overlap" vol=48.48.48.96 MPI_GEOM=3.6.4.4 diff --git a/systems/Frontier/schur2d_ladder.job b/systems/Frontier/schur2d_ladder.job index 891c8184b..9cd40ecde 100644 --- a/systems/Frontier/schur2d_ladder.job +++ b/systems/Frontier/schur2d_ladder.job @@ -131,14 +131,14 @@ export NRHS=4 echo "----- F4a: 1D baseline (DENSE_SCHUR2D unset) -----" unset DENSE_SCHUR2D -srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ +srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_3level_DenseCoarseMatrix \ --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap if [ $F3RC -eq 0 ] then echo "----- F4b: 2D block-cyclic (DENSE_SCHUR2D=1) -----" export DENSE_SCHUR2D=1 -srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ +srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_3level_DenseCoarseMatrix \ --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap else echo "----- F4b SKIPPED: F3 failed (rc=$F3RC), not spending the setup on it -----" diff --git a/systems/Frontier/smoother_modes.job b/systems/Frontier/smoother_modes.job index cb0b7a9d1..aebc7effc 100644 --- a/systems/Frontier/smoother_modes.job +++ b/systems/Frontier/smoother_modes.job @@ -121,10 +121,10 @@ run_cell () { export FineSmootherMode=$5; export CoarseSmootherMode=$6 echo "----- $name : Fso=$FineSmootherOrder Fss=$FineSmootherShift Csn=$CoarseSmootherNstep fine=$FineSmootherMode coarse=$CoarseSmootherMode -----" fname=log.ladder.$name - srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \ + srun -N36 -n288 --kill-on-bad-exit=1 ./select_gpu $root/examples/Example_pvdagm_3level_DenseCoarseMatrix \ --mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1 echo " exit $?"; sleep 60 - echo " $(grep -h 'V2 3-level solve Nrhs' $fname | sed 's/.*V2/V2/' | tr '\n' ' ')" + echo " $(grep -h '3-level solve Nrhs' $fname | sed 's/.*3-level/3-level/' | tr '\n' ' ')" echo " $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | cut -c1-60 | tr '\n' ' ')" echo " replay per-call |r|/|r0| (Nrhs=1 solve): $(awk '/THREE-level solve, Nrhs = 1/{s=1} s && /Fsmoother replay \|r\|/{v=$NF; n++; t+=v; if(v>mx)mx=v} END{if(n) printf "mean %.4f max %.4f over %d calls", t/n, mx, n}' $fname)" grep -h "SCHUR fp64 distributed invert took\|GB/s/rank\|BIG LEAVES\|ring histogram\|^Grid : Message : [0-9.]* s : >=" $fname | sed 's/^Grid : Message : [0-9.]* s : //' | cut -c1-150 | head -12 diff --git a/tests/Test_fft_memory.cc b/tests/Test_fft_memory.cc index 40cd6085f..f7aee04a8 100644 --- a/tests/Test_fft_memory.cc +++ b/tests/Test_fft_memory.cc @@ -8,7 +8,7 @@ then repeatedly applies FFT_all_dim to the same propagator 400 times. If PlannedFFT is working correctly the RSS should remain flat after the first - iteration — no new plans, no new deviceVector allocations beyond the per-call + iteration -- no new plans, no new deviceVector allocations beyond the per-call pencil buffer which is freed at the end of each FFT_dim_execute call. Build exactly like any other Grid test, e.g.: @@ -63,7 +63,7 @@ static long getGPUUsedMb() } // ============================================================ -// Convenience struct — one snapshot of both sides +// Convenience struct -- one snapshot of both sides // ============================================================ struct MemSnapshot { long cpu_rss_kb; // host RSS in kB (-1 if unavailable) @@ -113,7 +113,7 @@ int main(int argc, char **argv) << "Grid is setup to use " << threads << " threads" << std::endl; // ------------------------------------------------------------------ - // Grid setup — use whatever lattice/mpi/simd was passed on the CLI, + // Grid setup -- use whatever lattice/mpi/simd was passed on the CLI, // e.g. --grid 8.8.8.8 --mpi 1.1.1.1 // ------------------------------------------------------------------ Coordinate latt_size = GridDefaultLatt(); @@ -174,13 +174,13 @@ int main(int argc, char **argv) << std::endl; // ------------------------------------------------------------------ - // Create the PlannedFFT — plans are allocated here ONCE for all + // Create the PlannedFFT -- plans are allocated here ONCE for all // dimensions and stored inside the object. // ------------------------------------------------------------------ PlannedFFT> plannedFFT(&GRID); // ------------------------------------------------------------------ - // Snapshot AFTER plan construction — this is the true baseline + // Snapshot AFTER plan construction -- this is the true baseline // for the loop, because cufftPlanMany itself grabs device memory. // ------------------------------------------------------------------ MemSnapshot snap_after_plan = takeSnapshot(); @@ -197,7 +197,7 @@ int main(int argc, char **argv) // ------------------------------------------------------------------ // 400-iteration loop. // Each iteration computes the full 4d forward FFT of `prop`. - // We deliberately do NOT cache the result — we always start from + // We deliberately do NOT cache the result -- we always start from // the same `prop` so the FFT is recomputed identically each time. // The point is to watch memory, not correctness. // ------------------------------------------------------------------ @@ -229,7 +229,7 @@ int main(int argc, char **argv) } // cudaMemGetInfo reflects the state *after* any pooled frees have - // been committed, so this is accurate without an explicit sync — + // been committed, so this is accurate without an explicit sync -- // FFT_dim_execute already calls accelerator_barrier() internally. MemSnapshot snap_now = takeSnapshot(); printRow(i, snap_now, snap_prev); diff --git a/tests/core/Test_planned_fft.cc b/tests/core/Test_planned_fft.cc index 4778aa1e3..e3a0c305f 100644 --- a/tests/core/Test_planned_fft.cc +++ b/tests/core/Test_planned_fft.cc @@ -156,7 +156,7 @@ int main (int argc, char ** argv) } //////////////////////////////////////////////////// - // Dwf matrix — verify Fourier representation using PlannedFFT + // Dwf matrix -- verify Fourier representation using PlannedFFT //////////////////////////////////////////////////// { std::cout<<"****************************************"< - - 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 */ - -// -// MultiGeneralCoarsenedOperatorV2 against the existing mrhs coarse operator. -// -// V1 is constructed on the D+1 grid as now; V2 on the D dimensional grid, -// with SetGrid() adopting the caller owned D+1 grid and building its padded -// cell and neighbour table from the D dimensional stencil, with the Nrhs -// factor multiplied in. -// -// Both are given identical matrix elements, so any difference in the apply is -// the restructured neighbour table. The same geometry object is passed to -// both: V1 adds one to skip for the rhs direction, V2 uses it as is on the D -// dimensional grid, so both describe the same stencil over the D dimensions. -// -#include - -using namespace Grid; - -const int nbasis = 8; - -typedef vSpinColourVector FineObj; -typedef sTComplexD CComplexT; // unvectorised coarse space - -typedef MultiGeneralCoarsenedMatrix MrhsV1; -typedef MultiGeneralCoarsenedOperatorV2 MrhsV2; - -//////////////////////////////////////////////////////////////////////// -// Identical random matrix elements into both operators -//////////////////////////////////////////////////////////////////////// -template -void SeedMatrixElements(OpA &A,OpB &B,int npoint,GridSerialRNG &sRNG) -{ - typedef typename OpA::calcMatrix calcMatrix; - - for(int p=0;p host(sites); - - ComplexD *w = (ComplexD *)&host[0]; - int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); - for(int64_t i=0;i -RealD MatrixChecksum(Op &O,int npoint) -{ - typedef typename Op::calcMatrix calcMatrix; - RealD sum=0.0; - for(int p=0;p host(sites); - acceleratorCopyFromDevice(&O.BLAS_A[p][0],&host[0],sites*sizeof(calcMatrix)); - ComplexD *w = (ComplexD *)&host[0]; - int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); - for(int64_t i=0;iNsimd() << std::endl; - - //////////////////////////////////////////////// - // One geometry object for both, and one D+1 grid - // owned here and shared by both operators: fields - // conform only across a shared grid object. - //////////////////////////////////////////////// - NextToNearestStencilGeometry4D geom(CoarseD); - - Coordinate mlatt(1,nrhs), msimd(1,1), mmpi(1,1); - for(int d=0;dNsimd() << std::endl; - - std::cout << GridLogMessage << "npoint V1 " << OpV1.geom.npoint - << " npoint V2 " << OpV2.geom.npoint << std::endl; - GRID_ASSERT(OpV1.geom.npoint == OpV2.geom.npoint); - - int npoint = OpV1.geom.npoint; - - //////////////////////////////////////////////// - // Identical matrix elements - //////////////////////////////////////////////// - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(std::vector({7,8,9,10})); - SeedMatrixElements(OpV1,OpV2,npoint,sRNG); - - RealD ckV1 = MatrixChecksum(OpV1,npoint); - RealD ckV2 = MatrixChecksum(OpV2,npoint); - std::cout << GridLogMessage << "matrix element checksum V1 " << ckV1 - << " V2 " << ckV2 << std::endl; - GRID_ASSERT( ckV1 == ckV2 ); - - //////////////////////////////////////////////// - // Same input, compare the applies - //////////////////////////////////////////////// - typedef MrhsV1::CoarseVector CoarseVector; - // RNG on the D dimensional grid fills any D+1 field: the rhs direction is - // undistributed and divides cleanly. One RNG serves every Nrhs. - GridParallelRNG pRNG(CoarseD); pRNG.SeedFixedIntegers(std::vector({1,2,3,4})); - - CoarseVector in (CoarseMulti); random(pRNG,in); - CoarseVector out1(CoarseMulti); - CoarseVector out2(CoarseMulti); - CoarseVector err (CoarseMulti); - - OpV1.M(in,out1); - OpV2.M(in,out2); - - err = out1 - out2; - std::cout << GridLogMessage << "|V1 out|^2 = " << norm2(out1) - << " |V2 out|^2 = " << norm2(out2) << std::endl; - std::cout << GridLogMessage << "|V1 - V2|^2 = " << norm2(err) << std::endl; - GRID_ASSERT( norm2(out1) > 0.0 ); - GRID_ASSERT( norm2(err) == 0.0 ); - - //////////////////////////////////////////////// - // SetGrid is idempotent on pointer identity - //////////////////////////////////////////////// - OpV2.SetGrid(CoarseMulti); - GRID_ASSERT( MatrixChecksum(OpV2,npoint) == ckV2 ); - OpV2.M(in,out2); - err = out1 - out2; - GRID_ASSERT( norm2(err) == 0.0 ); - std::cout << GridLogMessage << "SetGrid idempotent on identity" << std::endl; - - //////////////////////////////////////////////// - // Move to a different Nrhs and back. The matrix - // elements are Nrhs independent and must survive - // both the release and the rebuild. - //////////////////////////////////////////////// - // Nrhs 1 is the single RHS case through the multiRHS path, and each slice - // of the Nrhs 4 apply must come back unchanged. - OpV2.M(in,out2); - for(int nr=2;nr>=1;nr--){ - Coordinate latt2(1,nr), simd2(1,1), mpi2(1,1); - for(int d=0;d 2 -> 1 -> release -> 4, |V1 - V2|^2 = " - << norm2(err) << std::endl; - GRID_ASSERT( norm2(err) == 0.0 ); - - std::cout << GridLogMessage << "Test_coarse_v2: ALL PASS" << std::endl; - - Grid_finalize(); -} diff --git a/tests/debug/Test_coarse_v2_coarsen.cc b/tests/debug/Test_coarse_v2_coarsen.cc deleted file mode 100644 index 37f1bdff2..000000000 --- a/tests/debug/Test_coarse_v2_coarsen.cc +++ /dev/null @@ -1,273 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./tests/debug/Test_coarse_v2_coarsen.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. - - 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 */ - -// -// Coarsen the same fine operator two ways and compare the matrix elements. -// -// V1 : existing mrhs CoarsenOperator, vectorised coarse space matched to -// the fine SIMD layout, single RHS fine applications -// V2 : D+1 CoarsenOperator, unvectorised sComplexD coarse space, the batch -// of phased basis vectors carried in the rhs direction and applied -// through MrhsPromotedOperator -// -// BLAS_A is written by GridtoBLAS in lSite order, which does not depend on -// the SIMD layout, so the two are directly comparable. -// -#include - -using namespace Grid; - -const int nbasis = 8; -const int batch = 9; - -typedef vSpinColourVector FineObj; -typedef vTComplex CComplexV; // vectorised coarse space, V1 reference -typedef sTComplexD CComplexS; // unvectorised coarse space, V2 - -typedef MultiGeneralCoarsenedMatrix MrhsV1; -typedef MultiGeneralCoarsenedOperatorV2 MrhsV2; - -template -void ReadMatrix(Op &O,int npoint,std::vector > &host) -{ - typedef typename Op::calcMatrix calcMatrix; - host.resize(npoint); - for(int p=0;pNsimd() - << " coarse V1 "<Nsimd() - << " coarse V2 "<Nsimd()<({1,2,3,4})); - LatticeGaugeFieldD Umu(FineGrid); SU::HotConfiguration(pRNG,Umu); - - RealD mass = 0.1; - WilsonFermionD Dw(Umu,*FineGrid,*FrbGrid,mass); - MdagMLinearOperator HermOp(Dw); - - //////////////////////////////////////////////// - // One random subspace, two copies: CoarsenOperator orthogonalises in place - //////////////////////////////////////////////// - Aggregation Subspace(CoarseV,FineGrid,0); - std::vector subspace(nbasis,FineGrid); - for(int i=0;i MrhsHermOp(HermOp,FineGrid,batch); - - std::cout << GridLogMessage << "V2 CoarsenOperator (D+1)" << std::endl; - OpV2.CoarsenOperator(MrhsHermOp,FineGridMulti,subspace,CoarseS); - - //////////////////////////////////////////////// - // Compare matrix elements - //////////////////////////////////////////////// - int npoint = OpV1.geom.npoint; - GRID_ASSERT(npoint == OpV2.geom.npoint); - typedef MrhsV1::calcMatrix calcMatrix; - std::vector > A1,A2; - ReadMatrix(OpV1,npoint,A1); - ReadMatrix(OpV2,npoint,A2); - - RealD num=0.0, den=0.0; - for(int p=0;p > A3; - ReadMatrix(OpV2s,npoint,A3); - - RealD nums=0.0; - for(int p=0;p LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -223,7 +223,7 @@ int main (int argc, char ** argv) GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); #if 0 - MultiGeneralCoarsenedMatrix mrhs(LittleDiracOp,CoarseMrhs); + DeprecatedMultiGeneralCoarsenedMatrix mrhs(LittleDiracOp,CoarseMrhs); typedef decltype(mrhs) MultiGeneralCoarsenedMatrix_t; ////////////////////////////////////////// diff --git a/tests/debug/Test_general_coarse_hdcg.cc b/tests/debug/Test_general_coarse_hdcg.cc index 73c3b2721..80e30ea87 100644 --- a/tests/debug/Test_general_coarse_hdcg.cc +++ b/tests/debug/Test_general_coarse_hdcg.cc @@ -96,7 +96,7 @@ int main (int argc, char ** argv) //////////////////////////////////////////////////////////// ///////////// Coarse basis and Little Dirac Operator /////// //////////////////////////////////////////////////////////// - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -140,7 +140,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); mrhs.CoarsenOperator(FineHermOp,Aggregates,Coarse5d); diff --git a/tests/debug/Test_general_coarse_hdcg_phys.cc b/tests/debug/Test_general_coarse_hdcg_phys.cc index 6e2145c50..eb854f080 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys.cc @@ -189,7 +189,7 @@ int main (int argc, char ** argv) //////////////////////////////////////////////////////////// ///////////// Coarse basis and Little Dirac Operator /////// //////////////////////////////////////////////////////////// - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -371,7 +371,7 @@ slurm-1482367.out:Grid : Message : 6169.469330 s : HDCG: Pcg converged in 487 it ////////////////////////////////////////// // Build a HDCG solver ////////////////////////////////////////// - TwoLevelADEF2 + DeprecatedTwoLevelADEF2 HDCG(1.0e-8, 700, FineHermOp, CGsmooth, diff --git a/tests/debug/Test_general_coarse_hdcg_phys48.cc b/tests/debug/Test_general_coarse_hdcg_phys48.cc index 910e404f6..4529aac5a 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys48.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys48.cc @@ -218,7 +218,7 @@ int main (int argc, char ** argv) //////////////////////////////////////////////////////////// ///////////// Coarse basis and Little Dirac Operator /////// //////////////////////////////////////////////////////////// - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -269,7 +269,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); std::cout << "**************************************"< LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -456,7 +456,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); std::cout << "**************************************"< LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -455,7 +455,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); std::cout << "**************************************"< LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); diff --git a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc index e57fcc050..2e275d84b 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc @@ -161,7 +161,7 @@ int main (int argc, char ** argv) //////////////////////////////////////////////////////////// ///////////// Coarse basis and Little Dirac Operator /////// //////////////////////////////////////////////////////////// - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -210,7 +210,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); std::cout << "**************************************"< LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); @@ -217,7 +217,7 @@ int main (int argc, char ** argv) Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1}); GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); - typedef MultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; + typedef DeprecatedMultiGeneralCoarsenedMatrix MultiGeneralCoarsenedMatrix_t; MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs); /////////////////////// diff --git a/tests/debug/Test_general_coarse_pvdagm.cc b/tests/debug/Test_general_coarse_pvdagm.cc index d0ea894c1..220f3ec90 100644 --- a/tests/debug/Test_general_coarse_pvdagm.cc +++ b/tests/debug/Test_general_coarse_pvdagm.cc @@ -254,7 +254,7 @@ int main (int argc, char ** argv) const int cb = 0 ; LatticeFermion prom(FGrid); - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNearestStencilGeometry5D geom(Coarse5d); diff --git a/tests/debug/Test_general_coarse_pvdagm_svd.cc b/tests/debug/Test_general_coarse_pvdagm_svd.cc index 85589ebd4..015bac143 100644 --- a/tests/debug/Test_general_coarse_pvdagm_svd.cc +++ b/tests/debug/Test_general_coarse_pvdagm_svd.cc @@ -376,10 +376,10 @@ int main (int argc, char ** argv) std::cout < LittleDiracOperatorV; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorV; typedef LittleDiracOperatorV::CoarseVector CoarseVectorV; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; V.Orthogonalise(); diff --git a/tests/debug/Test_general_coarse_pvdagm_svd_cg.cc b/tests/debug/Test_general_coarse_pvdagm_svd_cg.cc index 06f7632ef..2d5c06ff9 100644 --- a/tests/debug/Test_general_coarse_pvdagm_svd_cg.cc +++ b/tests/debug/Test_general_coarse_pvdagm_svd_cg.cc @@ -375,10 +375,10 @@ int main (int argc, char ** argv) std::cout < LittleDiracOperatorV; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperatorV; typedef LittleDiracOperatorV::CoarseVector CoarseVectorV; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; V.Orthogonalise(); diff --git a/tests/debug/Test_general_coarse_pvdagm_svd_uv.cc b/tests/debug/Test_general_coarse_pvdagm_svd_uv.cc index 299178feb..20b8466be 100644 --- a/tests/debug/Test_general_coarse_pvdagm_svd_uv.cc +++ b/tests/debug/Test_general_coarse_pvdagm_svd_uv.cc @@ -369,8 +369,8 @@ int main (int argc, char ** argv) // } - // typedef GeneralCoarsenedMatrix LittleDiracOperator; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + // typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM,V,V); diff --git a/tests/debug/Test_general_coarse_wilson.cc b/tests/debug/Test_general_coarse_wilson.cc index d2845cb86..d2b294102 100644 --- a/tests/debug/Test_general_coarse_wilson.cc +++ b/tests/debug/Test_general_coarse_wilson.cc @@ -185,7 +185,7 @@ int main (int argc, char ** argv) const int cb = 0 ; LatticeFermion prom(FGrid); - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NearestStencilGeometry4D geom(Coarse4d); diff --git a/tests/debug/Test_general_coarse_wilson_nog5.cc b/tests/debug/Test_general_coarse_wilson_nog5.cc index c3382c25f..c444c84e6 100644 --- a/tests/debug/Test_general_coarse_wilson_nog5.cc +++ b/tests/debug/Test_general_coarse_wilson_nog5.cc @@ -185,7 +185,7 @@ int main (int argc, char ** argv) const int cb = 0 ; LatticeFermion prom(FGrid); - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NearestStencilGeometry4D geom(Coarse4d); diff --git a/tests/debug/Test_general_coarse_wilson_svd.cc b/tests/debug/Test_general_coarse_wilson_svd.cc index 1d0588783..ff0109449 100644 --- a/tests/debug/Test_general_coarse_wilson_svd.cc +++ b/tests/debug/Test_general_coarse_wilson_svd.cc @@ -182,7 +182,7 @@ int main (int argc, char ** argv) const int cb = 0 ; LatticeFermion prom(FGrid); - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NearestStencilGeometry4D geom(Coarse4d); diff --git a/tests/debug/Test_general_coarse_wilson_svd_no5g.cc b/tests/debug/Test_general_coarse_wilson_svd_no5g.cc index ad6e59faf..70d25b4cf 100644 --- a/tests/debug/Test_general_coarse_wilson_svd_no5g.cc +++ b/tests/debug/Test_general_coarse_wilson_svd_no5g.cc @@ -182,7 +182,7 @@ int main (int argc, char ** argv) const int cb = 0 ; LatticeFermion prom(FGrid); - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NearestStencilGeometry4D geom(Coarse4d); diff --git a/tests/debug/Test_hipfft_bug_fail.cc b/tests/debug/Test_hipfft_bug_fail.cc index 990250d16..9cc26e7d2 100644 --- a/tests/debug/Test_hipfft_bug_fail.cc +++ b/tests/debug/Test_hipfft_bug_fail.cc @@ -3,9 +3,9 @@ * * Tests three orderings with an empty rocFFT cache to find which GPU * operation before plan creation triggers the failure: - * A) hipMalloc only — hypothesis: passes (no async GPU work) - * B) hipMalloc + hipMemset — hypothesis: fails (async GPU work in flight) - * C) hipMalloc + hipMemset — hypothesis: passes (work completed before plan) + * A) hipMalloc only -- hypothesis: passes (no async GPU work) + * B) hipMalloc + hipMemset -- hypothesis: fails (async GPU work in flight) + * C) hipMalloc + hipMemset -- hypothesis: passes (work completed before plan) * + hipDeviceSynchronize * * Compile: @@ -59,13 +59,13 @@ int main(void) { hipfftResult rvB = makePlan(G, howmany); printf("G=%-4d B) hipMalloc + hipMemset : %s\n", G, res(rvB)); - // C: hipMalloc + hipMemset + sync — does syncing before plan creation fix it? + // C: hipMalloc + hipMemset + sync -- does syncing before plan creation fix it? hipMemset(buf, 0, nelems * sizeof(hipfftDoubleComplex)); hipDeviceSynchronize(); hipfftResult rvC = makePlan(G, howmany); printf("G=%-4d C) hipMalloc + hipMemset + sync: %s\n", G, res(rvC)); - // A last: hipMalloc only, no async GPU work — should always pass + // A last: hipMalloc only, no async GPU work -- should always pass hipfftResult rvA = makePlan(G, howmany); printf("G=%-4d A) hipMalloc only : %s\n\n", G, res(rvA)); diff --git a/tests/debug/Test_hipfft_bug_pass.cc b/tests/debug/Test_hipfft_bug_pass.cc index 3cd7c2103..43974cde2 100644 --- a/tests/debug/Test_hipfft_bug_pass.cc +++ b/tests/debug/Test_hipfft_bug_pass.cc @@ -33,7 +33,7 @@ int main(void) { int n[] = {G}; long nelems = (long)G * howmany; - // Plan created BEFORE hipMalloc — succeeds for all G + // Plan created BEFORE hipMalloc -- succeeds for all G hipfftHandle p; size_t workSize = 0; hipfftCreate(&p); diff --git a/tests/debug/Test_hipfft_minimal.cc b/tests/debug/Test_hipfft_minimal.cc index bbccff3c0..629793a30 100644 --- a/tests/debug/Test_hipfft_minimal.cc +++ b/tests/debug/Test_hipfft_minimal.cc @@ -38,8 +38,8 @@ static const char *hipfftResultString(hipfftResult r) { // Plan creation + execution for (G, howmany). // Tests two orderings to isolate whether a prior hipMalloc poisons hipfft // plan creation for small G on ROCm 7: -// A) plan BEFORE hipMalloc — hypothesis: succeeds -// B) hipMalloc BEFORE plan — hypothesis: fails for G < 32 +// A) plan BEFORE hipMalloc -- hypothesis: succeeds +// B) hipMalloc BEFORE plan -- hypothesis: fails for G < 32 static void tryPlanAndExec(int G, long howmany) { int n[] = {G}; long nelems = (long)G * howmany; diff --git a/tests/debug/Test_hipfft_repro.cc b/tests/debug/Test_hipfft_repro.cc index 12f0d3d61..b6c1d1a37 100644 --- a/tests/debug/Test_hipfft_repro.cc +++ b/tests/debug/Test_hipfft_repro.cc @@ -58,12 +58,12 @@ static void tryPlanAndExec(int G, long howmany) { hipfftDoubleComplex *dbuf = nullptr; hipError_t herr = hipMalloc(&dbuf, nelems * sizeof(hipfftDoubleComplex)); if (herr != hipSuccess) { - printf(" hipMalloc failed (%d) for %ld elems — skipping\n\n", (int)herr, nelems); + printf(" hipMalloc failed (%d) for %ld elems -- skipping\n\n", (int)herr, nelems); return; } hipMemset(dbuf, 0, nelems * sizeof(hipfftDoubleComplex)); - // 1. hipfftPlanMany (one-step, nullptr embed) — current Grid path + // 1. hipfftPlanMany (one-step, nullptr embed) -- current Grid path { hipfftHandle p; hipfftResult rv = hipfftPlanMany(&p, 1, n, @@ -79,7 +79,7 @@ static void tryPlanAndExec(int G, long howmany) { } } - // 2. hipfftCreate + hipfftMakePlanMany (two-step) — also current Grid path + // 2. hipfftCreate + hipfftMakePlanMany (two-step) -- also current Grid path { hipfftHandle p; size_t workSize = 0; diff --git a/tests/debug/Test_schur_dense_coarse.cc b/tests/debug/Test_schur_dense_coarse.cc index 506b30eb4..1eb98af8c 100644 --- a/tests/debug/Test_schur_dense_coarse.cc +++ b/tests/debug/Test_schur_dense_coarse.cc @@ -30,7 +30,7 @@ Author: Peter Boyle // The DenseCoarseMatrix GLUE test, on a real (tiny) lattice coarse // operator, CPU laptop build. // -// Builds a genuine GeneralCoarsenedMatrix (DWF MdagM + 0.5 shift for a +// Builds a genuine DeprecatedGeneralCoarsenedMatrix (DWF MdagM + 0.5 shift for a // guaranteed-invertible Galerkin coarse op, random aggregation basis) // and runs the whole Import certificate chain through the glue: // - fresh ImportDense + IMPORT CERTIFICATE vs Op.M @@ -143,7 +143,7 @@ int main (int argc, char ** argv) Aggregates.CreateSubspaceRandom(RNG5); std::cout << GridLogMessage << "Coarsening shifted MdagM" << std::endl; - typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef DeprecatedGeneralCoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d); LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d); diff --git a/tests/solver/Test_multigrid_common.h b/tests/solver/Test_multigrid_common.h index 7f4c14be0..75eafe617 100644 --- a/tests/solver/Test_multigrid_common.h +++ b/tests/solver/Test_multigrid_common.h @@ -397,7 +397,7 @@ public: MdagMLinearOperator coarseMdagMOp(_CoarseMatrix); std::cout << GridLogMG << " Level " << _CurrentLevel << ": **************************************************" << std::endl; - std::cout << GridLogMG << " Level " << _CurrentLevel << ": MG correctness check: 0 == (M - (Mdiag + Σ_μ Mdir_μ)) * v" << std::endl; + std::cout << GridLogMG << " Level " << _CurrentLevel << ": MG correctness check: 0 == (M - (Mdiag + Sum_mu Mdir_mu)) * v" << std::endl; std::cout << GridLogMG << " Level " << _CurrentLevel << ": **************************************************" << std::endl; random(_LevelInfo.PRNGs[_CurrentLevel], fineTmps[0]); @@ -406,22 +406,22 @@ public: fineMdagMOp.OpDiag(fineTmps[0], fineTmps[2]); // Mdiag * v fineTmps[4] = zero; - for(int dir = 0; dir < 4; dir++) { // Σ_μ Mdir_μ * v + for(int dir = 0; dir < 4; dir++) { // Sum_mu Mdir_mu * v for(auto disp : {+1, -1}) { fineMdagMOp.OpDir(fineTmps[0], fineTmps[3], dir, disp); fineTmps[4] = fineTmps[4] + fineTmps[3]; } } - fineTmps[5] = fineTmps[2] + fineTmps[4]; // (Mdiag + Σ_μ Mdir_μ) * v + fineTmps[5] = fineTmps[2] + fineTmps[4]; // (Mdiag + Sum_mu Mdir_mu) * v fineTmps[6] = fineTmps[1] - fineTmps[5]; auto deviation = std::sqrt(norm2(fineTmps[6]) / norm2(fineTmps[1])); std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2(M * v) = " << norm2(fineTmps[1]) << std::endl; std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2(Mdiag * v) = " << norm2(fineTmps[2]) << std::endl; - std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2(Σ_μ Mdir_μ * v) = " << norm2(fineTmps[4]) << std::endl; - std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2((Mdiag + Σ_μ Mdir_μ) * v) = " << norm2(fineTmps[5]) << std::endl; + std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2(Sum_mu Mdir_mu * v) = " << norm2(fineTmps[4]) << std::endl; + std::cout << GridLogMG << " Level " << _CurrentLevel << ": norm2((Mdiag + Sum_mu Mdir_mu) * v) = " << norm2(fineTmps[5]) << std::endl; std::cout << GridLogMG << " Level " << _CurrentLevel << ": relative deviation = " << deviation; if(deviation > tolerance) { @@ -664,7 +664,7 @@ createMGInstance(WilsonMGParams &mgParams, LevelInfo &levelInfo, Matrix &FineMat CASE_FOR_N_LEVELS(3); CASE_FOR_N_LEVELS(4); default: - std::cout << GridLogError << "We currently only support nLevels ∈ {2, 3, 4}" << std::endl; + std::cout << GridLogError << "We currently only support nLevels in {2, 3, 4}" << std::endl; exit(EXIT_FAILURE); break; }