diff --git a/Grid/GridCore.h b/Grid/GridCore.h index 51bcecd3a..251d32691 100644 --- a/Grid/GridCore.h +++ b/Grid/GridCore.h @@ -56,6 +56,7 @@ Author: paboyle #include #include #include +#include #include #include #include diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 6ebb1c58d..fb26692e3 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -33,6 +33,7 @@ NAMESPACE_CHECK(blas); #include NAMESPACE_CHECK(algorithms); +#include #include #include #include diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 986a46d02..e0eacc420 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -52,9 +52,31 @@ public: virtual void AdjOp (const Field &in, Field &out) = 0; // Abstract base virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0; virtual void HermOp(const Field &in, Field &out)=0; + + // As CheckerBoardedSparseMatrixBase::SplitClone; the returned bundle's Linop is the + // same kind of wrapper as this one. Collective. nullptr if not supported. + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + return nullptr; + } + virtual ~LinearOperatorBase(){}; }; +///////////////////////////////////////////////////////////////////////////////////////////// +// The wrappers below are templates on an arbitrary Matrix, and their virtual SplitClone is +// instantiated with the class, so it must compile for Matrix types with no SplitClone. +///////////////////////////////////////////////////////////////////////////////////////////// +template +SplitOperator *SplitCloneMatrix(Matrix &Mat,const Coordinate &mpi_split) +{ + if constexpr ( std::is_base_of,Matrix>::value ) { + return Mat.SplitClone(mpi_split); + } else { + return nullptr; + } +} + ///////////////////////////////////////////////////////////////////////////////////////////// // By sharing the class for Sparse Matrix across multiple operator wrappers, we can share code @@ -102,6 +124,17 @@ public: void HermOp(const Field &in, Field &out){ _Mat.MdagM(in,out); } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new MdagMLinearOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->Grid(); + return split; + } }; template class MMdagLinearOperator : public LinearOperatorBase { @@ -401,6 +434,17 @@ template _Mat.MooeeDag(in,out); axpy(out,-1.0,tmp,out); } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new SchurDiagMooeeOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template class SchurDiagOneOperator : public SchurOperatorBase { @@ -427,6 +471,17 @@ template _Mat.MeooeDag(out,tmp); axpy(out,-1.0,tmp,in); } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new SchurDiagOneOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template class SchurDiagTwoOperator : public SchurOperatorBase { @@ -455,6 +510,17 @@ template axpy(out,-1.0,tmp,in); } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new SchurDiagTwoOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template @@ -522,6 +588,17 @@ class NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase axpy(out, -1.0, tmp, out); } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new NonHermitianSchurDiagMooeeOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template @@ -552,6 +629,17 @@ class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new NonHermitianSchurDiagOneOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template @@ -583,6 +671,17 @@ class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new NonHermitianSchurDiagTwoOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; /////////////////////////////////////////////////////////////////////////////////////////////////// @@ -636,6 +735,17 @@ class SchurStaggeredOperator : public SchurOperatorBase { virtual void MpcDagMpc(const Field &in, Field &out) { GRID_ASSERT(0);// Never need with staggered } + + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + SplitOperator *split = SplitCloneMatrix(_Mat,mpi_split); + if ( split == nullptr ) { + return nullptr; + } + split->Linop = new SchurStaggeredOperator,Field>(*split->Matrix); + split->FieldGrid = split->Matrix->RedBlackGrid(); + return split; + } }; template using SchurStagOperator = SchurStaggeredOperator; diff --git a/Grid/algorithms/SparseMatrix.h b/Grid/algorithms/SparseMatrix.h index 2922b3e38..ccf750b04 100644 --- a/Grid/algorithms/SparseMatrix.h +++ b/Grid/algorithms/SparseMatrix.h @@ -78,6 +78,17 @@ public: virtual void MeooeDag (const Field &in, Field &out)=0; virtual void MooeeDag (const Field &in, Field &out)=0; virtual void MooeeInvDag (const Field &in, Field &out)=0; + + ////////////////////////////////////////////////////////////////////// + // Copy of this operator on grids whose communicator is split into + // partitions of MPI layout mpi_split. Collective over the full + // communicator. nullptr if this operator cannot be cloned. + ////////////////////////////////////////////////////////////////////// + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + return nullptr; + } + virtual ~CheckerBoardedSparseMatrixBase() {}; }; diff --git a/Grid/algorithms/SplitOperator.h b/Grid/algorithms/SplitOperator.h new file mode 100644 index 000000000..ba9e1589c --- /dev/null +++ b/Grid/algorithms/SplitOperator.h @@ -0,0 +1,172 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/SplitOperator.h + + 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 */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +template class LinearOperatorBase; +template class CheckerBoardedSparseMatrixBase; + +///////////////////////////////////////////////////////////////////////////////////////////// +// A copy of an operator living on grids whose communicator is split into independent +// partitions, together with those grids. Produced by SplitClone(). +// +// Owns everything it points to. Destruction is the reverse of creation: linear operator, +// matrix, then grids, so no object outlives a grid it references. For 4d operators the +// fermion grids are the gauge grids and are deleted once. +///////////////////////////////////////////////////////////////////////////////////////////// +template +class SplitOperator +{ +public: + GridCartesian *GaugeGrid = nullptr; + GridRedBlackCartesian *GaugeRBGrid = nullptr; + GridCartesian *FermionGrid = nullptr; + GridRedBlackCartesian *FermionRBGrid = nullptr; + + int Partition = 0; // Grid_split vector index held by this rank's partition + int Partitions = 1; // number of partitions + + CheckerBoardedSparseMatrixBase *Matrix = nullptr; + LinearOperatorBase *Linop = nullptr; + GridBase *FieldGrid = nullptr; // grid of Linop's fields + + SplitOperator(void) {}; + + SplitOperator(const SplitOperator &) = delete; + + SplitOperator &operator=(const SplitOperator &) = delete; + + ~SplitOperator(void) + { + delete Linop; + delete Matrix; + if ( FermionRBGrid != GaugeRBGrid ) { + delete FermionRBGrid; + } + if ( FermionGrid != GaugeGrid ) { + delete FermionGrid; + } + delete GaugeRBGrid; + delete GaugeGrid; + } +}; + +///////////////////////////////////////////////////////////////////////////////////////////// +// Index of the vector that Grid_split(std::vector full, Field split) delivers to this +// rank's partition. Grid_split orders partitions lexicographically with the first +// dimension fastest; the split communicator's own rank (srank) uses the reversed MPI +// convention, so the two differ whenever more than one dimension is split. +///////////////////////////////////////////////////////////////////////////////////////////// +inline int GridSplitVectorIndex(GridBase *full,GridBase *split) +{ + int nd = full->_ndimension; + GRID_ASSERT(split->_ndimension == nd); + + Coordinate scoor(nd); + Coordinate ssize(nd); + for(int d=0;dThisProcessorCoor()[d] / split->ProcessorGrid()[d]; + ssize[d] = full->ProcessorGrid()[d] / split->ProcessorGrid()[d]; + } + int index; + Lexicographic::IndexFromCoor(scoor,index,ssize); + return index; +} + +///////////////////////////////////////////////////////////////////////////////////////////// +// Partition MPI layout for a batched solve on grid, from a request as given on the command +// line (--batched-solver-split). Uses the trailing dimensions of grid's processor and shm +// layouts, so it works for 4d and 5d grids. Returns grid's own processor layout (one +// partition) when no split is requested. Asserts divisibility; warns when partitions +// straddle nodes or when node boundaries are not visible. +///////////////////////////////////////////////////////////////////////////////////////////// +inline Coordinate BatchedSolverSplitLayout(GridBase *grid, + const Coordinate &request, + bool node, + int &partitions) +{ + int nd = GridDefaultMpi().size(); + int pad = grid->_ndimension - nd; + GRID_ASSERT(pad >= 0); + GRID_ASSERT(grid->ShmGrid().size() == grid->_ndimension); + + Coordinate processors(nd); + Coordinate shm(nd); + for(int d=0;dProcessorGrid()[pad+d]; + shm[d] = grid->ShmGrid()[pad+d]; + } + + Coordinate split(nd); + if ( node ) { + split = shm; + } else if ( request.size() == 0 ) { + split = processors; + } else { + GRID_ASSERT(request.size() == nd); + split = request; + } + + partitions = 1; + for(int d=0;d *guesser; bool updateResidual; + // Inner solves on independent partitions of the communicator (--batched-solver-split). + // BatchedSplit is the partition MPI layout, empty for no split; BatchedSplitNode selects + // one partition per node. Fixed at construction, default from the command line. + const Coordinate BatchedSplit; + const bool BatchedSplitNode; + + // Copy of Linop_f on the partitions, made once at construction and owned; nullptr when + // there is no split or the operator cannot be split. It holds the gauge field as it was at + // construction: a solver must not outlive a change of the operator's gauge field. The + // outer defect correction uses Linop_d, so a stale copy would cost convergence, not + // correctness. + SplitOperator *split = nullptr; + int Partitions = 1; + MixedPrecisionConjugateGradientBatched(RealD tol, Integer maxinnerit, Integer maxouterit, @@ -58,10 +72,37 @@ public: GridBase* _sp_grid, LinearOperatorBase &_Linop_f, LinearOperatorBase &_Linop_d, - bool _updateResidual=true) : + bool _updateResidual=true, + const Coordinate &_split=GridDefaultBatchedSolverSplit(), + bool _split_node=GridDefaultBatchedSolverSplitNode()) : Linop_f(_Linop_f), Linop_d(_Linop_d), Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), MaxPatchupIterations(maxpatchit), SinglePrecGrid(_sp_grid), - OuterLoopNormMult(100.), guesser(NULL), updateResidual(_updateResidual) { }; + OuterLoopNormMult(100.), guesser(NULL), updateResidual(_updateResidual), + BatchedSplit(_split), BatchedSplitNode(_split_node) + { + Coordinate layout = BatchedSolverSplitLayout(SinglePrecGrid,BatchedSplit,BatchedSplitNode,Partitions); + if ( Partitions > 1 ) { + double t0 = usecond(); + split = Linop_f.SplitClone(layout); + double t1 = usecond(); + if ( split == nullptr ) { + std::cout << GridLogMessage << "MixedPrecisionConjugateGradientBatched: operator cannot be split; serial inner solves" << std::endl; + } else { + std::cout << GridLogMessage << "MixedPrecisionConjugateGradientBatched: split clone " << (t1-t0)/1.0e6 << " s" << std::endl; + HostMemoryReport(SinglePrecGrid,GridLogMessage,"MixedPrecisionConjugateGradientBatched: after split clone"); + } + } + }; + + // Owns split + MixedPrecisionConjugateGradientBatched(const MixedPrecisionConjugateGradientBatched &) = delete; + + MixedPrecisionConjugateGradientBatched &operator=(const MixedPrecisionConjugateGradientBatched &) = delete; + + ~MixedPrecisionConjugateGradientBatched(void) + { + delete split; + } void useGuesser(LinearFunction &g){ guesser = &g; @@ -91,6 +132,8 @@ public: GridStopWatch InnerCGtimer; GridStopWatch PrecChangeTimer; + GridStopWatch OuterResidualTimer; + GridStopWatch PatchupTimer; int cb = src_d_in[0].Checkerboard(); @@ -130,18 +173,27 @@ public: ConjugateGradient CG_f(inner_tol, MaxInnerIterations); CG_f.ErrorOnNoConverge = false; + SplitTimers splitTimers; + if ( split != nullptr && split->Partitions > NBatch ) { + std::cout << GridLogWarning << "MixedPrecisionConjugateGradientBatched: " << split->Partitions << " partitions for " + << NBatch << " right-hand sides; the extra partitions only solve zero padding" << std::endl; + } + Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){ std::cout << GridLogMessage << std::endl; std::cout << GridLogMessage << "Outer iteration " << outer_iter << std::endl; + HostMemoryReport(DoublePrecGrid,GridLogMessage,"MixedPrecisionConjugateGradientBatched: outer iteration "+std::to_string(outer_iter)); bool allConverged = true; for (int i=0; i CG_d(Tolerance, MaxPatchupIterations); CG_d(Linop_d, src_d_in[i], sol_d[i]); TotalFinalStepIterations[i] += CG_d.IterationsToComplete; } + PatchupTimer.Stop(); TotalTimer.Stop(); @@ -204,8 +264,118 @@ public: } std::cout << GridLogMessage << std::endl; std::cout< &split, + ConjugateGradient &CG_f, + std::vector &src_f, + std::vector &sol_f, + std::vector &TotalInnerIterations, + GridStopWatch &InnerCGtimer, + SplitTimers &timers) + { + int NBatch = src_f.size(); + int P = split.Partitions; + int cb = src_f[0].Checkerboard(); + + FieldF s_src(split.FieldGrid); + FieldF s_sol(split.FieldGrid); + std::vector group_src(P,SinglePrecGrid); + std::vector group_sol(P,SinglePrecGrid); + std::vector iters(P); + + for(int g=0;gThisRank() == 0 ) { + iters[split.Partition] = CG_f.IterationsToComplete; + } + SinglePrecGrid->GlobalSumVector(&iters[0],P); + timers.Reduce.Stop(); + + // Includes the host to device copy of group_sol left by Grid_unsplit + timers.Staging.Start(); + for(int p=0;p

&out, const Lattice &in){ //////////////////////////////////////////////////////////////////////////////// // Communicate between grids //////////////////////////////////////////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////////////////////////// -// SIMPLE CASE: -/////////////////////////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////////// +// Grid_split / Grid_unsplit // -// Mesh of nodes (2x2) ; subdivide to 1x1 subdivisions +// The split grid divides the full grid's processor layout Pf into partitions of layout +// Ps; ratio R = Pf/Ps per dimension, nvector = prod R partitions. A rank at full processor +// coordinate c belongs to partition q = c/Ps at split coordinate s = c%Ps, and partition q +// holds vector index lex(q,R) (GridSplitVectorIndex). Local volumes satisfy +// Ls = Lf * R per dimension. // -// Lex ord: -// N0 va0 vb0 vc0 vd0 N1 va1 vb1 vc1 vd1 -// N2 va2 vb2 vc2 vd2 N3 va3 vb3 vc3 vd3 +// Every full-grid rank's local block of field v lies wholly inside one split-grid rank's +// local volume, in the sub-box at offset o = c%R (units of Lf), on the rank of partition +// q(v) at split coordinate c/R. The redistribution is therefore nvector whole-block +// point-to-point exchanges, nothing finer, and no data moves twice. Step k (a coordinate in +// the R box) pairs every rank with one sender and one receiver: // -// Ratio = full[dim] / split[dim] +// send field q = (c%R + k)%R to full coordinate q*Ps + c/R +// recv sub-box o = (c/Ps - k)%R from full coordinate (c%Ps)*R + o // -// For each dimension do an all to all; get Nvec -> Nvec / ratio -// Ldim -> Ldim * ratio -// LocalVol -> LocalVol * ratio -// full AllToAll(0) -// N0 va0 vb0 va1 vb1 N1 vc0 vd0 vc1 vd1 -// N2 va2 vb2 va3 vb3 N3 vc2 vd2 vc3 vd3 +// Within a step the send map is a permutation of the ranks and the receive map is its +// inverse, so each step is one symmetric SendToRecvFrom per rank; a rank paired with +// itself copies locally. Grid_unsplit runs the same steps with the roles exchanged. // -// REARRANGE -// N0 va01 vb01 N1 vc01 vd01 -// N2 va23 vb23 N3 vc23 vd23 -// -// full AllToAll(1) // Not what is wanted. FIXME -// N0 va01 va23 N1 vc01 vc23 -// N2 vb01 vb23 N3 vd01 vd23 -// -// REARRANGE -// N0 va0123 N1 vc0123 -// N2 vb0123 N3 vd0123 -// -// Must also rearrange data to get into the NEW lex order of grid at each stage. Some kind of "insert/extract". -// NB: Easiest to programme if keep in lex order. -/* - * Let chunk = (fvol*nvec)/sP be size of a chunk. ( Divide lexico vol * nvec into fP/sP = M chunks ) - * - * 2nd A2A (over sP nodes; subdivide the fP into sP chunks of M) - * - * node 0 1st chunk of node 0M..(1M-1); 2nd chunk of node 0M..(1M-1).. data chunk x M x sP = fL / sP * M * sP = fL * M growth - * node 1 1st chunk of node 1M..(2M-1); 2nd chunk of node 1M..(2M-1).. - * node 2 1st chunk of node 2M..(3M-1); 2nd chunk of node 2M..(3M-1).. - * node 3 1st chunk of node 3M..(3M-1); 2nd chunk of node 2M..(3M-1).. - * etc... - */ +// Data stays on the device. A block on the full grid is the local volume in lexicographic +// order; the split grid's local volume is staged box-major (sub-box index, then +// lexicographic within the sub-box), so every exchanged block is contiguous. Packing and +// unpacking are accelerator kernels. Buffers are deviceVectors: the staged split volume +// plus one block. Without ACCELERATOR_AWARE_MPI each block is staged through host memory. +///////////////////////////////////////////////////////////////////////////////////////////// + +// Scalar sites of in, written to buf in box-major order: sub-boxes of extent bdims taken +// lexicographically, sites lexicographically within each. bdims = local dimensions gives +// plain lexicographic order. +template +void GridSplitUnvectorise(const Lattice &in,typename vobj::scalar_object *buf,const Coordinate &bdims) +{ + GridBase *grid = in.Grid(); + int nd = grid->_ndimension; + Coordinate rdims = grid->_rdimensions; + Coordinate simd = grid->_simd_layout; + Coordinate nbox(nd); + int64_t bvol = 1; + for(int d=0;d_ldimensions[d] % bdims[d]) == 0 ); + nbox[d] = grid->_ldimensions[d] / bdims[d]; + bvol *= bdims[d]; + } + const int Nsimd = grid->Nsimd(); + autoView( in_v , in, AcceleratorRead); + accelerator_for(ss, grid->oSites(), 1, { + Coordinate ocoor(nd); + Coordinate icoor(nd); + Coordinate bcoor(nd); + Coordinate wcoor(nd); + Lexicographic::CoorFromIndex(ocoor,ss,rdims); + for(int lane=0;lane +void GridSplitVectorise(const typename vobj::scalar_object *buf,Lattice &out,const Coordinate &bdims) +{ + GridBase *grid = out.Grid(); + int nd = grid->_ndimension; + Coordinate rdims = grid->_rdimensions; + Coordinate simd = grid->_simd_layout; + Coordinate nbox(nd); + int64_t bvol = 1; + for(int d=0;d_ldimensions[d] % bdims[d]) == 0 ); + nbox[d] = grid->_ldimensions[d] / bdims[d]; + bvol *= bdims[d]; + } + const int Nsimd = grid->Nsimd(); + autoView( out_v , out, AcceleratorWriteDiscard); + accelerator_for(ss, grid->oSites(), 1, { + Coordinate ocoor(nd); + Coordinate icoor(nd); + Coordinate bcoor(nd); + Coordinate wcoor(nd); + Lexicographic::CoorFromIndex(ocoor,ss,rdims); + for(int lane=0;lane Field; + std::vector FieldRank; + std::vector Box; + std::vector BoxRank; + + GridSplitSchedule(GridBase *_full_grid,GridBase *split_grid) : full_grid(_full_grid) + { + nd = full_grid->_ndimension; + GRID_ASSERT(split_grid->_ndimension == nd); + ratio.resize(nd); + fldims = full_grid->_ldimensions; + Coordinate sprocs = split_grid->_processors; + Coordinate fcoor = full_grid->_processor_coor; + nvector = 1; + for(int d=0;d_gdimensions[d] == split_grid->_gdimensions[d]); + GRID_ASSERT(full_grid->_fdimensions[d] == split_grid->_fdimensions[d]); + GRID_ASSERT( (full_grid->_processors[d] % sprocs[d]) == 0 ); + ratio[d] = full_grid->_processors[d] / sprocs[d]; + GRID_ASSERT(split_grid->_ldimensions[d] == fldims[d]*ratio[d]); + nvector *= ratio[d]; + } + GRID_ASSERT(nvector*split_grid->_Nprocessors == full_grid->_Nprocessors); + + Field.resize(nvector); + FieldRank.resize(nvector); + Box.resize(nvector); + BoxRank.resize(nvector); + for(int k=0;kRankFromProcessorCoor(dest); + BoxRank[k] = full_grid->RankFromProcessorCoor(from); + } + } + + // One step's exchange; returns the bytes that left this rank + uint64_t Exchange(void *xmit,int dest,void *recv,int from,uint64_t bytes) + { + int me = full_grid->ThisRank(); + // The send map is a permutation and the receive map its inverse: self pairs with self + GRID_ASSERT( (dest == me) == (from == me) ); + if ( dest == me ) { + acceleratorCopyDeviceToDevice(xmit,recv,bytes); + return 0; + } + GRID_ASSERT( (bytes % sizeof(int32_t)) == 0 ); + GRID_ASSERT( bytes/sizeof(int32_t) < (1ULL<<31) ); // SendToRecvFrom counts int32 words +#ifdef ACCELERATOR_AWARE_MPI + full_grid->SendToRecvFrom(xmit,dest,recv,from,bytes); +#else + std::vector hxmit(bytes); + std::vector hrecv(bytes); + acceleratorCopyFromDevice(xmit,&hxmit[0],bytes); + full_grid->SendToRecvFrom(&hxmit[0],dest,&hrecv[0],from,bytes); + acceleratorCopyToDevice(&hrecv[0],recv,bytes); +#endif + return bytes; + } +}; + +// Common body: full[v] may alias one another (the single field overload) template -void Grid_split(std::vector > & full,Lattice & split) +void Grid_split(const std::vector *> &full,Lattice &split) { typedef typename Vobj::scalar_object Sobj; - int full_vecs = full.size(); - - GRID_ASSERT(full_vecs>=1); - - GridBase * full_grid = full[0].Grid(); + GRID_ASSERT(full.size()>=1); + GridBase *full_grid = full[0]->Grid(); GridBase *split_grid = split.Grid(); - int ndim = full_grid->_ndimension; - int full_nproc = full_grid->_Nprocessors; - int split_nproc =split_grid->_Nprocessors; - - //////////////////////////////// - // Checkerboard management - //////////////////////////////// - int cb = full[0].Checkerboard(); + int cb = full[0]->Checkerboard(); split.Checkerboard() = cb; - - ////////////////////////////// - // Checks - ////////////////////////////// - GRID_ASSERT(full_grid->_ndimension==split_grid->_ndimension); - for(int n=0;n_gdimensions[d]==split.Grid()->_gdimensions[d]); - GRID_ASSERT(full[n].Grid()->_fdimensions[d]==split.Grid()->_fdimensions[d]); - } + for(int n=0;nGrid() == full_grid); + GRID_ASSERT(full[n]->Checkerboard() == cb); } - int nvector =full_nproc/split_nproc; - GRID_ASSERT(nvector*split_nproc==full_nproc); - GRID_ASSERT(nvector == full_vecs); + double t_start = usecond(); + GridSplitSchedule sched(full_grid,split_grid); + GRID_ASSERT(sched.nvector == full.size()); - Coordinate ratio(ndim); - for(int d=0;d_processors[d]/ split_grid->_processors[d]; + uint64_t fvol = full_grid->lSites(); + uint64_t bytes = fvol*sizeof(Sobj); + deviceVector staged(split_grid->lSites()); + deviceVector block(fvol); + Sobj *staged_p = &staged[0]; + Sobj *block_p = &block[0]; + double t_alloc = usecond() - t_start; + + double t_pack = 0; + double t_comm = 0; + uint64_t offrank = 0; + for(int k=0;k_ldimensions); + double t1 = usecond(); + offrank += sched.Exchange((void *)block_p,sched.FieldRank[k], + (void *)&staged_p[sched.Box[k]*fvol],sched.BoxRank[k],bytes); + double t2 = usecond(); + t_pack += t1-t0; + t_comm += t2-t1; } - uint64_t lsites = full_grid->lSites(); - uint64_t sz = lsites * nvector; - std::vector tmpdata(sz); - std::vector alldata(sz); - std::vector scalardata(lsites); + double t0 = usecond(); + GridSplitVectorise(staged_p,split,sched.fldims); + double t_unpack = usecond() - t0; + double t_total = usecond() - t_start; - for(int v=0;v_ldimensions; - - for(int d=ndim-1;d>=0;d--){ - - if ( ratio[d] != 1 ) { - - full_grid ->AllToAll(d,alldata,tmpdata); - if ( split_grid->_processors[d] > 1 ) { - alldata=tmpdata; - split_grid->AllToAll(d,alldata,tmpdata); - } - - auto rdims = ldims; - auto M = ratio[d]; - auto rsites= lsites*M;// increases rsites by M - nvec /= M; // Reduce nvec by subdivision factor - rdims[d] *= M; // increase local dim by same factor - - int sP = split_grid->_processors[d]; - int fP = full_grid->_processors[d]; - - int fvol = lsites; - - int chunk = (nvec*fvol)/sP; GRID_ASSERT(chunk*sP == nvec*fvol); - - // Loop over reordered data post A2A - thread_for(c, chunk, { - Coordinate coor(ndim); - for(int m=0;m +void Grid_split(std::vector > & full,Lattice & split) +{ + std::vector *> full_p(full.size()); + for(int n=0;n void Grid_split(Lattice &full,Lattice & split) { int nvector = full.Grid()->_Nprocessors / split.Grid()->_Nprocessors; - std::vector > full_v(nvector,full.Grid()); - for(int n=0;n *> full_p(nvector,&full); + Grid_split(full_p,split); } template @@ -1725,119 +1828,57 @@ void Grid_unsplit(std::vector > & full,Lattice & split) { typedef typename Vobj::scalar_object Sobj; - int full_vecs = full.size(); - - GRID_ASSERT(full_vecs>=1); - - GridBase * full_grid = full[0].Grid(); + GRID_ASSERT(full.size()>=1); + GridBase *full_grid = full[0].Grid(); GridBase *split_grid = split.Grid(); - int ndim = full_grid->_ndimension; - int full_nproc = full_grid->_Nprocessors; - int split_nproc =split_grid->_Nprocessors; - - //////////////////////////////// - // Checkerboard management - //////////////////////////////// int cb = full[0].Checkerboard(); split.Checkerboard() = cb; - - ////////////////////////////// - // Checks - ////////////////////////////// - GRID_ASSERT(full_grid->_ndimension==split_grid->_ndimension); - for(int n=0;n_gdimensions[d]==split.Grid()->_gdimensions[d]); - GRID_ASSERT(full[n].Grid()->_fdimensions[d]==split.Grid()->_fdimensions[d]); - } } - int nvector =full_nproc/split_nproc; - GRID_ASSERT(nvector*split_nproc==full_nproc); - GRID_ASSERT(nvector == full_vecs); + double t_start = usecond(); + GridSplitSchedule sched(full_grid,split_grid); + GRID_ASSERT(sched.nvector == full.size()); - Coordinate ratio(ndim); - for(int d=0;d_processors[d]/ split_grid->_processors[d]; + uint64_t fvol = full_grid->lSites(); + uint64_t bytes = fvol*sizeof(Sobj); + deviceVector staged(split_grid->lSites()); + deviceVector block(fvol); + Sobj *staged_p = &staged[0]; + Sobj *block_p = &block[0]; + double t_alloc = usecond() - t_start; + + double t0 = usecond(); + GridSplitUnvectorise(split,staged_p,sched.fldims); + double t_pack = usecond() - t0; + + double t_comm = 0; + double t_unpack = 0; + uint64_t offrank = 0; + for(int k=0;k_ldimensions); + double t3 = usecond(); + t_comm += t2-t1; + t_unpack += t3-t2; } + double t_total = usecond() - t_start; - uint64_t lsites = full_grid->lSites(); - uint64_t sz = lsites * nvector; - std::vector tmpdata(sz); - std::vector alldata(sz); - std::vector scalardata(lsites); - - unvectorizeToLexOrdArray(alldata,split); - - ///////////////////////////////////////////////////////////////// - // Start from split grid and work towards full grid - ///////////////////////////////////////////////////////////////// - - int nvec = 1; - uint64_t rsites = split_grid->lSites(); - Coordinate rdims = split_grid->_ldimensions; - - for(int d=0;d_processors[d]; - int fP = full_grid->_processors[d]; - - auto ldims = rdims; ldims[d] /= M; // Decrease local dims by same factor - auto lsites= rsites/M; // Decreases rsites by M - - int fvol = lsites; - int chunk = (nvec*fvol)/sP; GRID_ASSERT(chunk*sP == nvec*fvol); - - { - // Loop over reordered data post A2A - thread_for(c, chunk,{ - Coordinate coor(ndim); - for(int m=0;m_processors[d] > 1 ) { - split_grid->AllToAll(d,tmpdata,alldata); - tmpdata=alldata; - } - full_grid ->AllToAll(d,tmpdata,alldata); - rdims[d]/= M; - rsites /= M; - nvec *= M; // Increase nvec by subdivision factor - } - } - - lsites = full_grid->lSites(); - for(int v=0;v + + 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 + +#include +#include +#include +#include +#ifdef __APPLE__ +#include +#endif + +NAMESPACE_BEGIN(Grid); + +///////////////////////////////////////////////////////////////////////////////////////////// +// Host memory as the OOM killer sees it, plus the device allocator footprint. +// Targets are Apple or Linux; on Apple only the RSS and device fields are filled. +// +// RSS : this process, now and its high-water mark (getrusage). +// Node : MemAvailable and MemTotal from /proc/meminfo (Linux). +// Cgroup : the nearest enclosing memory cgroup with a finite limit (the one the batch +// system enforces), its usage now, its high-water mark and the limit (Linux, +// v2 or v1 hierarchy). Cgroup usage includes page cache, tmpfs and kernel memory +// charged to the job, none of which appear in any process RSS. +// Device : MemoryManager Lattice footprint on the device and its free-block cache. Device +// allocations outside the MemoryManager (comms buffers) are not included. +// +// All values in GB; -1 where the source does not exist on this platform or kernel. +///////////////////////////////////////////////////////////////////////////////////////////// +struct HostMemoryStatus +{ + RealD RSS = -1; + RealD RSSPeak = -1; + RealD NodeAvailable = -1; + RealD NodeTotal = -1; + RealD CgroupCurrent = -1; + RealD CgroupPeak = -1; + RealD CgroupLimit = -1; + RealD DeviceLattice = -1; + RealD DeviceCache = -1; +}; + +#ifndef __APPLE__ // Linux + +// One number from a file; -1 if the file is absent or holds no number (e.g. "max") +inline RealD HostMemoryReadBytes(const std::string &file) +{ + std::ifstream f(file); + if ( !f.good() ) { + return -1; + } + std::string word; + f >> word; + if ( word.empty() || !isdigit(word[0]) ) { + return -1; + } + return std::stod(word); +} + +// Value of a "Key: N kB" line in /proc/meminfo, in bytes; -1 if absent +inline RealD HostMemoryMeminfo(const std::string &key) +{ + std::ifstream f("/proc/meminfo"); + std::string name; + RealD value; + std::string unit; + while ( f >> name >> value >> unit ) { + if ( name == key + ":" ) { + return value*1024.0; + } + } + return -1; +} + +// Fill the cgroup fields from the nearest enclosing memory cgroup with a finite limit, or +// from this process's own cgroup (limit left at -1) if no ancestor has one +inline void HostMemoryCgroup(HostMemoryStatus &s) +{ + std::ifstream f("/proc/self/cgroup"); + std::string line; + std::string root; + std::string path; + std::string current; + std::string peak; + std::string limit; + while ( std::getline(f,line) ) { + size_t c1 = line.find(':'); + size_t c2 = line.find(':',c1+1); + if ( c1 == std::string::npos || c2 == std::string::npos ) { + continue; + } + std::string controllers = line.substr(c1+1,c2-c1-1); + // v2 unified hierarchy: "0::/path" + if ( controllers.empty() && root.empty() ) { + root = "/sys/fs/cgroup"; + path = line.substr(c2+1); + current = "memory.current"; + peak = "memory.peak"; + limit = "memory.max"; + } + // v1 memory controller: "N:memory:/path", takes precedence + if ( (","+controllers+",").find(",memory,") != std::string::npos ) { + root = "/sys/fs/cgroup/memory"; + path = line.substr(c2+1); + current = "memory.usage_in_bytes"; + peak = "memory.max_usage_in_bytes"; + limit = "memory.limit_in_bytes"; + break; + } + } + if ( root.empty() ) { + return; + } + + // v1 reports "no limit" as a number near 2^63 + const RealD unlimited = 1.0e18; + std::string dir = path; + while ( true ) { + std::string base = root + dir + "/"; + RealD lim = HostMemoryReadBytes(base+limit); + if ( lim > 0 && lim < unlimited ) { + s.CgroupCurrent = HostMemoryReadBytes(base+current)/1.0e9; + s.CgroupPeak = HostMemoryReadBytes(base+peak)/1.0e9; + s.CgroupLimit = lim/1.0e9; + return; + } + if ( dir.empty() || dir == "/" ) { + break; + } + dir = dir.substr(0,dir.rfind('/')); + } + std::string base = root + path + "/"; + s.CgroupCurrent = HostMemoryReadBytes(base+current)/1.0e9; + s.CgroupPeak = HostMemoryReadBytes(base+peak)/1.0e9; +} + +#endif + +// This rank only; no communication +inline HostMemoryStatus HostMemoryQuery(void) +{ + HostMemoryStatus s; + + struct rusage ru; + getrusage(RUSAGE_SELF,&ru); +#ifdef __APPLE__ + s.RSSPeak = ru.ru_maxrss/1.0e9; // bytes on macOS + mach_task_basic_info_data_t info; + mach_msg_type_number_t count = MACH_TASK_BASIC_INFO_COUNT; + if ( task_info(mach_task_self(),MACH_TASK_BASIC_INFO,(task_info_t)&info,&count) == KERN_SUCCESS ) { + s.RSS = info.resident_size/1.0e9; + } +#else // Linux + s.RSSPeak = ru.ru_maxrss*1024.0/1.0e9; // kilobytes on Linux + std::ifstream statm("/proc/self/statm"); + long pages = 0; + long resident = 0; + if ( statm >> pages >> resident ) { + s.RSS = resident*(RealD)sysconf(_SC_PAGESIZE)/1.0e9; + } + RealD avail = HostMemoryMeminfo("MemAvailable"); + RealD total = HostMemoryMeminfo("MemTotal"); + if ( avail >= 0 ) { + s.NodeAvailable = avail/1.0e9; + } + if ( total >= 0 ) { + s.NodeTotal = total/1.0e9; + } + HostMemoryCgroup(s); +#endif + + s.DeviceLattice = MemoryManager::DeviceBytes/1.0e9; + s.DeviceCache = MemoryManager::DeviceCacheBytes()/1.0e9; + return s; +} + +///////////////////////////////////////////////////////////////////////////////////////////// +// Collective over grid. One line on log: maxima over ranks, except node MemAvailable which +// is the minimum over ranks (the tightest node). +///////////////////////////////////////////////////////////////////////////////////////////// +inline void HostMemoryReport(GridBase *grid,GridLogger &log,const std::string &label) +{ + HostMemoryStatus s = HostMemoryQuery(); + RealD minus_available = -s.NodeAvailable; + grid->GlobalMax(s.RSS); + grid->GlobalMax(s.RSSPeak); + grid->GlobalMax(minus_available); + grid->GlobalMax(s.NodeTotal); + grid->GlobalMax(s.CgroupCurrent); + grid->GlobalMax(s.CgroupPeak); + grid->GlobalMax(s.CgroupLimit); + grid->GlobalMax(s.DeviceLattice); + grid->GlobalMax(s.DeviceCache); + std::cout << log << "HOSTMEM " << label + << " : rank RSS " << s.RSS << " peak " << s.RSSPeak + << " ; node available " << -minus_available << " of " << s.NodeTotal + << " ; cgroup " << s.CgroupCurrent << " peak " << s.CgroupPeak << " limit " << s.CgroupLimit + << " ; device lattice " << s.DeviceLattice << " cache " << s.DeviceCache + << " GB" << std::endl; +} + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index 77491def2..11019b4b2 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -81,6 +81,15 @@ public: mass_minus=_mass_minus; SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs } ; + + // Everything SetMass needs, plus the boost, then rebuild the coefficient arrays in clone + void CloneCoefficientsInto(CayleyFermion5D &clone) + { + clone.mass_plus = mass_plus; + clone.mass_minus = mass_minus; + clone.qmu = qmu; + clone.SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); + } void P(const FermionField &psi, FermionField &chi); void Pdag(const FermionField &psi, FermionField &chi); diff --git a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h index 18287d8aa..aa130007f 100644 --- a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h @@ -189,6 +189,71 @@ public: void ImportGauge(const GaugeField& _Umu) override; + // Install already computed clover fields and inverses, and derive the checkerboarded parts + void ImportCloverTerms(const CloverDiagonalField& D, + const CloverTriangleField& T, + const CloverDiagonalField& Dinv, + const CloverTriangleField& Tinv) + { + Diagonal = D; + Triangle = T; + DiagonalInv = Dinv; + TriangleInv = Tinv; + pickCheckerboard(Even, DiagonalEven, Diagonal); + pickCheckerboard(Even, TriangleEven, Triangle); + pickCheckerboard(Odd, DiagonalOdd, Diagonal); + pickCheckerboard(Odd, TriangleOdd, Triangle); + pickCheckerboard(Even, DiagonalInvEven, DiagonalInv); + pickCheckerboard(Even, TriangleInvEven, TriangleInv); + pickCheckerboard(Odd, DiagonalInvOdd, DiagonalInv); + pickCheckerboard(Odd, TriangleInvOdd, TriangleInv); + } + + // Exact type only: derived operators must not inherit this + SplitOperator *SplitClone(const Coordinate& mpi_split) override + { + if ( typeid(*this) != typeid(CompactWilsonCloverFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field and clover fields are overwritten below. + // Boundary masks are rebuilt by the constructor from global coordinates. + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + CompactWilsonCloverFermion *clone = + new CompactWilsonCloverFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass, + csw_r, + csw_t, + cF, + this->anisotropyCoeff, + this->Params); + // The stored coefficients are post-scaled; the constructor would scale them again + clone->csw_r = csw_r; + clone->csw_t = csw_t; + + CloverDiagonalField sD(split->GaugeGrid); + CloverTriangleField sT(split->GaugeGrid); + CloverDiagonalField sDinv(split->GaugeGrid); + CloverTriangleField sTinv(split->GaugeGrid); + Grid_split(Diagonal,sD); + Grid_split(Triangle,sT); + Grid_split(DiagonalInv,sDinv); + Grid_split(TriangleInv,sTinv); + clone->ImportCloverTerms(sD,sT,sDinv,sTinv); + + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + ///////////////////////////////////////////// // Helpers ///////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/DWFSlow.h b/Grid/qcd/action/fermion/DWFSlow.h index 0ffc3ec63..b926dccff 100644 --- a/Grid/qcd/action/fermion/DWFSlow.h +++ b/Grid/qcd/action/fermion/DWFSlow.h @@ -109,10 +109,69 @@ public: //////////////////////// // Derivative interface //////////////////////// - // Interface calls an internal routine - void DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { GRID_ASSERT(0);}; - void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){ GRID_ASSERT(0);}; - void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){ GRID_ASSERT(0);}; + // Force on every link of every s-slice, in the form WilsonFermion5D uses: + // F5[mu](x,s) = tr_spin[ Btilde(x,s) A(x,s)^dag ] + // Btilde(x,s) = -1/2 U_mu(x,s) (1 -+ gamma_mu) B(x+mu,s), lower sign for dag + // Dhop5 carries no gauge links and contributes nothing. + void DhopDeriv5D(std::vector &F5,const FermionField &A,const FermionField &B,int dag) + { + RealD sgn= 1.0; + if (dag ) sgn=-1.0; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + FermionField tmp(B.Grid()); + FermionField Btilde(B.Grid()); + for(int mu=0;mu(outerProduct(Btilde,A)); + } + } + + // Every slice holds the same 4D gauge field, so the 4D force is the sum over slices + void DhopDeriv(GaugeField &mat,const FermionField &A,const FermionField &B,int dag) + { + std::vector F5(Nd,_grid); + DhopDeriv5D(F5,A,B,dag); + + GaugeLinkField slice(_grid4); + GaugeLinkField sum(_grid4); + for(int mu=0;mu(mat,sum,mu); + } + } + + // The checkerboard here is 5D, so a field on one checkerboard touches links at both 4D + // parities: the even-odd forces are returned on the full 4D gauge grid. + int CheckerboardedForce(void) { return 0; }; + + // A and B on opposite 5D checkerboards, placed in zeroed full fields: the full-lattice + // derivative then picks up exactly the links Meooe uses between them. + void DhopDerivCB(GaugeField &mat,const FermionField &A,const FermionField &B,int dag) + { + GRID_ASSERT(A.Checkerboard() != B.Checkerboard()); + FermionField Afull(_grid); + FermionField Bfull(_grid); + Afull = Zero(); + Bfull = Zero(); + setCheckerboard(Afull,A); + setCheckerboard(Bfull,B); + DhopDeriv(mat,Afull,Bfull,dag); + } + void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){ DhopDerivCB(mat,U,V,dag);}; + void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){ DhopDerivCB(mat,U,V,dag);}; /////////////////////////////////////////////////////////////// // non-hermitian hopping term; half cb or both @@ -200,6 +259,7 @@ public: _grid(&Fgrid), _cbgrid(&Hgrid), _grid4(_Umu.Grid()), + _cbgrid4(nullptr), // the checkerboard is 5D; there is no 4D red-black grid Umu(Nd,&Fgrid), UmuEven(Nd,&Hgrid), UmuOdd(Nd,&Hgrid), @@ -221,10 +281,6 @@ public: MassField =scalar(-mass); one =scalar(1.0); MassField =where(coor==Integer(Ls-1),MassField,one); - for(int mu=0;muLs;s++){ InsertSlice(U4,Umu[mu],s,0); } + pickCheckerboard(Even,UmuEven[mu],Umu[mu]); + pickCheckerboard(Odd ,UmuOdd[mu],Umu[mu]); } } diff --git a/Grid/qcd/action/fermion/DomainWallFermion.h b/Grid/qcd/action/fermion/DomainWallFermion.h index 680c203eb..38fbcdc36 100644 --- a/Grid/qcd/action/fermion/DomainWallFermion.h +++ b/Grid/qcd/action/fermion/DomainWallFermion.h @@ -132,6 +132,35 @@ public: Approx::zolotarev_free(zdata); } + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(DomainWallFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + DomainWallFermion *clone = new DomainWallFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index a2abaa781..802f1604b 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -59,6 +59,38 @@ public: virtual GridBase *GaugeGrid(void) =0; virtual GridBase *GaugeRedBlackGrid(void) =0; + //////////////////////////////////////////////////////////////// + // Split-communicator copies of this operator's gauge and fermion + // grids, returned in a new bundle with no Matrix or Linop yet. + // Partitions have MPI layout mpi_split. Collective. + //////////////////////////////////////////////////////////////// + SplitOperator *MakeSplitGrids(const Coordinate &mpi_split) + { + GridCartesian *U = dynamic_cast(GaugeGrid()); + GRID_ASSERT(U != nullptr); + GRID_ASSERT(mpi_split.size() == Nd); + + double t0 = usecond(); + + SplitOperator *split = new SplitOperator(); + + split->GaugeGrid = new GridCartesian(U->FullDimensions(),U->_simd_layout,mpi_split,*U); + split->GaugeRBGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(split->GaugeGrid); + split->Partition = GridSplitVectorIndex(U,split->GaugeGrid); + split->Partitions = U->ProcessorCount()/split->GaugeGrid->ProcessorCount(); + + if ( FermionGrid()->Nd() == Nd+1 ) { + int Ls = FermionGrid()->_fdimensions[0]; + split->FermionGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,split->GaugeGrid); + split->FermionRBGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,split->GaugeGrid); + } else { + split->FermionGrid = split->GaugeGrid; + split->FermionRBGrid = split->GaugeRBGrid; + } + std::cout << GridLogPerformance << "MakeSplitGrids: " << (usecond()-t0)/1.0e6 << " s" << std::endl; + return split; + } + // override multiply virtual void M (const FermionField &in, FermionField &out)=0; virtual void Mdag (const FermionField &in, FermionField &out)=0; @@ -88,6 +120,10 @@ public: virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; + // 1 if MoeDeriv/MeoDeriv return a force on a single 4D checkerboard (4D red-black), + // 0 if they return it on the full 4D gauge grid (a checkerboard that includes s) + virtual int CheckerboardedForce(void) { return 1; }; + virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac virtual void MdirAll(const FermionField &in, std::vector &out)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac diff --git a/Grid/qcd/action/fermion/MobiusFermion.h b/Grid/qcd/action/fermion/MobiusFermion.h index cb3c85fd6..e1c3c5141 100644 --- a/Grid/qcd/action/fermion/MobiusFermion.h +++ b/Grid/qcd/action/fermion/MobiusFermion.h @@ -70,6 +70,37 @@ public: } + // Exact type only: derived operators (e.g. ScaledShamirFermion) must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(MobiusFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + MobiusFermion *clone = new MobiusFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + this->_b, + this->_c, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h index 562f08f4a..ecad7f2c9 100644 --- a/Grid/qcd/action/fermion/WilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h @@ -91,6 +91,60 @@ public: // Derivative parts unpreconditioned pseudofermions void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag); + // Install already computed clover term and inverse, and derive the checkerboarded parts + void ImportCloverTerms(const CloverField &C, const CloverField &Cinv) + { + CloverTerm = C; + CloverTermInv = Cinv; + pickCheckerboard(Even, CloverTermEven, CloverTerm); + pickCheckerboard(Odd, CloverTermOdd, CloverTerm); + pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); + pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm)); + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); + pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv); + pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); + pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); + } + + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(WilsonCloverFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field and clover terms are overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + WilsonCloverFermion *clone = + new WilsonCloverFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass, + csw_r, + csw_t, + this->anisotropyCoeff, + this->Params); + // The stored coefficients are post-scaled; the constructor would scale them again + clone->csw_r = csw_r; + clone->csw_t = csw_t; + + CloverField sC(split->GaugeGrid); + CloverField sCinv(split->GaugeGrid); + Grid_split(CloverTerm,sC); + Grid_split(CloverTermInv,sCinv); + clone->ImportCloverTerms(sC,sCinv); + + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + public: // here fixing the 4 dimensions, make it more general? diff --git a/Grid/qcd/action/fermion/WilsonFermion.h b/Grid/qcd/action/fermion/WilsonFermion.h index 41a2471e2..0a2bb570e 100644 --- a/Grid/qcd/action/fermion/WilsonFermion.h +++ b/Grid/qcd/action/fermion/WilsonFermion.h @@ -147,6 +147,48 @@ public: // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); + // Install an already doubled gauge field (phases, twists, anisotropy and -1/2 applied) + void ImportDoubledGauge(const DoubledGaugeField &U) + { + Umu = U; + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd ,Umu); + } + + // Redistribute this operator's doubled gauge field onto clone's split grid. Collective. + void SplitDoubledGaugeInto(WilsonFermion &clone) + { + DoubledGaugeField sUmu(clone.GaugeGrid()); + Grid_split(Umu,sUmu); + clone.ImportDoubledGauge(sUmu); + } + + // Exact type only: derived operators (e.g. WilsonTMFermion) must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(WilsonFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + WilsonFermion *clone = new WilsonFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + mass, + this->Params, + anisotropyCoeff); + SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h index ea0eaa32e..4d5df4fbb 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonFermion5D.h @@ -181,6 +181,27 @@ public: // DoubleStore void ImportGauge(const GaugeField &_Umu); + // Install an already doubled gauge field (phases, twists and -1/2 applied) + void ImportDoubledGauge(const DoubledGaugeField &U) + { + Umu = U; + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd ,Umu); + } + + // Redistribute this operator's doubled gauge field onto clone's split grid. Collective. + void SplitDoubledGaugeInto(WilsonFermion5D &clone) + { + double t0 = usecond(); + DoubledGaugeField sUmu(clone.GaugeGrid()); + Grid_split(Umu,sUmu); + double t1 = usecond(); + clone.ImportDoubledGauge(sUmu); + double t2 = usecond(); + std::cout << GridLogPerformance << "SplitDoubledGaugeInto: Grid_split " << (t1-t0)/1.0e6 + << " s import " << (t2-t1)/1.0e6 << " s" << std::endl; + } + /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/ZMobiusFermion.h b/Grid/qcd/action/fermion/ZMobiusFermion.h index de832b7d7..2fa1b8d8f 100644 --- a/Grid/qcd/action/fermion/ZMobiusFermion.h +++ b/Grid/qcd/action/fermion/ZMobiusFermion.h @@ -67,6 +67,44 @@ public: this->SetCoefficientsInternal(1.0,zgamma,b,c); } + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(ZMobiusFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + // Constructor argument only; CloneCoefficientsInto restores _gamma at full precision + std::vector gamma(this->Ls); + for ( int s=0; sLs; s++ ) { + gamma[s] = this->_gamma[s]; + } + + ZMobiusFermion *clone = new ZMobiusFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + gamma, + this->_b, + this->_c, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h index 29408ea29..38e419bd8 100644 --- a/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h +++ b/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -64,6 +64,20 @@ public: GRID_ASSERT(U.Checkerboard()==Odd); GRID_ASSERT(V.Checkerboard()==U.Checkerboard()); + // Forces returned on the full 4D grid (checkerboard including s): same sequence, summed + if ( !this->_Mat.CheckerboardedForce() ) { + GaugeField ForceO(Force.Grid()); + GaugeField ForceE(Force.Grid()); + this->_Mat.Meooe (V,tmp1); + this->_Mat.MooeeInv(tmp1,tmp2); + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); + this->_Mat.MeooeDag (U,tmp1); + this->_Mat.MooeeInvDag(tmp1,tmp2); + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); + Force = -(ForceO+ForceE); + return; + } + // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // it is not conformable with the HMC force field // Case: Ls vectorised fields @@ -109,6 +123,20 @@ public: GRID_ASSERT(V.Checkerboard()==Odd); GRID_ASSERT(V.Checkerboard()==V.Checkerboard()); + // Forces returned on the full 4D grid (checkerboard including s): same sequence, summed + if ( !this->_Mat.CheckerboardedForce() ) { + GaugeField ForceO(Force.Grid()); + GaugeField ForceE(Force.Grid()); + this->_Mat.MeooeDag (V,tmp1); + this->_Mat.MooeeInvDag(tmp1,tmp2); + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); + this->_Mat.Meooe (U,tmp1); + this->_Mat.MooeeInv(tmp1,tmp2); + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); + Force = -(ForceO+ForceE); + return; + } + // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // it is not conformable with the HMC force field // INHERIT FROM THE Force field instead diff --git a/Grid/qcd/smearing/GaugeConfiguration.h b/Grid/qcd/smearing/GaugeConfiguration.h index 6f6fd016e..cc7fa8ec8 100644 --- a/Grid/qcd/smearing/GaugeConfiguration.h +++ b/Grid/qcd/smearing/GaugeConfiguration.h @@ -125,6 +125,16 @@ public: return SmearedSet[Level]; } + /*! @brief Level n of the smearing chain: 0 is the thin links, n is n smearing steps */ + const GaugeField &get_level(int n) const + { + GRID_ASSERT(n >= 0 && n <= (int)smearingLevels); + if ( n == 0 ) { + return *ThinLinks; + } + return SmearedSet[n-1]; + } + //==================================================================== void set_iLambda(GaugeLinkField& iLambda, GaugeLinkField& e_iQ, const GaugeLinkField& iQ, const GaugeLinkField& Sigmap, @@ -277,6 +287,45 @@ public: SigmaTilde=Gimpl::projectForce(SigmaTilde); // Ta } + + //==================================================================== + // Force on the thin links from an action that depends on every level of the chain. + // LevelForce[n], n = 0..smearingLevels, is the force on level n in the same + // "U times dS/dU" form smeared_force takes (zero fields allowed). The chain rule is + // applied from the top level down, adding each level's own force on the way: + // Sigma_L = U_L^dag F_L + // Sigma_n = AnalyticSmearedForce(Sigma_{n+1}, U_n) + U_n^dag F_n + // ThinForce = Ta( U_0 Sigma_0 ) + // With force on the top level only this is smeared_force. + // SmearedConfigurationMasked would need its own: its chain rule is level dependent. + virtual void smeared_force_levels(const std::vector &LevelForce, GaugeField &ThinForce) + { + GRID_ASSERT(LevelForce.size() == smearingLevels+1); + + GridBase *grid = LevelForce[0].Grid(); + GaugeField Sigma(grid); + GaugeLinkField tmp_mu(grid); + + // Sigma = U_L^dag F_L + for (int mu = 0; mu < Nd; mu++) { + tmp_mu = adj(peekLorentz(get_level(smearingLevels), mu)) * peekLorentz(LevelForce[smearingLevels], mu); + pokeLorentz(Sigma, tmp_mu, mu); + } + + for (int n = (int)smearingLevels - 1; n >= 0; n--) { + Sigma = AnalyticSmearedForce(Sigma, get_level(n)); + for (int mu = 0; mu < Nd; mu++) { + tmp_mu = peekLorentz(Sigma, mu) + adj(peekLorentz(get_level(n), mu)) * peekLorentz(LevelForce[n], mu); + pokeLorentz(Sigma, tmp_mu, mu); + } + } + + for (int mu = 0; mu < Nd; mu++) { + tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(Sigma, mu); + pokeLorentz(ThinForce, tmp_mu, mu); + } + ThinForce = Gimpl::projectForce(ThinForce); // Ta + } //==================================================================== virtual GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; } diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 0c5787879..e6f621cde 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -95,6 +95,8 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////// static Coordinate Grid_default_latt; static Coordinate Grid_default_mpi; +static Coordinate Grid_default_batched_solver_split; +static bool Grid_default_batched_solver_split_node = false; int GridThread::_threads =1; int GridThread::_hyperthreads=1; @@ -108,6 +110,8 @@ char *GridHostname(void) } const Coordinate &GridDefaultLatt(void) {return Grid_default_latt;}; const Coordinate &GridDefaultMpi(void) {return Grid_default_mpi;}; +const Coordinate &GridDefaultBatchedSolverSplit(void) {return Grid_default_batched_solver_split;}; +bool GridDefaultBatchedSolverSplitNode(void) {return Grid_default_batched_solver_split_node;}; const Coordinate GridDefaultSimd(int dims,int nsimd) { Coordinate layout(dims); @@ -562,6 +566,7 @@ void Grid_init(int *argc,char ***argv) std::cout< ~45 GB, hence --device-mem 40000 as in pvdagm_mixed_precision.job. +# Host memory: the device is an inclusive cache, so host RSS per rank is the whole +# ~45 GB working set (plus transient Grid_split staging): ~400 GB of a node's 512 GB. +# The driver prints MEMORY lines (host RSS, peak, allocator caches) at each phase. +# Grid_split stages 2 x 4.8 GB per rank on the host per transfer. +############################################################################## + +LUSTRE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle +RUNDIR=$LUSTRE/runs/${SLURM_JOB_NAME}_${SLURM_JOB_ID} +mkdir -p $RUNDIR && cd $RUNDIR && echo "RUNDIR $RUNDIR" + +cat << EOF > select_gpu +#!/bin/bash +export GPU_MAP=(0 1 2 3 7 6 5 4) +export NUMA_MAP=(3 3 1 1 2 2 0 0) +export GPU=\${GPU_MAP[\$SLURM_LOCALID]} +export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]} +export HIP_VISIBLE_DEVICES=\$GPU +unset ROCR_VISIBLE_DEVICES +if [ \$SLURM_PROCID = "0" ]; then echo \$*; fi +exec numactl -m \$NUMA -N \$NUMA \$* +EOF +chmod +x ./select_gpu + +# Source tree and build of THIS branch (feature/splitGridCGBatch) on Frontier +root=$LUSTRE/SplitCG/Grid/systems/Frontier +source $root/sourceme-rocm7.2.sh + +export OMP_NUM_THREADS=7 +ulimit -c 0 # no 22 GB GPU core dumps +export FI_MR_CACHE_MONITOR=kdreg2 # site default; device-buffer MPI on Slingshot +export MPICH_GPU_SUPPORT_ENABLED=1 +export MPICH_SMP_SINGLE_COPY_MODE=CMA +export MPICH_OFI_NIC_POLICY=GPU +# No GRID_ALLOC_NCACHE_LARGE here. With --enable-unified=no every Lattice lives in host +# memory (the device is an inclusive cache), so a deep large-allocation cache keeps freed +# host blocks on top of the full working set: measured on the laptop at 64 entries it +# held about as much again as the working set, enough to OOM a 4-node run of this job. +module load libfabric + +BIN=$root/tests/solver/Test_split_mobius_batched +OPTS="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 40000 --comms-overlap" +vol=64.64.64.128 +MPI_GEOM=2.2.2.4 + +# A thermalised 64^3x128 NERSC configuration. Leave empty for a hot start: +# fine for plumbing, but mass 0.026 on a hot field is not a physical +# iteration count and may run long. +CONFIG= + +PHYS="--Ls 12 --mass 0.026 --M5 1.8 --b 1.5 --c 0.5 --nbatch 4 --tol 1e-8" +if [ -n "$CONFIG" ]; then PHYS="$PHYS --config $CONFIG"; fi + +export GRID_STDOUT_ROOT=$RUNDIR/split +srun -N4 -n32 --kill-on-bad-exit=1 ./select_gpu $BIN --mpi $MPI_GEOM --grid $vol $OPTS \ + $PHYS --batched-solver-split node \ + --debug-stdout --log Error,Warning,Message,Performance \ + > log.split 2>&1 +echo "exit $?" + +############################################################################## +# Readout. Device OOM appears only on stderr as "hipMalloc failed". +############################################################################## +f=$(grep -l "Memory access fault\|NO_TRANSLATION\|hipMalloc failed\|out of memory\|illegal memory access\|GRID_ASSERT" $GRID_STDOUT_ROOT/*/Grid.stderr.* 2>/dev/null | head -1) +if [ -n "$f" ]; then + echo "FAULT in $f"; tail -20 $f | cut -c1-160 +fi +grep -h "MEMORY\|BatchedSolverSplit\|Inner CG iterations\|Total time\|true residual\|SUMMARY" \ + $GRID_STDOUT_ROOT/0/Grid.stdout.0 2>/dev/null | cut -c1-200 + +############################################################################## +# What to read. +# +# BatchedSolverSplit line: the partition layout must equal the node layout +# (ShmGrid), e.g. [1 2 2 2] of [2 2 2 4] : 4 partitions. A "straddle node +# boundaries" warning means the rank reordering did not give node-aligned +# blocks; the result is still correct but the inner CG talks off node. +# +# UNSPLIT vs SPLIT: per-rhs inner iterations should agree closely (the same +# solves, different reduction order); the SUMMARY lines give the wall clock. +# The gain is the strong-scaling loss of the 32-GCD inner CG against four +# 8-GCD CGs with node-local halos. "Split setup and transfer" in the Total +# time line prices the clone (grids, gauge split) and the Grid_split / +# Grid_unsplit traffic, which is host-staged. +############################################################################## diff --git a/tests/forces/ForceTest.h b/tests/forces/ForceTest.h new file mode 100644 index 000000000..9531358f0 --- /dev/null +++ b/tests/forces/ForceTest.h @@ -0,0 +1,138 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/forces/ForceTest.h + + Copyright (C) 2022 + +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 */ +#pragma once + +///////////////////////////////////////////////////////////////////////////////////////////// +// Finite-difference check of an action's force: refresh, S1 = S(U), step the links by eps +// along a random momentum P, take the force at the midpoint, step again, S2 = S(U''). +// The force predicts S2 - S1 to O(eps^3). +// +// Runs on a copy of U, so the caller's field is unchanged. Returns S2 - S1 - dSpred for +// the caller to assert on; halving eps should reduce it by about 8. +///////////////////////////////////////////////////////////////////////////////////////////// +template +RealD ForceTest(Action &action,LatticeGaugeField & Uin,MomentumFilterBase &Filter,RealD eps=0.005) +{ + GridBase *UGrid = Uin.Grid(); + + std::vector seeds({1,2,3,5}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + + LatticeColourMatrix Pmu(UGrid); + LatticeGaugeField P(UGrid); + LatticeGaugeField UdSdU(UGrid); + LatticeGaugeField U(UGrid); + U = Uin; + + std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); + Pmu= PeekIndex(P,mu); + dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; + } + ComplexD dSpred = sum(dS); + RealD diff = S2-S1-dSpred.real(); + + std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"< +RealD ForceTestScaling(Action &action,LatticeGaugeField &U,MomentumFilterBase &Filter,const std::string &name, + RealD eps0=0.01,int neps=4) +{ + std::vector eps(neps); + std::vector d(neps); + for(int i=0;i(action,U,Filter,eps[i]); + } + for(int i=0;i 0 ) { + std::cout << " ratio to previous " << d[i-1]/d[i]; + } + std::cout << std::endl; + } + RealD scale = d[neps-2]/d[neps-1]; + std::cout << GridLogMessage << name << ": smallest-eps ratio " << scale << " (expect ~8)" << std::endl; + return scale; +} diff --git a/tests/forces/Test_bdy.cc b/tests/forces/Test_bdy.cc index d328b49fb..106a73f31 100644 --- a/tests/forces/Test_bdy.cc +++ b/tests/forces/Test_bdy.cc @@ -99,91 +99,12 @@ Author: Peter Boyle using namespace std; using namespace Grid; +#include "ForceTest.h" + typedef MobiusFermionD FermionAction; typedef WilsonImplD FimplD; typedef WilsonImplD FermionImplPolicy; -template -void ForceTest(Action &action,LatticeGaugeField & U,MomentumFilterBase &Filter) -{ - GridBase *UGrid = U.Grid(); - - std::vector seeds({1,2,3,5}); - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); - - LatticeColourMatrix Pmu(UGrid); - LatticeGaugeField P(UGrid); - LatticeGaugeField UdSdU(UGrid); - - std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); - Pmu= PeekIndex(P,mu); - dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; - } - ComplexD dSpred = sum(dS); - RealD diff = S2-S1-dSpred.real(); - - std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"< -void ForceTest(Action &action,LatticeGaugeField & U,MomentumFilterBase &Filter) -{ - GridBase *UGrid = U.Grid(); - - std::vector seeds({1,2,3,5}); - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); - - LatticeColourMatrix Pmu(UGrid); - LatticeGaugeField P(UGrid); - LatticeGaugeField UdSdU(UGrid); - - std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); - Pmu= PeekIndex(P,mu); - dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; - } - ComplexD dSpred = sum(dS); - RealD diff = S2-S1-dSpred.real(); - - std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"< // composition identity [P M^-1 V Pdag][P V^-1 M Pdag] = 1 for the chosen // (P,Pdag) wall pair. NO PREDICTION is made about which convention // passes -- that is what the test decides. -// T2 ForceTest (idiom from Test_double_ratio.cc) : midpoint-derivative +// T2 ForceTest (tests/forces/ForceTest.h) : midpoint-derivative // check of deriv against S. Should PASS for BOTH conventions (S and // deriv use the same literal-adjoint pair by construction). // T3 Trivial-ratio control (V == M) : T1 with NumOp = DenOp. The solve @@ -124,69 +124,7 @@ RealD HeatbathIdentityTest(TwoFlavourRatio4DPseudoFermionAction &action, return defect; } -//////////////////////////////////////////////////////////////////// -// T2 : ForceTest idiom from Test_double_ratio.cc (midpoint derivative) -//////////////////////////////////////////////////////////////////// -template -void ForceTest(Action &action,LatticeGaugeField & U,MomentumFilterBase &Filter) -{ - GridBase *UGrid = U.Grid(); - - std::vector seeds({1,2,3,5}); - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); - - LatticeColourMatrix Pmu(UGrid); - LatticeGaugeField P(UGrid); - LatticeGaugeField UdSdU(UGrid); - - std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); - Pmu= PeekIndex(P,mu); - dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; - } - ComplexD dSpred = sum(dS); - RealD diff = S2-S1-dSpred.real(); - - std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"< // BOTH classes (E0a also validates the twin-eta capture). // E1 : S_classic == S_leftprec (relative, ~1e-8) // E2 : deriv_classic == deriv_leftprec (pointwise field norm, ~1e-8) -// F1 : ForceTest (Test_double_ratio.cc idiom) on the LeftPrec class. +// F1 : ForceTest (tests/forces/ForceTest.h) on the LeftPrec class. // // All asserts are hard: this is the regression gate for the new class. // Run small, e.g.: ./Test_dwf_ratio_leftprec --grid 8.8.8.8 @@ -138,67 +138,7 @@ public: } }; -//////////////////////////////////////////////////////////////////// -// ForceTest idiom from Test_double_ratio.cc (midpoint derivative) -//////////////////////////////////////////////////////////////////// -template -void ForceTest(Action &action,LatticeGaugeField & U,MomentumFilterBase &Filter) -{ - GridBase *UGrid = U.Grid(); - - std::vector seeds({1,2,3,5}); - GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); - - LatticeColourMatrix Pmu(UGrid); - LatticeGaugeField P(UGrid); - LatticeGaugeField UdSdU(UGrid); - - std::cout << GridLogMessage << "*********************************************************"<(UdSdU,mu); - Pmu= PeekIndex(P,mu); - dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0; - } - ComplexD dSpred = sum(dS); - RealD diff = S2-S1-dSpred.real(); - - std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<> mpi_split[k]; - } - break; + if ( GridCmdOptionExists(argv,argv+argc,"--split") ) { + std::vector split; + GridCmdOptionIntVector(GridCmdOptionPayload(argv,argv+argc,"--split"),split); + GRID_ASSERT(split.size() == mpi_layout.size()); + for(int k=0;k> mpi_split[k]; - } - break; + if ( GridCmdOptionExists(argv,argv+argc,"--split") ) { + std::vector split; + GridCmdOptionIntVector(GridCmdOptionPayload(argv,argv+argc,"--split"),split); + GRID_ASSERT(split.size() == mpi_layout.size()); + for(int k=0;k> mpi_split[k]; - } - break; + if ( GridCmdOptionExists(argv,argv+argc,"--split") ) { + std::vector split; + GridCmdOptionIntVector(GridCmdOptionPayload(argv,argv+argc,"--split"),split); + GRID_ASSERT(split.size() == mpi_layout.size()); + for(int k=0;k> mpi_split[k]; - } - break; + if ( GridCmdOptionExists(argv,argv+argc,"--split") ) { + std::vector split; + GridCmdOptionIntVector(GridCmdOptionPayload(argv,argv+argc,"--split"),split); + GRID_ASSERT(split.size() == mpi_layout.size()); + for(int k=0;k + + 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 */ +///////////////////////////////////////////////////////////////////////////////////////////// +// MixedPrecisionConjugateGradientBatched with inner solves on split-communicator +// partitions (--batched-solver-split), against the same solver unsplit. +// +// mpirun -n 2 ./Test_split_mixedprec_batched --grid 8.8.8.8 --mpi 1.1.1.2 --batched-solver-split 1.1.1.1 +// +// With no --batched-solver-split the partitions are single ranks (1.1.1.1). +// NBATCH is deliberately not a multiple of the partition count, so the last group is +// zero-padded. +///////////////////////////////////////////////////////////////////////////////////////////// +#include + +using namespace std; +using namespace Grid; + +const int NBATCH = 3; +const RealD TOLERANCE = 1.0e-8; + +template +RealD RelativeDifference(const Field &a,const Field &b) +{ + Field diff(a.Grid()); + diff = a - b; + return std::sqrt(norm2(diff)/norm2(b)); +} + +///////////////////////////////////////////////////////////////////////////////////////////// +// The solver assumes vector p of a group lands in the partition whose +// GridSplitVectorIndex is p: fill vector p with the constant p+1, split, and check +// every partition sees its own index. +///////////////////////////////////////////////////////////////////////////////////////////// +void CheckPartitionOrder(GridCartesian *UGrid,const Coordinate &layout) +{ + GridCartesian SGrid(UGrid->FullDimensions(),UGrid->_simd_layout,layout,*UGrid); + int P = UGrid->ProcessorCount()/SGrid.ProcessorCount(); + int index = GridSplitVectorIndex(UGrid,&SGrid); + + std::vector full(P,UGrid); + for(int p=0;p +void CheckSplitBatched(const std::string &name, + LinearOperatorBase &Linop_d, + LinearOperatorBase &Linop_f, + GridBase *grid_f, + std::vector &src, + const Coordinate &layout) +{ + std::cout << GridLogMessage << "==================================================" << std::endl; + std::cout << GridLogMessage << name << std::endl; + std::cout << GridLogMessage << "==================================================" << std::endl; + + int cb = src[0].Checkerboard(); + + SplitOperator *split = Linop_f.SplitClone(layout); + GRID_ASSERT(split != nullptr); + int P = split->Partitions; + + std::vector full(P,grid_f); + std::vector back(P,grid_f); + std::vector Mfull(P,grid_f); + for(int p=0;pFieldGrid); + FieldF s_out(split->FieldGrid); + Grid_split(full,s_in); + Grid_unsplit(back,s_in); + for(int p=0;pLinop->HermOp(s_in,s_out); + Grid_unsplit(back,s_out); + for(int p=0;p sol_ref(NBatch,src[0].Grid()); + std::vector sol_split(NBatch,src[0].Grid()); + for(int i=0;i mCG(TOLERANCE,10000,50,1000,grid_f,Linop_f,Linop_d, + true,Coordinate(),false); + mCG(src,sol_ref); + } + { + std::cout << GridLogMessage << name << ": split batched solve, partition layout " << layout << std::endl; + MixedPrecisionConjugateGradientBatched mCG(TOLERANCE,10000,50,1000,grid_f,Linop_f,Linop_d, + true,layout,false); + GRID_ASSERT( (mCG.Partitions == 1) || (mCG.split != nullptr) ); + mCG(src,sol_split); + } + + FieldD Msol(src[0].Grid()); + Msol.Checkerboard() = cb; + for(int i=0;i seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG4(UGrid_d); + GridParallelRNG RNG5(FGrid_d); + RNG4.SeedFixedIntegers(seeds4); + RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeFieldD Umu_d(UGrid_d); + LatticeGaugeFieldF Umu_f(UGrid_f); + SU::HotConfiguration(RNG4,Umu_d); + precisionChange(Umu_f,Umu_d); + + // Antiperiodic in time for some actions, so boundary phases must travel with the links + WilsonImplParams antiperiodic; + antiperiodic.boundary_phases[Nd-1] = -1.0; + + // Sources: odd checkerboard (Schur) and full lattice (MdagM), 4d and 5d + std::vector src4_o(NBATCH,UrbGrid_d); + std::vector src5_o(NBATCH,FrbGrid_d); + std::vector src5(NBATCH,FGrid_d); + LatticeFermionD tmp4(UGrid_d); + LatticeFermionD tmp5(FGrid_d); + for(int i=0;i Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("Wilson SchurDiagMooee",Ld,Lf,UrbGrid_f,src4_o,layout); + } + + ////////////////////////////////////////// + // Wilson clover + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + WilsonCloverFermionD Dd(Umu_d,*UGrid_d,*UrbGrid_d,mass,csw_r,csw_t); + WilsonCloverFermionF Df(Umu_f,*UGrid_f,*UrbGrid_f,mass,csw_r,csw_t); + SchurDiagMooeeOperator Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("WilsonClover SchurDiagMooee",Ld,Lf,UrbGrid_f,src4_o,layout); + } + + ////////////////////////////////////////// + // Compact Wilson clover, antiperiodic + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + RealD cF = 1.0; + WilsonAnisotropyCoefficients anis; + CompactWilsonCloverFermionD Dd(Umu_d,*UGrid_d,*UrbGrid_d,mass,csw_r,csw_t,cF,anis,antiperiodic); + CompactWilsonCloverFermionF Df(Umu_f,*UGrid_f,*UrbGrid_f,mass,csw_r,csw_t,cF,anis,antiperiodic); + SchurDiagMooeeOperator Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("CompactWilsonClover SchurDiagMooee antiperiodic",Ld,Lf,UrbGrid_f,src4_o,layout); + } + + ////////////////////////////////////////// + // Domain wall + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD M5 = 1.8; + DomainWallFermionD Dd(Umu_d,*FGrid_d,*FrbGrid_d,*UGrid_d,*UrbGrid_d,mass,M5); + DomainWallFermionF Df(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5); + SchurDiagMooeeOperator Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("DomainWall SchurDiagMooee",Ld,Lf,FrbGrid_f,src5_o,layout); + } + + ////////////////////////////////////////// + // Mobius, antiperiodic, with unequal + // masses; all wrapper kinds + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD M5 = 1.8; + RealD b = 1.5; + RealD c = 0.5; + MobiusFermionD Dd(Umu_d,*FGrid_d,*FrbGrid_d,*UGrid_d,*UrbGrid_d,mass,M5,b,c,antiperiodic); + MobiusFermionF Df(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5,b,c,antiperiodic); + Dd.SetMass(0.1,0.12); + Df.SetMass(0.1,0.12); + { + SchurDiagMooeeOperator Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("Mobius SchurDiagMooee antiperiodic",Ld,Lf,FrbGrid_f,src5_o,layout); + } + { + SchurDiagOneOperator Ld(Dd); + SchurDiagOneOperator Lf(Df); + CheckSplitBatched("Mobius SchurDiagOne antiperiodic",Ld,Lf,FrbGrid_f,src5_o,layout); + } + { + SchurDiagTwoOperator Ld(Dd); + SchurDiagTwoOperator Lf(Df); + CheckSplitBatched("Mobius SchurDiagTwo antiperiodic",Ld,Lf,FrbGrid_f,src5_o,layout); + } + { + MdagMLinearOperator Ld(Dd); + MdagMLinearOperator Lf(Df); + CheckSplitBatched("Mobius MdagM antiperiodic",Ld,Lf,FGrid_f,src5,layout); + } + } + + ////////////////////////////////////////// + // ZMobius + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD M5 = 1.8; + RealD b = 1.0; + RealD c = 0.0; + std::vector gamma(Ls); + for(int s=0;s Ld(Dd); + SchurDiagMooeeOperator Lf(Df); + CheckSplitBatched("ZMobius SchurDiagMooee",Ld,Lf,FrbGrid_f,src5_o,layout); + } + + ////////////////////////////////////////// + // A derived operator must not inherit + // its parent's SplitClone + ////////////////////////////////////////// + { + RealD mass = 0.1; + RealD mu = 0.1; + WilsonTMFermionF Df(Umu_f,*UGrid_f,*UrbGrid_f,mass,mu); + SchurDiagMooeeOperator Lf(Df); + SplitOperator *split = Lf.SplitClone(layout); + std::cout << GridLogMessage << "WilsonTM SplitClone refused: " << (split == nullptr) << std::endl; + GRID_ASSERT(split == nullptr); + } + + std::cout << GridLogMessage << "Test_split_mixedprec_batched: all checks passed" << std::endl; + + Grid_finalize(); +} diff --git a/tests/solver/Test_split_mobius_batched.cc b/tests/solver/Test_split_mobius_batched.cc new file mode 100644 index 000000000..6f3c07b76 --- /dev/null +++ b/tests/solver/Test_split_mobius_batched.cc @@ -0,0 +1,236 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/solver/Test_split_mobius_batched.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 */ +///////////////////////////////////////////////////////////////////////////////////////////// +// Production-size timing of MixedPrecisionConjugateGradientBatched for Mobius with the +// Hadrons default SchurDiagMooeeOperator: the same batch solved without and then with +// split inner solves (--batched-solver-split), reporting wall clock, per-rhs iterations +// and true residuals for each. +// +// --Ls 12 --mass 0.026 --M5 1.8 --b 1.5 --c 0.5 --nbatch 4 --tol 1e-8 +// --config (omit for a hot configuration: timing only, not physics) +// --nounsplit (skip the reference unsplit solve) +// --repeat N (split solve N times; host RSS must not grow between them) +// +// MEMORY lines report host RSS (current and peak) and allocator cache sizes, maximum over +// ranks, at each phase: with --enable-unified=no every Lattice lives in host memory. +// +// Only one solution vector is kept, so the driver's own footprint is two batches of +// double red-black 5d fields; the solver adds about as much again. +///////////////////////////////////////////////////////////////////////////////////////////// +#include +#include +#ifdef __APPLE__ +#include +#endif + +using namespace std; +using namespace Grid; + +// Host memory of this process in GB: current resident set and its high-water mark +void HostRSS(RealD ¤t,RealD &peak) +{ + struct rusage ru; + getrusage(RUSAGE_SELF,&ru); +#ifdef __APPLE__ + peak = ru.ru_maxrss/1.0e9; // bytes on macOS + mach_task_basic_info_data_t info; + mach_msg_type_number_t count = MACH_TASK_BASIC_INFO_COUNT; + task_info(mach_task_self(),MACH_TASK_BASIC_INFO,(task_info_t)&info,&count); + current = info.resident_size/1.0e9; +#else + peak = ru.ru_maxrss*1024.0/1.0e9; // kilobytes on Linux + long pages = 0; + long resident = 0; + FILE *f = fopen("/proc/self/statm","r"); + if ( f ) { + if ( fscanf(f,"%ld %ld",&pages,&resident) != 2 ) { + resident = 0; + } + fclose(f); + } + current = resident*(RealD)sysconf(_SC_PAGESIZE)/1.0e9; +#endif +} + +// Largest values over ranks: host RSS now and at peak, and the allocator caches +void ReportMemory(GridBase *grid,const std::string &phase) +{ + RealD rss; + RealD peak; + HostRSS(rss,peak); + RealD hostcache = MemoryManager::HostCacheBytes()/1.0e9; + RealD devcache = MemoryManager::DeviceCacheBytes()/1.0e9; + grid->GlobalMax(rss); + grid->GlobalMax(peak); + grid->GlobalMax(hostcache); + grid->GlobalMax(devcache); + std::cout << GridLogMessage << "MEMORY " << phase + << " : host RSS " << rss << " GB, peak " << peak + << " GB; allocator cache host " << hostcache << " GB, device " << devcache + << " GB (max over ranks)" << std::endl; + HostMemoryReport(grid,GridLogMessage,phase); +} + +typedef LatticeFermionD FieldD; +typedef LatticeFermionF FieldF; + +template +T CmdOption(int argc,char **argv,const std::string &name,T def) +{ + T val = def; + if ( GridCmdOptionExists(argv,argv+argc,name) ) { + std::stringstream ss(GridCmdOptionPayload(argv,argv+argc,name)); + ss >> val; + } + return val; +} + +void SolveAndReport(const std::string &label, + MixedPrecisionConjugateGradientBatched &mCG, + LinearOperatorBase &Linop_d, + std::vector &src, + std::vector &sol) +{ + int nbatch = src.size(); + for(int i=0;i (argc,argv,"--Ls",12); + RealD mass = CmdOption (argc,argv,"--mass",0.026); + RealD M5 = CmdOption (argc,argv,"--M5",1.8); + RealD b = CmdOption (argc,argv,"--b",1.5); + RealD c = CmdOption (argc,argv,"--c",0.5); + int nbatch = CmdOption (argc,argv,"--nbatch",4); + RealD tol = CmdOption (argc,argv,"--tol",1.0e-8); + std::string config = CmdOption(argc,argv,"--config",std::string("")); + bool unsplit = !GridCmdOptionExists(argv,argv+argc,"--nounsplit"); + int repeat = CmdOption (argc,argv,"--repeat",1); + + std::cout << GridLogMessage << "Mobius Ls " << Ls << " mass " << mass << " M5 " << M5 + << " b " << b << " c " << c << " nbatch " << nbatch << " tol " << tol << std::endl; + + GridCartesian *UGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_d); + GridCartesian *FGrid_d = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_d); + GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_d); + + GridCartesian *UGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f); + GridCartesian *FGrid_f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f); + GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG4(UGrid_d); + GridParallelRNG RNG5(FGrid_d); + RNG4.SeedFixedIntegers(seeds4); + RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeFieldD Umu_d(UGrid_d); + LatticeGaugeFieldF Umu_f(UGrid_f); + if ( config.size() ) { + FieldMetaData header; + NerscIO::readConfiguration(Umu_d,header,config); + } else { + std::cout << GridLogMessage << "No --config: hot configuration, timing only" << std::endl; + SU::HotConfiguration(RNG4,Umu_d); + } + precisionChange(Umu_f,Umu_d); + ReportMemory(UGrid_d,"gauge field ready"); + + // Antiperiodic in time, as in production + WilsonImplParams params; + params.boundary_phases[Nd-1] = -1.0; + + MobiusFermionD Dd(Umu_d,*FGrid_d,*FrbGrid_d,*UGrid_d,*UrbGrid_d,mass,M5,b,c,params); + MobiusFermionF Df(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5,b,c,params); + SchurDiagMooeeOperator Linop_d(Dd); + SchurDiagMooeeOperator Linop_f(Df); + + std::vector src(nbatch,FrbGrid_d); + std::vector sol(nbatch,FrbGrid_d); + { + FieldD tmp(FGrid_d); + for(int i=0;i mCG(tol,10000,50,10000,FrbGrid_f,Linop_f,Linop_d, + true,Coordinate(),false); + SolveAndReport("UNSPLIT",mCG,Linop_d,src,sol); + ReportMemory(UGrid_d,"after unsplit solve"); + } + + if ( repeat > 0 ) { + // Split layout from the command line; the split operator is cloned here, once + MixedPrecisionConjugateGradientBatched mCG(tol,10000,50,10000,FrbGrid_f,Linop_f,Linop_d); + ReportMemory(UGrid_d,"after split solver construction"); + // Repeated split solves reuse the clone and expose allocations not released between calls + for(int r=0;r