diff --git a/CLAUDE.md b/CLAUDE.md index 741262340..ebd2bffae 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: `CoarsenedMatrix.h` (coarse operator), `GeneralCoarsenedMatrix.h` and `GeneralCoarsenedMatrixMultiRHS.h` (general coarsening supporting multi-RHS solves), `Aggregates.h` (near-null vector construction), `Geometry.h` (coarse-grid geometry). `MultiGrid.h` is the top-level include. +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. ### GPU acceleration diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index d99c6dd54..575a03a30 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -214,6 +214,29 @@ public: } }; +//////////////////////////////////////////////////////////////////// +// Present any LinearOperatorBase as Hermitian: Op = AdjOp = HermOp. +// For coarsening (CoarsenOperator applies Op) an HPD operator whose +// Op is a factor and HermOp the product, e.g. SchurDiagMooee or MdagM. +//////////////////////////////////////////////////////////////////// +template +class HermOpAdaptor : public LinearOperatorBase { + LinearOperatorBase &_Mat; +public: + HermOpAdaptor(LinearOperatorBase &Mat): _Mat(Mat){}; + 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){ _Mat.HermOp(in,out); } + void AdjOp (const Field &in, Field &out){ _Mat.HermOp(in,out); } + void HermOp (const Field &in, Field &out){ _Mat.HermOp(in,out); } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + HermOp(in,out); + ComplexD dot = innerProduct(in,out); + n1=real(dot); + n2=norm2(out); + } +}; //////////////////////////////////////////////////////////////////// // Wrap an already herm matrix diff --git a/Grid/algorithms/blas/BatchedBlas.h b/Grid/algorithms/blas/BatchedBlas.h index a4de5d696..4f245c64e 100644 --- a/Grid/algorithms/blas/BatchedBlas.h +++ b/Grid/algorithms/blas/BatchedBlas.h @@ -1942,6 +1942,10 @@ public: Cs); synchronise(); + // Synchronise ONCE, after the whole train of calls. A synchronise inside + // the loop measures launch and completion latency, which at small shapes + // exceeds the kernel: at M=K=60, N=1, batch 1024 it hid a factor of two + // between the precisions and reported them level. RealD t0 = usecond(); for(int i=0;i ipiv((uint64_t)N); + deviceVector info(1); + double t0 = usecond(); + auto st1 = rocsolver_cgetrf_64(handle, N, N, (rocblas_float_complex *)A, N, &ipiv[0], &info[0]); + GRID_ASSERT(st1 == rocblas_status_success); + accelerator_barrier(); + lastGetrfUs = usecond()-t0; + int64_t info_h = -1; acceleratorCopyFromDevice(&info[0], &info_h, sizeof(int64_t)); + GRID_ASSERT(info_h == 0); + deviceVector X((uint64_t)N*N); + { ComplexF *x = &X[0]; const int64_t NN = N; + accelerator_for(idx, (uint64_t)N*N, 1, { int64_t j = idx/NN, i = idx - j*NN; x[idx] = (i==j) ? ComplexF(1.0,0.0) : ComplexF(0.0,0.0); }); + accelerator_barrier(); } + auto st2 = rocsolver_cgetrs_64(handle, rocblas_operation_none, N, N, + (rocblas_float_complex *)A, N, &ipiv[0], + (rocblas_float_complex *)&X[0], N); + GRID_ASSERT(st2 == rocblas_status_success); + accelerator_barrier(); + lastGetrsUs = usecond()-t0-lastGetrfUs; + acceleratorCopyDeviceToDevice((void *)&X[0], (void *)A, (uint64_t)N*N*sizeof(ComplexF)); +#else + deviceVector bp(1); std::vector ptr(1); ptr[0] = A; + acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexF*)); + inverseBatched(N, bp); +#endif + } void inverseBatched(int64_t N, deviceVector &Amat) { diff --git a/Grid/algorithms/deflation/MultiRHSBlockProject.h b/Grid/algorithms/deflation/MultiRHSBlockProject.h index 1afceba80..a264ff0cf 100644 --- a/Grid/algorithms/deflation/MultiRHSBlockProject.h +++ b/Grid/algorithms/deflation/MultiRHSBlockProject.h @@ -93,6 +93,25 @@ public: deviceVector BLAS_F; // nrhs x fine_vol * words -- the sources deviceVector BLAS_C; // nrhs x coarse_vol * nbasis -- the coarse coeffs + // Resident device memory this object holds. deviceVector is not + // evictable, so this counts against the hard budget, not the cache. + // Capacity, not size: a call with fewer right-hand sides shrinks the size + // and frees nothing. + uint64_t DeviceBytes(void) + { + return (uint64_t)(BLAS_V.capacity()+BLAS_F.capacity()+BLAS_C.capacity())*sizeof(scalar); + } + // The import/export scratch grows to the largest right-hand-side count + // any call has used and keeps that capacity. Setup projects the whole + // basis in one call, so unless it is released the scratch stays as large + // as the basis store for the rest of the run. The next call regrows it + // to the size that call needs. + void ReleaseScratch(void) + { + deviceVector().swap(BLAS_F); + deviceVector().swap(BLAS_C); + } + RealD blasNorm2(deviceVector &blas) { scalar ss(0.0); @@ -117,9 +136,14 @@ public: block_vol=0; coarse_vol=0; words=0; - BLAS_V.resize(0); - BLAS_F.resize(0); - BLAS_C.resize(0); + // deviceVector is a std::vector with a device allocator, so resize(0) + // drops the SIZE and keeps the CAPACITY: it returns no device memory at + // all. The basis store is nbasis fine vectors (10.9 GB at nbasis 64 on + // a 48^3x96 Ls=24 rank), so that has to be a real free. Swap with an + // empty vector, which destroys the buffer. + deviceVector().swap(BLAS_V); + deviceVector().swap(BLAS_F); + deviceVector().swap(BLAS_C); } void Allocate(int _nbasis,GridBase *_fgrid,GridBase *_cgrid) { @@ -153,25 +177,39 @@ public: BLAS_V.resize (fine_vol * words * nbasis ); } - void ImportFineGridVectors(std::vector &vecs, deviceVector &blas) + //////////////////////////////////////////////////////////////////////////// + // The fine side is templated on the INCOMING vector object, not fixed to + // Field. Field fixes only the STORAGE: the BLAS scalar, the word count and + // the full local geometry. Import is already a layout transformation -- + // (oSite,lane) is unpacked to full local coordinates and regathered in + // block order -- so a precision change costs nothing here, and neither does + // a different SIMD layout. That is what lets one projector serve an fp64 + // outer level and an fp32 coarse sector with no second basis store. + // + // The full LOCAL dimensions must agree; the SIMD layout and the precision + // need not. + //////////////////////////////////////////////////////////////////////////// + template + void ImportFineGridVectors(std::vector > &vecs, deviceVector &blas) { GRID_TRACE("ImportFineGridVectors"); int nvec = vecs.size(); - typedef typename Field::vector_object vobj; - // std::cout << GridLogMessage <<" BlockProjector importing "<_ndimension; + for(int d=0;d<_ndimension;d++) GRID_ASSERT(fgrid->_ldimensions[d] == fine_ldimensions[d]); + GRID_ASSERT(sizeof(fine_scalar_object)/sizeof(fine_scalar) == words); uint64_t sz = blas.size(); acceleratorMemSet(&blas[0],0,blas.size()*sizeof(scalar)); - Coordinate fine_rdimensions = fine_grid->_rdimensions; + Coordinate fine_rdimensions = fgrid->_rdimensions; Coordinate coarse_l = coarse_ldimensions; Coordinate block_l = block_ldimensions; - Coordinate fsimd = fine_simd; + Coordinate fsimd = fgrid->_simd_layout; int64_t bv= block_vol; for(int v=0;voSites(); + int64_t osites = fgrid->oSites(); // loop over fine sites const int Nsimd = vobj::Nsimd(); @@ -213,7 +251,7 @@ public: Lexicographic::IndexFromCoor(coor_c,sc,coarse_l); Lexicographic::IndexFromCoor(coor_b,sb,block_l); - scalar_object data = extractLane(lane,fineData[sf]); + fine_scalar_object data = extractLane(lane,fineData[sf]); // BLAS_F[coarse_vol][nvec][block_vol][words] int64_t site = (sc*nvec + v)*bv @@ -221,9 +259,9 @@ public: // GRID_ASSERT(site*lwords &blas) + template + void ImportFineGridMrhsVectors(Lattice &vec_mrhs, deviceVector &blas) { - typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object fine_scalar_object; + typedef typename vobj::scalar_type fine_scalar; GridBase *fine_mrhs_grid = vec_mrhs.Grid(); int _ndimension = coarse_grid->_ndimension; @@ -255,27 +295,31 @@ public: GRID_ASSERT(fine_mrhs_grid->_ndimension == _ndimension+1); GRID_ASSERT(fine_mrhs_grid->_simd_layout[0] == 1); GRID_ASSERT(fine_mrhs_grid->_processors[0] == 1); + // Full LOCAL dimensions must agree; the SIMD layout and precision need not for(int d=0;d<_ndimension;d++){ - GRID_ASSERT(fine_mrhs_grid->_rdimensions[d+1] == fine_grid->_rdimensions[d]); - GRID_ASSERT(fine_mrhs_grid->_simd_layout[d+1] == fine_grid->_simd_layout[d]); + GRID_ASSERT(fine_mrhs_grid->_ldimensions[d+1] == fine_ldimensions[d]); } + GRID_ASSERT(sizeof(fine_scalar_object)/sizeof(fine_scalar) == words); int nvec = fine_mrhs_grid->_rdimensions[0]; // nrhs uint64_t sz = blas.size(); acceleratorMemSet(&blas[0],0,blas.size()*sizeof(scalar)); Coordinate fine_mrhs_rdimensions = fine_mrhs_grid->_rdimensions; - Coordinate fine_rdimensions = fine_grid->_rdimensions; + Coordinate fine_rdimensions(_ndimension); + Coordinate fsimd(_ndimension); + for(int d=0;d<_ndimension;d++){ + fine_rdimensions[d] = fine_mrhs_grid->_rdimensions[d+1]; + fsimd[d] = fine_mrhs_grid->_simd_layout[d+1]; + } Coordinate coarse_l = coarse_ldimensions; Coordinate block_l = block_ldimensions; - Coordinate fsimd = fine_simd; int64_t bv= block_vol; autoView( fineData , vec_mrhs, AcceleratorRead); auto blasData_p = &blas[0]; auto fineData_p = &fineData[0]; - int64_t osites = fine_grid->oSites(); // D dimensional int64_t osites_hi = fine_mrhs_grid->oSites(); // nvec * osites const int Nsimd = vobj::Nsimd(); @@ -312,13 +356,14 @@ public: Lexicographic::IndexFromCoor(coor_c,sc,coarse_l); Lexicographic::IndexFromCoor(coor_b,sb,block_l); - scalar_object data = extractLane(lane,fineData[sfr]); + fine_scalar_object data = extractLane(lane,fineData[sfr]); int64_t site = (sc*lnvec + v)*bv + sb; - scalar_object * ptr = (scalar_object *)&blasData_p[site*lwords]; - *ptr = data; + // element-wise: the store may differ in precision from the field + const fine_scalar *dp = (const fine_scalar *)&data; + for(uint64_t w=0;w &vec_mrhs, deviceVector &blas) { typedef typename vobj::scalar_object coarse_scalar_object; + typedef typename vobj::scalar_type coarse_scalar; // may differ in precision from the BLAS scalar GridBase *coarse_mrhs_grid = vec_mrhs.Grid(); int _ndimension = coarse_grid->_ndimension; @@ -365,7 +411,7 @@ public: int64_t osites_hi = coarse_mrhs_grid->oSites(); // nvec * osites const int Nsimd = vobj::Nsimd(); - uint64_t cwords=sizeof(typename vobj::scalar_object)/sizeof(scalar); + uint64_t cwords=sizeof(coarse_scalar_object)/sizeof(coarse_scalar); GRID_ASSERT(cwords==nbasis); int64_t lnvec = nvec; @@ -391,8 +437,11 @@ public: Lexicographic::IndexFromCoor(coor_c,sc,coarse_l); int64_t blas_site = (sc*lnvec + v)*cwords; - coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site]; - coarse_scalar_object data = *ptr; + // Element-wise, converting: the BLAS buffer is in the FINE precision, + // the coarse field in its own. + coarse_scalar_object data; + coarse_scalar *dp = (coarse_scalar *)&data; + for(uint64_t b=0;b fine mrhs field, coarse mrhs field -> BLAS_C //////////////////////////////////////////////////////////////////////////// - void ExportFineGridMrhsVectors(Field &vec_mrhs, deviceVector &blas) + template + void ExportFineGridMrhsVectors(Lattice &vec_mrhs, deviceVector &blas) { - typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object fine_scalar_object; + typedef typename vobj::scalar_type fine_scalar; GridBase *fine_mrhs_grid = vec_mrhs.Grid(); int _ndimension = coarse_grid->_ndimension; GRID_ASSERT(fine_mrhs_grid->_ndimension == _ndimension+1); GRID_ASSERT(fine_mrhs_grid->_simd_layout[0] == 1); + // Full LOCAL dimensions must agree; the SIMD layout and precision need not for(int d=0;d<_ndimension;d++){ - GRID_ASSERT(fine_mrhs_grid->_rdimensions[d+1] == fine_grid->_rdimensions[d]); - GRID_ASSERT(fine_mrhs_grid->_simd_layout[d+1] == fine_grid->_simd_layout[d]); + GRID_ASSERT(fine_mrhs_grid->_ldimensions[d+1] == fine_ldimensions[d]); } + GRID_ASSERT(sizeof(fine_scalar_object)/sizeof(fine_scalar) == words); int nvec = fine_mrhs_grid->_rdimensions[0]; Coordinate fine_mrhs_rdimensions = fine_mrhs_grid->_rdimensions; - Coordinate fine_rdimensions = fine_grid->_rdimensions; + Coordinate fine_rdimensions(_ndimension); + Coordinate fsimd(_ndimension); + for(int d=0;d<_ndimension;d++){ + fine_rdimensions[d] = fine_mrhs_grid->_rdimensions[d+1]; + fsimd[d] = fine_mrhs_grid->_simd_layout[d+1]; + } Coordinate coarse_l = coarse_ldimensions; Coordinate block_l = block_ldimensions; - Coordinate fsimd = fine_simd; int64_t bv= block_vol; autoView( fineData , vec_mrhs, AcceleratorWrite); auto blasData_p = &blas[0]; auto fineData_p = &fineData[0]; - int64_t osites = fine_grid->oSites(); int64_t osites_hi = fine_mrhs_grid->oSites(); const int Nsimd = vobj::Nsimd(); @@ -468,8 +523,10 @@ public: int64_t site = (sc*lnvec + v)*bv + sb; - scalar_object * ptr = (scalar_object *)&blasData_p[site*lwords]; - scalar_object data = *ptr; + // element-wise: the store may differ in precision from the field + fine_scalar_object data; + fine_scalar *dp = (fine_scalar *)&data; + for(uint64_t w=0;w &vec_mrhs, deviceVector &blas) { typedef typename vobj::scalar_object coarse_scalar_object; + typedef typename vobj::scalar_type coarse_scalar; // may differ in precision from the BLAS scalar GridBase *coarse_mrhs_grid = vec_mrhs.Grid(); int _ndimension = coarse_grid->_ndimension; @@ -508,7 +566,7 @@ public: int64_t osites_hi = coarse_mrhs_grid->oSites(); const int Nsimd = vobj::Nsimd(); - uint64_t cwords=sizeof(typename vobj::scalar_object)/sizeof(scalar); + uint64_t cwords=sizeof(coarse_scalar_object)/sizeof(coarse_scalar); GRID_ASSERT(cwords==nbasis); int64_t lnvec = nvec; @@ -536,8 +594,8 @@ public: coarse_scalar_object data = extractLane(lane,coarseData[scr]); int64_t blas_site = (sc*lnvec + v)*cwords; - coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site]; - *ptr = data; + const coarse_scalar *dp = (const coarse_scalar *)&data; + for(uint64_t b=0;b &vecs, deviceVector &blas) + template + void ExportFineGridVectors(std::vector > &vecs, deviceVector &blas) { GRID_TRACE("ExportFineGridVectors"); - typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_object fine_scalar_object; + typedef typename vobj::scalar_type fine_scalar; int nvec = vecs.size(); - GRID_ASSERT(vecs[0].Grid()==fine_grid); - + GridBase *fgrid = vecs[0].Grid(); int _ndimension = coarse_grid->_ndimension; + // Full LOCAL dimensions must agree; the SIMD layout and precision need not + for(int d=0;d<_ndimension;d++) GRID_ASSERT(fgrid->_ldimensions[d] == fine_ldimensions[d]); + GRID_ASSERT(sizeof(fine_scalar_object)/sizeof(fine_scalar) == words); - Coordinate fine_rdimensions = fine_grid->_rdimensions; + Coordinate fine_rdimensions = fgrid->_rdimensions; Coordinate coarse_l = coarse_ldimensions; Coordinate block_l = block_ldimensions; - Coordinate fsimd = fine_simd; + Coordinate fsimd = fgrid->_simd_layout; // std::cout << " export fine Blas norm "<oSites(); + int64_t osites = fgrid->oSites(); uint64_t lwords = words; // std::cout << " Nsimd is "< &vecs) + template + void ImportBasis(std::vector < Lattice > &vecs) { // std::cout << " BlockProjector Import basis size "< - void blockProject(std::vector &fine,std::vector< Lattice > & coarse) + template + void blockProject(std::vector > &fine,std::vector< Lattice > & coarse) { GRID_TRACE("BlockProject"); int nrhs=fine.size(); - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); // std::cout << "blockProject nbasis " < - void blockPromote(std::vector &fine,std::vector > & coarse) + template + void blockPromote(std::vector > &fine,std::vector > & coarse) { GRID_TRACE("BlockPromote"); int nrhs=fine.size(); - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); BLAS_F.resize (fine_vol * words * nrhs ); @@ -876,12 +942,12 @@ public: // multiRHS ordered interfaces. The GEMM is identical; only the import and // export differ, so those are the whole of the layout question. //////////////////////////////////////////////////////////////////////////// - template - void blockProject(Field &fine_mrhs,Lattice &coarse_mrhs) + template + void blockProject(Lattice &fine_mrhs,Lattice &coarse_mrhs) { GRID_TRACE("BlockProjectMrhs"); int nrhs = fine_mrhs.Grid()->_rdimensions[0]; - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs); @@ -893,12 +959,12 @@ public: ExportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C); } - template - void blockPromote(Field &fine_mrhs,Lattice &coarse_mrhs) + template + void blockPromote(Lattice &fine_mrhs,Lattice &coarse_mrhs) { GRID_TRACE("BlockPromoteMrhs"); int nrhs = fine_mrhs.Grid()->_rdimensions[0]; - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs); @@ -914,12 +980,12 @@ public: // Mixed orderings. A single RHS fine operator produces a vector of fine // fields with no packing; the coarse side is still wanted in mrhs order. //////////////////////////////////////////////////////////////////////////// - template - void blockProject(std::vector &fine,Lattice &coarse_mrhs) + template + void blockProject(std::vector > &fine,Lattice &coarse_mrhs) { GRID_TRACE("BlockProjectMixed"); int nrhs = fine.size(); - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs); @@ -931,12 +997,12 @@ public: ExportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C); } - template - void blockProject(Field &fine_mrhs,std::vector< Lattice > &coarse) + template + void blockProject(Lattice &fine_mrhs,std::vector< Lattice > &coarse) { GRID_TRACE("BlockProjectMixed"); int nrhs = fine_mrhs.Grid()->_rdimensions[0]; - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse.size()==nrhs); @@ -948,12 +1014,12 @@ public: ExportCoarseGridVectors(coarse,BLAS_C); } - template - void blockPromote(std::vector &fine,Lattice &coarse_mrhs) + template + void blockPromote(std::vector > &fine,Lattice &coarse_mrhs) { GRID_TRACE("BlockPromoteMixed"); int nrhs = fine.size(); - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs); @@ -965,12 +1031,12 @@ public: ExportFineGridVectors(fine,BLAS_F); } - template - void blockPromote(Field &fine_mrhs,std::vector< Lattice > &coarse) + template + void blockPromote(Lattice &fine_mrhs,std::vector< Lattice > &coarse) { GRID_TRACE("BlockPromoteMixed"); int nrhs = fine_mrhs.Grid()->_rdimensions[0]; - int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar); + int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(typename cobj::scalar_type); GRID_ASSERT(nbasis==_nbasis); GRID_ASSERT(coarse.size()==nrhs); diff --git a/Grid/algorithms/deflation/MultiRHSDeflation.h b/Grid/algorithms/deflation/MultiRHSDeflation.h index 069390f4a..3931d2d25 100644 --- a/Grid/algorithms/deflation/MultiRHSDeflation.h +++ b/Grid/algorithms/deflation/MultiRHSDeflation.h @@ -59,6 +59,7 @@ public: typedef typename Field::scalar_type scalar; typedef typename Field::scalar_object scalar_object; + typedef typename Field::vector_object vobj; int nev; std::vector eval; @@ -80,10 +81,18 @@ public: grid=nullptr; vol=0; words=0; - BLAS_E.resize(0); - BLAS_R.resize(0); - BLAS_C.resize(0); - BLAS_G.resize(0); + // deviceVector is a std::vector with a device allocator: resize(0) drops + // the size and keeps the capacity, returning no device memory. Swapping + // with an empty vector destroys the buffer. + deviceVector().swap(BLAS_E); + deviceVector().swap(BLAS_R); + deviceVector().swap(BLAS_C); + deviceVector().swap(BLAS_G); + } + // Resident (non-evictable) device memory held by this deflator. + uint64_t DeviceBytes(void) + { + return (BLAS_E.capacity()+BLAS_R.capacity()+BLAS_C.capacity()+BLAS_G.capacity())*sizeof(scalar); } void Allocate(int _nev,GridBase *_grid) { @@ -123,31 +132,98 @@ public: ImportEigenVector(evec[_ev0+e],_eval[_ev0+e],e); } } + ///////////////////////////////////////////////////////////////////////// + // Sources as a vector of D-dimensional fields: each is one contiguous + // column of BLAS_R already, so import/export are straight copies. + // The eigenvectors may live on a D grid or on the Nrhs=1 D+1 grid of the + // same local volume; the two are the same bytes, so only the site count + // is checked, not grid identity. + ///////////////////////////////////////////////////////////////////////// void DeflateSources(std::vector &source,std::vector & guess) { int nrhs = source.size(); GRID_ASSERT(source.size()==guess.size()); - GRID_ASSERT(grid == guess[0].Grid()); + GRID_ASSERT(grid->lSites() == guess[0].Grid()->lSites()); conformable(guess[0],source[0]); int64_t vw = vol * words; - - RealD t0 = usecond(); BLAS_R.resize(nrhs * vw); // cost free if size doesn't change BLAS_G.resize(nrhs * vw); // cost free if size doesn't change BLAS_C.resize(nev * nrhs);// cost free if size doesn't change - ///////////////////////////////////////////// - // Copy in the multi-rhs sources - ///////////////////////////////////////////// - // for(int r=0;r }; -template -class TwoLevelADEF2 : public TwoLevelCG -{ - public: - /////////////////////////////////////////////////////////////////////////////////// - // Need something that knows how to get from Coarse to fine and back again - // void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){ - // void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){ - /////////////////////////////////////////////////////////////////////////////////// - GridBase *coarsegrid; - Aggregation &_Aggregates; - LinearFunction &_CoarseSolver; - LinearFunction &_CoarseSolverPrecise; - /////////////////////////////////////////////////////////////////////////////////// - - // more most opertor functions - TwoLevelADEF2(RealD tol, - Integer maxit, - LinearOperatorBase &FineLinop, - LinearFunction &Smoother, - LinearFunction &CoarseSolver, - LinearFunction &CoarseSolverPrecise, - Aggregation &Aggregates - ) : - TwoLevelCG(tol,maxit,FineLinop,Smoother,Aggregates.FineGrid), - _CoarseSolver(CoarseSolver), - _CoarseSolverPrecise(CoarseSolverPrecise), - _Aggregates(Aggregates) - { - coarsegrid = Aggregates.CoarseGrid; - }; - - virtual void PcgM1(Field & in, Field & out) - { - GRID_TRACE("MultiGridPreconditioner "); - // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] - - Field tmp(this->grid); - Field Min(this->grid); - CoarseField PleftProj(this->coarsegrid); - CoarseField PleftMss_proj(this->coarsegrid); - - this->SmoothTimer.Start(); - this->_Smoother(in,Min); - this->SmoothTimer.Stop(); - this->SmoothCalls++; - - this->MatrixTimer.Start(); - this->_FineLinop.HermOp(Min,out); - this->MatrixTimer.Stop(); - this->MatrixCalls++; - axpy(tmp,-1.0,out,in); // tmp = in - A Min - - this->ProjectTimer.Start(); - this->_Aggregates.ProjectToSubspace(PleftProj,tmp); - this->ProjectTimer.Stop(); - this->ProjectCalls++; - this->CoarseTimer.Start(); - this->_CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s - this->CoarseTimer.Stop(); - this->CoarseCalls++; - this->PromoteTimer.Start(); - this->_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min] - this->PromoteTimer.Stop(); - this->PromoteCalls++; - - axpy(out,1.0,Min,tmp); // Min+tmp - } - - virtual void Vstart(Field & x,const Field & src) - { - std::cout << GridLogMessage<<"HDCG: fPcg Vstart "<grid); - Field mmp(this->grid); - CoarseField PleftProj(this->coarsegrid); - CoarseField PleftMss_proj(this->coarsegrid); - - std::cout << GridLogMessage<<"HDCG: fPcg Vstart projecting "<_Aggregates.ProjectToSubspace(PleftProj,src); - std::cout << GridLogMessage<<"HDCG: fPcg Vstart coarse solve "<_CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s - std::cout << GridLogMessage<<"HDCG: fPcg Vstart promote "<_Aggregates.PromoteFromSubspace(PleftMss_proj,x); - - } - -}; template diff --git a/Grid/algorithms/iterative/AdefMrhs.h b/Grid/algorithms/iterative/AdefMrhs.h index e20900099..b18db9894 100644 --- a/Grid/algorithms/iterative/AdefMrhs.h +++ b/Grid/algorithms/iterative/AdefMrhs.h @@ -27,9 +27,10 @@ Author: Peter Boyle /* END LEGAL */ #pragma once +#include /* - * Compared to Tang-2009: P=Pleft. P^T = PRight Q=MssInv. + * Compared to Tang-2009: P=Pleft. P^T = PRight Q=MssInv. * Script A = SolverMatrix * Script P = Preconditioner * @@ -42,58 +43,64 @@ Author: Peter Boyle NAMESPACE_BEGIN(Grid); +////////////////////////////////////////////////////////////////////// +// Two-level CG family on a vector of right-hand sides, with the +// preconditioner (M1 and Vstart) as an MrhsPreconditioner object: +// fPcg flexible PCG, one Krylov per rhs, shared M1 call +// PrecBlockCGrQ preconditioned BlockCGrQ (arXiv:2409.03904 s2.2); +// assumes a stationary preconditioner; not the +// production choice (its linalg grows as nrhs^2) +// Also a LinearFunction: one rhs is the nrhs=1 case. +////////////////////////////////////////////////////////////////////// +enum class MrhsCGAlgorithm { fPcg, PrecBlockCGrQ }; + template -class TwoLevelCGmrhs +class TwoLevelCGmrhs : public LinearFunction { public: + using LinearFunction::operator(); RealD Tolerance; Integer MaxIterations; GridBase *grid; + MrhsCGAlgorithm Algorithm; - // Fine operator, Smoother, CoarseSolver - LinearOperatorBase &_FineLinop; - LinearFunction &_Smoother; - MultiRHSBlockCGLinalg _BlockCGLinalg; + LinearOperatorBase &_FineLinop; + MrhsPreconditioner &_Precon; + MultiRHSBlockCGLinalg _BlockCGLinalg; - GridStopWatch ProjectTimer; - GridStopWatch PromoteTimer; - GridStopWatch DeflateTimer; - GridStopWatch CoarseTimer; - GridStopWatch FineTimer; - GridStopWatch SmoothTimer; - GridStopWatch InsertTimer; - - /* - Field rrr; - Field sss; - Field qqq; - Field zzz; - */ - // more most opertor functions TwoLevelCGmrhs(RealD tol, Integer maxit, - LinearOperatorBase &FineLinop, - LinearFunction &Smoother, - GridBase *fine) : - Tolerance(tol), + LinearOperatorBase &FineLinop, + MrhsPreconditioner &Precon, + GridBase *fine, + MrhsCGAlgorithm alg = MrhsCGAlgorithm::fPcg) : + Tolerance(tol), MaxIterations(maxit), + Algorithm(alg), _FineLinop(FineLinop), - _Smoother(Smoother) - /* - rrr(fine), - sss(fine), - qqq(fine), - zzz(fine) -*/ + _Precon(Precon) { - grid = fine; + grid = fine; }; - - // Vector case + + // mrhs entry point virtual void operator() (std::vector &src, std::vector &x) { - SolveSingleSystem(src,x); - // SolvePrecBlockCG(src,x); + if ( Algorithm == MrhsCGAlgorithm::PrecBlockCGrQ ) SolvePrecBlockCG(src,x); + else SolveSingleSystem(src,x); + } + // LinearFunction entry points + virtual void operator() (const Field &src, Field &x) + { + std::vector S(1,src); + std::vector X(1,x); + (*this)(S,X); + x = X[0]; + } + virtual void operator() (const std::vector &src, std::vector &x) + { + std::vector S(src); + (*this)(S,x); } //////////////////////////////////////////////////////////////////////////////////////////////////// @@ -220,7 +227,7 @@ class TwoLevelCGmrhs ////////////////////////// // x0 = Vstart -- possibly modify guess ////////////////////////// - Vstart(X,src); + _Precon.Vstart(X,src); ////////////////////////// // R = B-AX @@ -234,7 +241,7 @@ class TwoLevelCGmrhs ////////////////////////////////// // Compute MZ = M1 Z = M1 B - M1 A x0 ////////////////////////////////// - PcgM1(Z,MZ); + _Precon(Z,MZ); ////////////////////////////////// // QC = Z @@ -248,13 +255,7 @@ class TwoLevelCGmrhs std::cout << GridLogMessage<<"PrecBlockCGrQ vec computed initial residual and QR fact " <_M)^{-1} inner prod, generalising Saad derivation of Precon CG //////////////////// @@ -326,7 +326,6 @@ class TwoLevelCGmrhs //////////////////////////// m_rr = m_C.adjoint() * m_C; - FineTimer.Stop(); RealD max_resid=0; RealD rrsum=0; @@ -350,14 +349,7 @@ class TwoLevelCGmrhs std::cout< ssq(nrhs); std::vector rsq(nrhs); @@ -449,13 +441,7 @@ class TwoLevelCGmrhs // std::cout << GridLogMessage<<"mrhs HDCG: "< & in,std::vector & out) = 0; - virtual void Vstart(std::vector & x,std::vector & src) = 0; - virtual void PcgM2(const Field & in, Field & out) { - out=in; - } - - virtual RealD PcgM3(const Field & p, Field & mmp){ + RealD PcgM3(const Field & p, Field & mmp){ RealD dd; _FineLinop.HermOp(p,mmp); ComplexD dot = innerProduct(p,mmp); dd=real(dot); return dd; } - }; -template -class TwoLevelADEF2mrhs : public TwoLevelCGmrhs +////////////////////////////////////////////////////////////////////// +// ADEF-2 as an mrhs preconditioner object (Tang, Nabben, Vuik, Erlangga): +// M1 = [1 - Q A] M + Q with Q = P A_c^{-1} P^dag (deflated coarse solve) +// Vstart = Q b the coarse-corrected start (precise coarse solve) +// The coarse space is the D+1 multiRHS field throughout: the projector's +// mixed overloads take the vector of fine fields straight to it, and the +// deflator and coarse solver act on it. No per-rhs slices anywhere. +// The D+1 grid follows the solver's Nrhs through SetCoarseGridMrhs. +////////////////////////////////////////////////////////////////////// +// Projector_t defaults to a transfer operator whose STORE matches Field. +template > +class MrhsADEF2Preconditioner : public MrhsPreconditioner { public: - GridBase *coarsegrid; - GridBase *coarsegridmrhs; - LinearFunction &_CoarseSolverMrhs; - LinearFunction &_CoarseSolverPreciseMrhs; - MultiRHSBlockProject &_Projector; - MultiRHSDeflation &_Deflator; + GridBase *grid; // fine + GridBase *coarsegridmrhs; // D+1, rhs innermost + LinearOperatorBase &_FineLinop; + LinearFunction &_Smoother; + LinearFunction &_CoarseSolver; // in M1 + LinearFunction &_CoarseSolverPrecise; // in Vstart + Projector_t &_Projector; // store precision need not match Field + MultiRHSDeflation &_Deflator; // nev==0: no deflation, zero coarse guess - - TwoLevelADEF2mrhs(RealD tol, - Integer maxit, - LinearOperatorBase &FineLinop, - LinearFunction &Smoother, - LinearFunction &CoarseSolverMrhs, - LinearFunction &CoarseSolverPreciseMrhs, - MultiRHSBlockProject &Projector, - MultiRHSDeflation &Deflator, - GridBase *_coarsemrhsgrid) : - TwoLevelCGmrhs(tol, maxit,FineLinop,Smoother,Projector.fine_grid), - _CoarseSolverMrhs(CoarseSolverMrhs), - _CoarseSolverPreciseMrhs(CoarseSolverPreciseMrhs), - _Projector(Projector), - _Deflator(Deflator) + MrhsADEF2Preconditioner(LinearOperatorBase &FineLinop, + LinearFunction &Smoother, + LinearFunction &CoarseSolver, + LinearFunction &CoarseSolverPrecise, + Projector_t &Projector, + MultiRHSDeflation &Deflator, + GridBase *CoarseGridMrhs, + GridBase *FineGrid) : + _FineLinop(FineLinop), _Smoother(Smoother), + _CoarseSolver(CoarseSolver), _CoarseSolverPrecise(CoarseSolverPrecise), + _Projector(Projector), _Deflator(Deflator) { - coarsegrid = Projector.coarse_grid; - coarsegridmrhs = _coarsemrhsgrid;// Thi could be in projector + // The fine grid is given, NOT taken from the projector: the projector's + // grid carries its STORE's layout, which need not be this Field's (an + // fp32 preconditioner may project through an fp64 store, or the reverse). + grid = FineGrid; + coarsegridmrhs = CoarseGridMrhs; }; + // Store matches Field: the projector's grid IS this Field's grid. + MrhsADEF2Preconditioner(LinearOperatorBase &FineLinop, + LinearFunction &Smoother, + LinearFunction &CoarseSolver, + LinearFunction &CoarseSolverPrecise, + Projector_t &Projector, + MultiRHSDeflation &Deflator, + GridBase *CoarseGridMrhs) : + MrhsADEF2Preconditioner(FineLinop,Smoother,CoarseSolver,CoarseSolverPrecise, + Projector,Deflator,CoarseGridMrhs,Projector.fine_grid) {}; - // Override Vstart - virtual void Vstart(std::vector & x,std::vector & src) + void SetCoarseGridMrhs(GridBase *CoarseGridMrhs) { coarsegridmrhs = CoarseGridMrhs; } + + // The coarse solver's starting guess: deflated, or zero + void CoarseGuess(CoarseField &src,CoarseField &guess) { - int nrhs=x.size(); - /////////////////////////////////// - // Choose x_0 such that - // x_0 = guess + (A_ss^inv) r_s = guess + Ass_inv [src -Aguess] - // = [1 - Ass_inv A] Guess + Assinv src - // = P^T guess + Assinv src - // = Vstart [Tang notation] - // This gives: - // W^T (src - A x_0) = src_s - A guess_s - r_s - // = src_s - (A guess)_s - src_s + (A guess)_s - // = 0 - /////////////////////////////////// - std::vector PleftProj(nrhs,this->coarsegrid); - std::vector PleftMss_proj(nrhs,this->coarsegrid); - CoarseField PleftProjMrhs(this->coarsegridmrhs); - CoarseField PleftMss_projMrhs(this->coarsegridmrhs); - - this->_Projector.blockProject(src,PleftProj); - this->_Deflator.DeflateSources(PleftProj,PleftMss_proj); - for(int rhs=0;rhs_CoarseSolverPreciseMrhs(PleftProjMrhs,PleftMss_projMrhs); // Ass^{-1} r_s - - for(int rhs=0;rhs_Projector.blockPromote(x,PleftMss_proj); - } - - virtual void PcgM1(std::vector & in,std::vector & out){ - - int nrhs=in.size(); - - // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] - std::vector tmp(nrhs,this->grid); - std::vector Min(nrhs,this->grid); - - std::vector PleftProj(nrhs,this->coarsegrid); - std::vector PleftMss_proj(nrhs,this->coarsegrid); - - CoarseField PleftProjMrhs(this->coarsegridmrhs); - CoarseField PleftMss_projMrhs(this->coarsegridmrhs); - - // this->rrr=in[0]; - -#undef SMOOTHER_BLOCK_SOLVE -#if SMOOTHER_BLOCK_SOLVE - this->SmoothTimer.Start(); - this->_Smoother(in,Min); - this->SmoothTimer.Stop(); -#else - for(int rhs=0;rhsSmoothTimer.Start(); - this->_Smoother(in[rhs],Min[rhs]); - this->SmoothTimer.Stop(); - } -#endif - // this->sss=Min[0]; - - for(int rhs=0;rhsFineTimer.Start(); - this->_FineLinop.HermOp(Min[rhs],out[rhs]); - axpy(tmp[rhs],-1.0,out[rhs],in[rhs]); // resid = in - A Min - this->FineTimer.Stop(); - - } - - this->ProjectTimer.Start(); - this->_Projector.blockProject(tmp,PleftProj); - this->ProjectTimer.Stop(); this->DeflateTimer.Start(); - this->_Deflator.DeflateSources(PleftProj,PleftMss_proj); + if ( _Deflator.nev > 0 ) _Deflator.DeflateSources(src,guess); + else guess = Zero(); this->DeflateTimer.Stop(); - this->InsertTimer.Start(); - for(int rhs=0;rhsInsertTimer.Stop(); + } + // x_0 = Q b + virtual void Vstart(std::vector &x,std::vector &src) + { + CoarseField Csrc(coarsegridmrhs), Csol(coarsegridmrhs); + this->ProjectTimer.Start(); _Projector.blockProject(src,Csrc); this->ProjectTimer.Stop(); + CoarseGuess(Csrc,Csol); + this->CoarseTimer.Start(); _CoarseSolverPrecise(Csrc,Csol); this->CoarseTimer.Stop(); + this->PromoteTimer.Start(); _Projector.blockPromote(x,Csol); this->PromoteTimer.Stop(); + } + // [1 - Q A] M in + Q in = Min + Q [in - A Min] + virtual void operator()(std::vector &in,std::vector &out) + { + int nrhs=in.size(); + std::vector tmp(nrhs,grid), Min(nrhs,grid); + CoarseField Csrc(coarsegridmrhs), Csol(coarsegridmrhs); - this->CoarseTimer.Start(); - this->_CoarseSolverMrhs(PleftProjMrhs,PleftMss_projMrhs); // Ass^{-1} [in - A Min]_s - this->CoarseTimer.Stop(); + this->SmoothTimer.Start(); + for(int r=0;rSmoothTimer.Stop(); - this->InsertTimer.Start(); - for(int rhs=0;rhsInsertTimer.Stop(); - this->PromoteTimer.Start(); - this->_Projector.blockPromote(tmp,PleftMss_proj);// tmp= Q[in - A Min] - this->PromoteTimer.Stop(); this->FineTimer.Start(); - // this->qqq=tmp[0]; - for(int rhs=0;rhszzz=out[0]; + this->FineTimer.Stop(); + + this->ProjectTimer.Start(); _Projector.blockProject(tmp,Csrc); this->ProjectTimer.Stop(); + CoarseGuess(Csrc,Csol); + this->CoarseTimer.Start(); _CoarseSolver(Csrc,Csol); this->CoarseTimer.Stop(); + this->PromoteTimer.Start(); _Projector.blockPromote(tmp,Csol); this->PromoteTimer.Stop(); + + this->FineTimer.Start(); + for(int r=0;rFineTimer.Stop(); } }; - NAMESPACE_END(Grid); - - diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h index c7f61b217..1feef6691 100644 --- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h +++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h @@ -57,6 +57,7 @@ public: GridStopWatch MatTimer; GridStopWatch LinalgTimer; std::string name; + std::string trace_step = "PGCR_step"; // name + step; see Name() int ZeroGuess = 0; // caller contract: guess is always zero => first-cycle r0 = src, skip the apply // persistent GCR history (see GCRnStep) GridBase *hist_grid = nullptr; @@ -68,7 +69,10 @@ public: LinearFunction &Preconditioner; LinearOperatorBase &Linop; - void Name(std::string _name) { name = _name; }; + // 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; } @@ -222,7 +226,7 @@ public: GCRLogLevel<< "PGCR true residual "<< sqrt(cp/SSQ) < &hermop, - TwoLevelADEF2mrhs & theHDCG, + TwoLevelCGmrhs & theHDCG, int nrhs) { std::vector src_mrhs(nrhs,FineGrid); diff --git a/Grid/algorithms/multigrid/BlockCyclic.h b/Grid/algorithms/multigrid/BlockCyclic.h index c8c7879be..c27f53768 100644 --- a/Grid/algorithms/multigrid/BlockCyclic.h +++ b/Grid/algorithms/multigrid/BlockCyclic.h @@ -21,6 +21,20 @@ Author: Peter Boyle NAMESPACE_BEGIN(Grid); +/////////////////////////////////////////////////////////////////////////////// +// Scalar of the distributed dense inversion: block-cyclic redistribution, +// SUMMA, the Schur recursion and its LU leaf. A configure-time choice +// (--enable-dense-inverse-precision=double|single), INDEPENDENT of the +// coarse-space precision: the matrix elements arrive in the coarse scalar +// and the apply slab is fp32 either way; only the factorisation's +// arithmetic and footprint follow this type. +/////////////////////////////////////////////////////////////////////////////// +#ifdef GRID_DENSE_INVERSE_SINGLE +typedef ComplexF DenseInverseScalar; +#else +typedef ComplexD DenseInverseScalar; +#endif + /////////////////////////////////////////////////////////////////////////////// // BlockCyclicLayout: the index arithmetic of a 2D block-cyclic distribution // of an N x N matrix over a Pr x Pc logical process grid with block size nb. diff --git a/Grid/algorithms/multigrid/BlockCyclicRedistribute.h b/Grid/algorithms/multigrid/BlockCyclicRedistribute.h index 69eb24d45..7de843cd5 100644 --- a/Grid/algorithms/multigrid/BlockCyclicRedistribute.h +++ b/Grid/algorithms/multigrid/BlockCyclicRedistribute.h @@ -34,7 +34,7 @@ NAMESPACE_BEGIN(Grid); class BlockRows { public: - deviceVector data; + deviceVector data; int64_t rows; int64_t cols; int64_t ld; @@ -52,7 +52,7 @@ public: ld = r; data.resize((uint64_t)r*c); } - ComplexD *ColumnWindow(int64_t col0) + DenseInverseScalar *ColumnWindow(int64_t col0) { GRID_ASSERT( col0 >= 0 ); GRID_ASSERT( col0 <= cols ); @@ -128,10 +128,10 @@ public: // buffer(a,b) = elem(rows[a], cols[b]), a fastest. ///////////////////////////////////////////////////////////////////////// static void MoveEdge(int toBuffer, - ComplexD *mat, int64_t ld, + DenseInverseScalar *mat, int64_t ld, const std::vector &roff, // per-row offset in mat const std::vector &coff, // per-col offset in mat - ComplexD *buf) + DenseInverseScalar *buf) { int64_t nr = roff.size(); int64_t nc = coff.size(); @@ -193,7 +193,7 @@ public: ///////////////////////////////////////////////////////////////////////// static void Redistribute(int dir, GridBase *grid, const std::vector &rowStart, - ComplexD *rows1d, int64_t myrows, + DenseInverseScalar *rows1d, int64_t myrows, BlockCyclicMatrix &A) { BlockCyclicLayout &L = A.layout; @@ -206,7 +206,7 @@ public: int64_t ld1 = myrows ? myrows : 1; std::vector rows, cols, roff, coff; - deviceVector sbuf(1), rbuf(1); + deviceVector sbuf(1), rbuf(1); /////////////////////////////////////////////////////////////////////// // Self edge: purely local, via a bounce buffer (shares all the code). @@ -261,7 +261,7 @@ public: } grid->SendToRecvFrom((void *)&sbuf[0], partner, (void *)&rbuf[0], partner, - nmax*sizeof(ComplexD)); + nmax*sizeof(DenseInverseScalar)); if ( nin ){ if ( dir > 0 ) { Offsets2D(L, irow, icol, roff, coff); MoveEdge(0, &A.data[0], L.mloc, roff, coff, &rbuf[0]); } @@ -272,11 +272,11 @@ public: } static void RowsToCyclic(GridBase *grid, const std::vector &rowStart, - ComplexD *rows1d, int64_t myrows, BlockCyclicMatrix &A) + DenseInverseScalar *rows1d, int64_t myrows, BlockCyclicMatrix &A) { Redistribute(+1, grid, rowStart, rows1d, myrows, A); } static void CyclicToRows(GridBase *grid, const std::vector &rowStart, - BlockCyclicMatrix &A, ComplexD *rows1d, int64_t myrows) + BlockCyclicMatrix &A, DenseInverseScalar *rows1d, int64_t myrows) { Redistribute(-1, grid, rowStart, rows1d, myrows, A); } }; diff --git a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h index 24451cc6e..718f4370a 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h +++ b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h @@ -115,7 +115,7 @@ public: // Window copy-scale: Dst[i0:i1, j0:j1] = alpha * Src[same window]. // Both share one layout, so the local bands coincide; pure local kernel. /////////////////////////////////////////////////////////////////////////// - void WindowCopyScale(ComplexD alpha, + void WindowCopyScale(DenseInverseScalar alpha, BlockCyclicMatrix &Src, BlockCyclicMatrix &Dst, int64_t i0, int64_t i1, int64_t j0, int64_t j1) { @@ -127,8 +127,8 @@ public: L.ColRange(j0,j1, lj0,lj1); int64_t m = li1-li0, n = lj1-lj0; if ( !(m && n) ) return; - ComplexD *src = Src.LocalWindow(li0,lj0); - ComplexD *dst = Dst.LocalWindow(li0,lj0); + DenseInverseScalar *src = Src.LocalWindow(li0,lj0); + DenseInverseScalar *dst = Dst.LocalWindow(li0,lj0); int64_t ldS = Src.layout.mloc; int64_t ldD = L.mloc; tCopy -= usecond(); @@ -164,10 +164,10 @@ public: GRID_ASSERT( lc1-lc0 == w ); // Pack the strided block dense (inverseBatched assumes lda == w). - deviceVector dense((uint64_t)w*w); + deviceVector dense((uint64_t)w*w); { - ComplexD *src = A.LocalWindow(lr0,lc0); - ComplexD *dst = &dense[0]; + DenseInverseScalar *src = A.LocalWindow(lr0,lc0); + DenseInverseScalar *dst = &dense[0]; int64_t ld = L.mloc; accelerator_for(idx, (uint64_t)(w*w), 1, { int64_t jj = idx / w; @@ -176,15 +176,15 @@ public: }); } { - deviceVector bp(1); - std::vector ptr(1); + deviceVector bp(1); + std::vector ptr(1); ptr[0] = &dense[0]; - acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD*)); + acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(DenseInverseScalar*)); INV.inverseBatched(w, bp); } { - ComplexD *src = &dense[0]; - ComplexD *dst = A.LocalWindow(lr0,lc0); + DenseInverseScalar *src = &dense[0]; + DenseInverseScalar *dst = A.LocalWindow(lr0,lc0); int64_t ld = L.mloc; accelerator_for(idx, (uint64_t)(w*w), 1, { int64_t jj = idx / w; @@ -194,8 +194,8 @@ public: } // Growth telemetry, local only. { - std::vector h((uint64_t)w*w); - acceleratorCopyFromDevice(&dense[0], &h[0], h.size()*sizeof(ComplexD)); + std::vector h((uint64_t)w*w); + acceleratorCopyFromDevice(&dense[0], &h[0], h.size()*sizeof(DenseInverseScalar)); double mx = 0.0; for(auto &z : h){ double re=z.real(), im=z.imag(); @@ -232,23 +232,23 @@ public: int64_t lr0,lr1,lc0,lc1; L.RowRange(c0,c1,lr0,lr1); L.ColRange(c0,c1,lc0,lc1); const int64_t mq = lr1-lr0, nq = lc1-lc0; - deviceVector dense; // root only: W x W column major - deviceVector piece, dummy; // piece: my mq x nq contiguous; dummy: reverse-direction filler + deviceVector dense; // root only: W x W column major + deviceVector piece, dummy; // piece: my mq x nq contiguous; dummy: reverse-direction filler if ( me == root ) dense.resize((uint64_t)W*W); - auto pack_piece = [&](ComplexD *dst, int64_t m, int64_t n, int64_t r0, int64_t cc0){ - ComplexD *src = A.LocalWindow(r0,cc0); const int64_t ld = L.mloc; + auto pack_piece = [&](DenseInverseScalar *dst, int64_t m, int64_t n, int64_t r0, int64_t cc0){ + DenseInverseScalar *src = A.LocalWindow(r0,cc0); const int64_t ld = L.mloc; accelerator_for(idx,(uint64_t)(m*n),1,{ int64_t jj=idx/m, ii=idx-jj*m; dst[ii+jj*m] = src[ii+jj*ld]; }); }; - auto unpack_piece = [&](ComplexD *src, int64_t m, int64_t n, int64_t r0, int64_t cc0){ - ComplexD *dst = A.LocalWindow(r0,cc0); const int64_t ld = L.mloc; + auto unpack_piece = [&](DenseInverseScalar *src, int64_t m, int64_t n, int64_t r0, int64_t cc0){ + DenseInverseScalar *dst = A.LocalWindow(r0,cc0); const int64_t ld = L.mloc; accelerator_for(idx,(uint64_t)(m*n),1,{ int64_t jj=idx/m, ii=idx-jj*m; dst[ii+jj*ld] = src[ii+jj*m]; }); }; // root: piece of rank q <-> dense, via the closed-form block map - auto root_place = [&](ComplexD *pc, int64_t m, int64_t n, int q, int to_dense){ + auto root_place = [&](DenseInverseScalar *pc, int64_t m, int64_t n, int q, int to_dense){ const int pq=q/Pc, cq=q%Pc; const int64_t brq0=FirstBlock(b0,pq,Pr), bcq0=FirstBlock(b0,cq,Pc); - ComplexD *dn = &dense[0]; const int64_t WW=W, NB=nb, PR=Pr, PC=Pc, C0=c0; + DenseInverseScalar *dn = &dense[0]; const int64_t WW=W, NB=nb, PR=Pr, PC=Pc, C0=c0; accelerator_for(idx,(uint64_t)(m*n),1,{ int64_t jj=idx/m, ii=idx-jj*m; int64_t gr = (brq0 + (ii/NB)*PR)*NB + ii%NB - C0; @@ -266,18 +266,18 @@ public: tBigGather -= usecond(); if ( mq*nq ) { piece.resize((uint64_t)mq*nq); dummy.resize((uint64_t)mq*nq); pack_piece(&piece[0],mq,nq,lr0,lc0); accelerator_barrier(); } if ( me == root ) { - deviceVector stage; + deviceVector stage; for(int q=0;q junk((uint64_t)m*n); - grid->SendToRecvFrom((void *)&junk[0], q, (void *)&stage[0], q, (uint64_t)m*n*sizeof(ComplexD)); + deviceVector junk((uint64_t)m*n); + grid->SendToRecvFrom((void *)&junk[0], q, (void *)&stage[0], q, (uint64_t)m*n*sizeof(DenseInverseScalar)); root_place(&stage[0],m,n,q,1); } accelerator_barrier(); } else if ( mq*nq ) { - grid->SendToRecvFrom((void *)&piece[0], root, (void *)&dummy[0], root, (uint64_t)mq*nq*sizeof(ComplexD)); + grid->SendToRecvFrom((void *)&piece[0], root, (void *)&dummy[0], root, (uint64_t)mq*nq*sizeof(DenseInverseScalar)); } tBigGather += usecond(); @@ -289,17 +289,17 @@ public: // ---- scatter ---- tBigScatter -= usecond(); if ( me == root ) { - deviceVector stage; + deviceVector stage; for(int q=0;q junk((uint64_t)m*n); + deviceVector junk((uint64_t)m*n); root_place(&stage[0],m,n,q,0); accelerator_barrier(); - grid->SendToRecvFrom((void *)&stage[0], q, (void *)&junk[0], q, (uint64_t)m*n*sizeof(ComplexD)); + grid->SendToRecvFrom((void *)&stage[0], q, (void *)&junk[0], q, (uint64_t)m*n*sizeof(DenseInverseScalar)); } } else if ( mq*nq ) { - grid->SendToRecvFrom((void *)&dummy[0], root, (void *)&piece[0], root, (uint64_t)mq*nq*sizeof(ComplexD)); + grid->SendToRecvFrom((void *)&dummy[0], root, (void *)&piece[0], root, (uint64_t)mq*nq*sizeof(DenseInverseScalar)); } if ( mq*nq ) { unpack_piece(&piece[0],mq,nq,lr0,lc0); accelerator_barrier(); } tBigScatter += usecond(); @@ -326,7 +326,7 @@ public: int64_t c0 = b0*L.nb; int64_t m = bm*L.nb; int64_t c1 = std::min(L.N, b1*L.nb); - ComplexD one (1.0,0.0), mone(-1.0,0.0), zero(0.0,0.0); + DenseInverseScalar one (1.0,0.0), mone(-1.0,0.0), zero(0.0,0.0); // 1. A11 -> A11inv SchurNode(A,Bt,Ct,Tt,Ut, b0,bm); diff --git a/Grid/algorithms/multigrid/BlockCyclicSumma.h b/Grid/algorithms/multigrid/BlockCyclicSumma.h index a8f92dac9..45c774518 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSumma.h +++ b/Grid/algorithms/multigrid/BlockCyclicSumma.h @@ -71,7 +71,7 @@ class BlockCyclicMatrix public: GridBase *grid; // borrowed, never owned BlockCyclicLayout layout; - deviceVector data; // column major, ld = layout.mloc + deviceVector data; // column major, ld = layout.mloc BlockCyclicMatrix(GridBase *g, int64_t N, int64_t nb, int Pr, int Pc) : grid(g), @@ -82,7 +82,7 @@ public: data.resize( sz ? sz : 1 ); } - ComplexD *LocalWindow(int64_t li, int64_t lj) + DenseInverseScalar *LocalWindow(int64_t li, int64_t lj) { return &data[0] + li + lj*layout.mloc; } @@ -93,32 +93,32 @@ public: // oracle only; production data enters through the direct block-cyclic // import, never through these. ///////////////////////////////////////////////////////////////////////// - void ImportGlobal(const std::vector &G) + void ImportGlobal(const std::vector &G) { int64_t N = layout.N; GRID_ASSERT( (int64_t)G.size() == N*N ); - std::vector h((uint64_t)layout.mloc*layout.nloc, ComplexD(0.0,0.0)); + std::vector h((uint64_t)layout.mloc*layout.nloc, DenseInverseScalar(0.0,0.0)); for(int64_t j=0;j &G) + void ExportGlobal(std::vector &G) { int64_t N = layout.N; - G.assign((uint64_t)N*N, ComplexD(0.0,0.0)); - std::vector h((uint64_t)layout.mloc*layout.nloc); + G.assign((uint64_t)N*N, DenseInverseScalar(0.0,0.0)); + std::vector h((uint64_t)layout.mloc*layout.nloc); if ( h.size() ) - acceleratorCopyFromDevice(&data[0], &h[0], h.size()*sizeof(ComplexD)); + acceleratorCopyFromDevice(&data[0], &h[0], h.size()*sizeof(DenseInverseScalar)); for(int64_t j=0;jGlobalSumVector((ComplexD *)&G[0], (int)(N*N)); // zero-fill: exact + if ( N ) grid->GlobalSumVector((DenseInverseScalar *)&G[0], (int)(N*N)); // zero-fill: exact } }; @@ -139,8 +139,8 @@ public: // per message: measured 1.26 GB/s/rank on 13 MB ring messages (production, // GRID_ALLOC_NCACHE_LARGE=64) against 62 s total for the same inverse when // hipMalloc returned a stable address. - deviceVector Abuf; - deviceVector Bbuf; + deviceVector Abuf; + deviceVector Bbuf; double tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0; uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0; // Per-message-size histogram (bucket = floor(log2 bytes)): count, bytes, @@ -156,10 +156,10 @@ public: static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1) { return (a0 < b1) && (b0 < a1); } - void Multiply(ComplexD alpha, + void Multiply(DenseInverseScalar alpha, BlockCyclicMatrix &A, BlockCyclicMatrix &B, - ComplexD beta, + DenseInverseScalar beta, BlockCyclicMatrix &C, int64_t i0, int64_t i1, int64_t j0, int64_t j1, @@ -220,8 +220,8 @@ public: if ( Bbuf.size() < std::max(slotB*Pr,1) ) Bbuf.resize( std::max(slotB*Pr,1) ); tAlloc += usecond(); - deviceVector ap(1), bp(1), cp(1); - std::vector ptr(1); + deviceVector ap(1), bp(1), cp(1); + std::vector ptr(1); int firstblock = 1; for(int64_t r0=kb0; r0= 0 ); GRID_ASSERT( idxs < S ); - ComplexD *src = B.LocalWindow(lr0, lj0); - ComplexD *dst = &Bbuf[0] + slotB*prow + slotB1*idxs; + DenseInverseScalar *src = B.LocalWindow(lr0, lj0); + DenseInverseScalar *dst = &Bbuf[0] + slotB*prow + slotB1*idxs; int64_t ld = B.layout.mloc; int64_t nn = nloc_j; accelerator_for(idx, (uint64_t)(nb_s*nn), 1, { @@ -284,9 +284,9 @@ public: double tm = usecond(); grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest, (void *)(&Abuf[0]+slotA*cr), src, - slotA*sizeof(ComplexD)); - HistAdd(slotA*sizeof(ComplexD), usecond()-tm); - bytesRing += slotA*sizeof(ComplexD); nRingMsg++; + slotA*sizeof(DenseInverseScalar)); + HistAdd(slotA*sizeof(DenseInverseScalar), usecond()-tm); + bytesRing += slotA*sizeof(DenseInverseScalar); nRingMsg++; } tRingA += usecond(); } @@ -304,9 +304,9 @@ public: double tm = usecond(); grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest, (void *)(&Bbuf[0]+slotB*rr), src, - slotB*sizeof(ComplexD)); - HistAdd(slotB*sizeof(ComplexD), usecond()-tm); - bytesRing += slotB*sizeof(ComplexD); nRingMsg++; + slotB*sizeof(DenseInverseScalar)); + HistAdd(slotB*sizeof(DenseInverseScalar), usecond()-tm); + bytesRing += slotB*sizeof(DenseInverseScalar); nRingMsg++; } tRingB += usecond(); } @@ -322,15 +322,15 @@ public: int cA = (int)(s%Pc); int rB = (int)(s%Pr); int64_t idxs = (s - r0 - ((rB - r0%Pr + Pr) % Pr)) / Pr; - ComplexD beta_use = firstblock ? beta : ComplexD(1.0,0.0); + DenseInverseScalar beta_use = firstblock ? beta : DenseInverseScalar(1.0,0.0); firstblock = 0; ptr[0] = &Abuf[0] + slotA*cA; - acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(ComplexD *)); + acceleratorCopyToDevice(&ptr[0], &ap[0], sizeof(DenseInverseScalar *)); ptr[0] = &Bbuf[0] + slotB*rB + slotB1*idxs; - acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(ComplexD *)); + acceleratorCopyToDevice(&ptr[0], &bp[0], sizeof(DenseInverseScalar *)); ptr[0] = C.LocalWindow(li0, lj0); - acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(ComplexD *)); + acceleratorCopyToDevice(&ptr[0], &cp[0], sizeof(DenseInverseScalar *)); BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N, (int)mloc_i, (int)nloc_j, (int)nb_s, diff --git a/Grid/algorithms/multigrid/DenseCoarseMatrix.h b/Grid/algorithms/multigrid/DenseCoarseMatrix.h index 67cd47fc2..88c3697d9 100644 --- a/Grid/algorithms/multigrid/DenseCoarseMatrix.h +++ b/Grid/algorithms/multigrid/DenseCoarseMatrix.h @@ -64,7 +64,8 @@ NAMESPACE_BEGIN(Grid); // VERIFY ||A Ainv x - x||/||x|| certifies the DEVICE slab + split-K path at the // end of Import, since the single-RHS apply routes through the same core. // -// Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD +// Tensor-depth agnostic: site scalar objects treated as contiguous coarse +// scalars (ComplexF or ComplexD, following the coefficient precision) // (iScalar wrappers add no data), so any MG level's coarse operator imports. ////////////////////////////////////////////////////////////////////////////////////// // @@ -85,6 +86,10 @@ public: typedef typename vobj::scalar_object sobj; typedef typename CoarseMatrix::vector_object Mvobj; typedef typename Mvobj::scalar_object Msobj; + // Scalar of the coarse site objects (ComplexF or ComplexD): the apply slab + // is always fp32 and the inversion always fp64, but the SOURCE of both -- + // the coarse operator's matrix elements -- carries this precision. + typedef typename GridTypeMapper::scalar_type CoarseScalar; GridBase *grid; int nd; @@ -115,11 +120,20 @@ public: std::vector hY; int NK; // split-K chunk count (divides N) + // Resident device memory: the apply slab and its staging. deviceVector is + // not evictable, so this counts against the hard budget, not the cache. + uint64_t DeviceBytes(void) + { + return (uint64_t)(dSlab.capacity()+dX.capacity()+dY.capacity()+dG.capacity()+dPartial.capacity())*sizeof(ComplexF) + + (uint64_t)dLex2Rank.capacity()*sizeof(int) + + (uint64_t)dRm2G.capacity()*sizeof(int64_t); + } + DenseCoarseMatrix(GridBase *g) : grid(g) { - GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(ComplexD) ); - GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(ComplexD) ); + GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(CoarseScalar) ); + GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(CoarseScalar) ); nd = grid->_ndimension; N = grid->gSites() * nbasis; lsites = grid->lSites(); @@ -228,7 +242,7 @@ public: //////////////////////////////////////////////////////////////////// { Field x(grid); Field y(grid); Field z(grid); - x = ComplexD(1.0,0.0); + x = CoarseScalar(1.0,0.0); double ta = usecond(); (*this)(x, y); double tb = usecond(); @@ -263,7 +277,7 @@ public: Field min(mgrid), mout(mgrid); for(int r=0;r= 1.0e-6 ) { + // Slices differ by rounding amplified by the operator's conditioning + // when the product cancels (A applied to A^-1 x): the tolerance follows + // the coarse precision (fp32 measured ~5e-6 here). A genuine rhs + // mix-up is O(1). + const RealD otol = (sizeof(CoarseScalar)==sizeof(ComplexF)) ? 1.0e-4 : 1.0e-6; + if ( rel >= otol ) { std::cout << GridLogMessage << "DenseCoarseMatrix: oracle rhs "<GlobalSumVector(&xh[0], (int)N); std::vector yh(nrows); @@ -403,7 +422,7 @@ public: }); for(int ss=0; ss= 1.0e-3 ) { std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE FAILED. If O(1), the " << "stencil shift-sign convention of the coarse operator has changed: " - << "the import in ImportDense/ImportDenseFP64 must change with it" + << "the import in ImportDense/ImportDenseForInversion must change with it" << std::endl; } GRID_ASSERT(rel < 1.0e-3); @@ -468,24 +487,25 @@ public: } //////////////////////////////////////////////////////////////////// - // 3b. Direct stencil -> fp64 rank-major import of MY ROWS of A (the - // end-to-end fp64 path: the stencil source IS ComplexD; nothing is - // rounded through fp32 on the way into the inversion). Same + // 3b. Direct stencil -> rank-major import of MY ROWS of A in the + // inversion precision (DenseInverseScalar): the stencil source is + // read once, straight into the buffer the Schur recursion factorises, + // never via the fp32 apply slab. Same // loop/sign/accumulate/transposed-contraction discipline as // ImportDense; output is column-major rows x N with columns in // rank-major order (g2rm). - // ALWAYS-ON CERTIFICATE: the fp64 import, rounded, must agree with + // ALWAYS-ON CERTIFICATE: this import, rounded, must agree with // the fp32 slab entry at the corresponding global column, over the // WHOLE of my rows (few ulp: wrapped-shift collisions accumulate in // different precision order). NaN-proof: non-finite entries are // counted explicitly since max() silently masks NaN. //////////////////////////////////////////////////////////////////// template - void ImportDenseFP64(CoarseOp &Op, BlockRows &S, std::vector &g2rm) + void ImportDenseForInversion(CoarseOp &Op, BlockRows &S, std::vector &g2rm) { Coordinate gdims = grid->GlobalDimensions(); - std::vector h((uint64_t)nrows*N, ComplexD(0.0,0.0)); + std::vector h((uint64_t)nrows*N, DenseInverseScalar(0.0,0.0)); for(int p=0; pGlobalMax(gmx); grid->GlobalSumVector(&gbad, 1); - std::cout << GridLogMessage << "DenseCoarseMatrix: fp64 import certificate " - << "max|A64 - A32| = " << gmx + std::cout << GridLogMessage << "DenseCoarseMatrix: inversion-source import certificate " + << "max|A_inv - A_slab| = " << gmx << " non-finite entries " << (int64_t)gbad << std::endl; GRID_ASSERT( gbad == 0 ); GRID_ASSERT( gmx < 1.0e-5 ); S.Resize(nrows, N); - acceleratorCopyToDevice(&h[0], &S.data[0], (uint64_t)nrows*N*sizeof(ComplexD)); + acceleratorCopyToDevice(&h[0], &S.data[0], (uint64_t)nrows*N*sizeof(DenseInverseScalar)); } //////////////////////////////////////////////////////////////////// - // 3c. The inverse: distributed recursive Schur, END-TO-END fp64. - // stencil (ComplexD) -> fp64 rank-major import -> fp64 recursion -> - // ONE terminal rounding into the fp32 apply slab. Everything - // downstream (device residency, split-K apply, VERIFY) is fp32. + // 3c. The inverse: distributed recursive Schur, end to end in the + // inversion precision (DenseInverseScalar, a configure-time choice): + // stencil -> rank-major import -> recursion -> ONE terminal rounding + // into the fp32 apply slab. Everything downstream (device + // residency, split-K apply, VERIFY) is fp32 regardless. //////////////////////////////////////////////////////////////////// template void InvertDense(CoarseOp &Op) @@ -578,7 +599,7 @@ public: } BlockRows S; - ImportDenseFP64(Op, S, g2rm); + ImportDenseForInversion(Op, S, g2rm); //////////////////////////////////////////////////////////////// // The 2D block-cyclic recursion (BlockCyclicSchurInverse). @@ -609,8 +630,8 @@ public: // The single terminal rounding: fp64 inverse -> fp32 apply slab // (row-major, global columns) { - std::vector h((uint64_t)nrows*N); - acceleratorCopyFromDevice(&S.data[0], &h[0], (uint64_t)nrows*N*sizeof(ComplexD)); + std::vector h((uint64_t)nrows*N); + acceleratorCopyFromDevice(&S.data[0], &h[0], (uint64_t)nrows*N*sizeof(DenseInverseScalar)); thread_for(gcol, N, { int64_t jj = g2rm[gcol]; for(int64_t i=0; i - - 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 // needed for Dagger(Yes|No), Inverse(Yes|No) - -#include -#include -#include - -NAMESPACE_BEGIN(Grid); - -// Fine Object == (per site) type of fine field -// nbasis == number of deflation vectors -template -class GeneralCoarsenedMatrix : public SparseMatrixBase > > { -public: - - typedef GeneralCoarsenedMatrix GeneralCoarseOp; - typedef iVector siteVector; - typedef iMatrix siteMatrix; - 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 Lattice FineComplexField; - typedef CoarseVector Field; - //////////////////// - // Data members - //////////////////// - int hermitian; - GridBase * _FineGrid; - GridCartesian * _CoarseGrid; - NonLocalStencilGeometry &geom; - PaddedCell Cell; - GeneralLocalStencil Stencil; - - std::vector _A; - std::vector _Adag; - std::vector MultTemporaries; - - int64_t MultCalls; - double MultFlopsAccum; - double MultUsecAccum; - - /////////////////////// - // Interface - /////////////////////// - ////////////////////////////////////////////////////////////////////////// - // Bilingual accessors: everything a consumer needs to read the operator - // without knowing which of the three coarse classes it holds. The D - // dimensional grid the elements live on, the geometry they are indexed by, - // and one unpadded point at a time (a whole npoint vector is too much - // memory at production nbasis). - ////////////////////////////////////////////////////////////////////////// - GridCartesian * CoarseGridD(void) { return _CoarseGrid; }; - NonLocalStencilGeometry & Geometry(void) { return geom; }; - void ExtractMatrix(int p,CoarseMatrix &A) { A = Cell.Extract(_A[p]); }; - - GridBase * Grid(void) { return _CoarseGrid; }; // this is all the linalg routines need to know - GridBase * FineGrid(void) { return _FineGrid; }; // this is all the linalg routines need to know - GridCartesian * CoarseGrid(void) { return _CoarseGrid; }; // this is all the linalg routines need to know - - /* void ShiftMatrix(RealD shift) - { - int Nd=_FineGrid->Nd(); - Coordinate zero_shift(Nd,0); - for(int p=0;p &A,const CoarseVector &in, CoarseVector &out) - { - RealD tviews=0; RealD ttot=0; RealD tmult=0; RealD texch=0; RealD text=0; RealD ttemps=0; RealD tcopy=0; - RealD tmult2=0; - - ttot=-usecond(); - conformable(CoarseGrid(),in.Grid()); - conformable(in.Grid(),out.Grid()); - out.Checkerboard() = in.Checkerboard(); - CoarseVector tin=in; - - texch-=usecond(); - CoarseVector pin = Cell.ExchangePeriodic(tin); - texch+=usecond(); - - CoarseVector pout(pin.Grid()); - - int npoint = geom.npoint; - typedef LatticeView Aview; - typedef LatticeView Vview; - - const int Nsimd = CComplex::Nsimd(); - - int64_t osites=pin.Grid()->oSites(); - - RealD flops = 1.0* npoint * nbasis * nbasis * 8.0 * osites * CComplex::Nsimd(); - RealD bytes = 1.0*osites*sizeof(siteMatrix)*npoint - + 2.0*osites*sizeof(siteVector)*npoint; - - { - tviews-=usecond(); - autoView( in_v , pin, AcceleratorRead); - autoView( out_v , pout, AcceleratorWriteDiscard); - autoView( Stencil_v , Stencil, AcceleratorRead); - tviews+=usecond(); - - // Static and prereserve to keep UVM region live and not resized across multiple calls - ttemps-=usecond(); - MultTemporaries.resize(npoint,pin.Grid()); - ttemps+=usecond(); - std::vector AcceleratorViewContainer_h; - std::vector AcceleratorVecViewContainer_h; - - tviews-=usecond(); - for(int p=0;p AcceleratorViewContainer; AcceleratorViewContainer.resize(npoint); - static deviceVector AcceleratorVecViewContainer; AcceleratorVecViewContainer.resize(npoint); - - auto Aview_p = &AcceleratorViewContainer[0]; - auto Vview_p = &AcceleratorVecViewContainer[0]; - tcopy-=usecond(); - acceleratorCopyToDevice(&AcceleratorViewContainer_h[0],&AcceleratorViewContainer[0],npoint *sizeof(Aview)); - acceleratorCopyToDevice(&AcceleratorVecViewContainer_h[0],&AcceleratorVecViewContainer[0],npoint *sizeof(Vview)); - tcopy+=usecond(); - - tmult-=usecond(); - accelerator_for(spb, osites*nbasis*npoint, Nsimd, { - typedef decltype(coalescedRead(in_v[0](0))) calcComplex; - int32_t ss = spb/(nbasis*npoint); - int32_t bp = spb%(nbasis*npoint); - int32_t point= bp/nbasis; - int32_t b = bp%nbasis; - auto SE = Stencil_v.GetEntry(point,ss); - auto nbr = coalescedReadGeneralPermute(in_v[SE->_offset],SE->_permute,Nd); - auto res = coalescedRead(Aview_p[point][ss](0,b))*nbr(0); - for(int bb=1;bbgSites() ;bidx++){ - Coordinate bcoor; - CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor); - - for(int p=0;pGlobalDimensions()[mu]; - scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic - } - // Flip to poke/peekLocalSite and not too bad - auto link = peekSite(_A[p],scoor); - int pp = geom.Reverse(p); - pokeSite(adj(link),_Adag[pp],bcoor); - } - } -#else - // Parallel: _Adag[pp](x) = adj( _A[p](x + s_pp) ), pp = Reverse(p), s_pp = -s_p. - // The neighbour fetch reuses the same padded-cell + stencil machinery as Mult, - // reading one matrix element per coalesced access so no whole site matrix - // (230KB at nbasis=60) ever lands on a GPU thread stack (HIP limit 128KB). - // Halo sites compute garbage neighbours; Cell.Extract discards them. - // Must run on the unpadded _A, i.e. before ExchangeCoarseLinks. - const int Nsimd = CComplex::Nsimd(); - for(int p=0;poSites(); - { - autoView( Apad_v , Apad, AcceleratorRead); - autoView( Dpad_v , Dpad, AcceleratorWriteDiscard); - autoView( Stencil_v, Stencil, AcceleratorRead); - accelerator_for(sj, osites*nbasis, Nsimd, { - int32_t ss = sj/nbasis; - int32_t j = sj%nbasis; - auto SE = Stencil_v.GetEntry(pp,ss); - for(int i=0;i_offset](i,j),SE->_permute,Nd); - coalescedWrite(Dpad_v[ss](j,i),conjugate(z)); - } - }); - } - _Adag[pp] = Cell.Extract(Dpad); - } -#endif - } - ///////////////////////////////////////////////////////////// - // - // A) Only reduced flops option is to use a padded cell of depth 4 - // and apply MpcDagMpc in the padded cell. - // - // Makes for ONE application of MpcDagMpc per vector instead of 30 or 80. - // With the effective cell size around (B+8)^4 perhaps 12^4/4^4 ratio - // Cost is 81x more, same as stencil size. - // - // But: can eliminate comms and do as local dirichlet. - // - // Local exchange gauge field once. - // Apply to all vectors, local only computation. - // Must exchange ghost subcells in reverse process of PaddedCell to take inner products - // - // B) Can reduce cost: pad by 1, apply Deo (4^4+6^4+8^4+8^4 )/ (4x 4^4) - // pad by 2, apply Doe - // pad by 3, apply Deo - // then break out 8x directions; cost is ~10x MpcDagMpc per vector - // - // => almost factor of 10 in setup cost, excluding data rearrangement - // - // Intermediates -- ignore the corner terms, leave approximate and force Hermitian - // Intermediates -- pad by 2 and apply 1+8+24 = 33 times. - ///////////////////////////////////////////////////////////// - - ////////////////////////////////////////////////////////// - // BFM HDCG style approach: Solve a system of equations to get Aij - ////////////////////////////////////////////////////////// - /* - * Here, k,l index which possible shift within the 3^Nd "ball" connected by MdagM. - * - * conj(phases[block]) proj[k][ block*Nvec+j ] = \sum_ball e^{i q_k . delta} < phi_{block,j} | MdagM | phi_{(block+delta),i} > - * = \sum_ball e^{iqk.delta} A_ji - * - * Must invert matrix M_k,l = e^[i q_k . delta_l] - * - * Where q_k = delta_k . (2*M_PI/global_nb[mu]) - */ -#if 0 - void CoarsenOperator(LinearOperatorBase > &linop, - Aggregation & Subspace) - { - std::cout << GridLogMessage<< "GeneralCoarsenMatrix "<< std::endl; - GridBase *grid = FineGrid(); - - RealD tproj=0.0; - RealD teigen=0.0; - RealD tmat=0.0; - RealD tphase=0.0; - RealD tinv=0.0; - - ///////////////////////////////////////////////////////////// - // Orthogonalise the subblocks over the basis - ///////////////////////////////////////////////////////////// - CoarseScalar InnerProd(CoarseGrid()); - blockOrthogonalise(InnerProd,Subspace.subspace); - - const int npoint = geom.npoint; - - Coordinate clatt = CoarseGrid()->GlobalDimensions(); - int Nd = CoarseGrid()->Nd(); - - /* - * Here, k,l index which possible momentum/shift within the N-points connected by MdagM. - * Matrix index i is mapped to this shift via - * geom.shifts[i] - * - * conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] - * = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > - * = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} - * = M_{kl} A_ji^{b.b+l} - * - * Must assemble and invert matrix M_k,l = e^[i q_k . delta_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} - */ - teigen-=usecond(); - Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint); - Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint); - ComplexD ci(0.0,1.0); - for(int k=0;k ComputeProj(npoint,CoarseGrid()); - std::vector FT(npoint,CoarseGrid()); - for(int i=0;ioSites(); - autoView( A_v , _A[k], AcceleratorWrite); - autoView( FT_v , FT[k], AcceleratorRead); - accelerator_for(sss, osites, nbasis, { -#ifdef GRID_SIMT - int j = acceleratorSIMTlane(nbasis); - A_v[sss](i,j) = FT_v[sss](j); -#else - // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor - // index writes ONLY j=0 and silently drops the other nbasis-1 - // columns (caught by Test_schur_dense_coarse import certificate, - // 2026-08-14). Loop explicitly. - for(int j=0;j > &linop, - Aggregation & Subspace) - { - CoarsenOperator(linop,Subspace,Subspace); - } - ////////////////////////////////////////////////////////////////////// - // Petrov - Galerkin projection of matrix - ////////////////////////////////////////////////////////////////////// - void CoarsenOperator(LinearOperatorBase > &linop, - Aggregation & U, - Aggregation & V) - { - std::cout << GridLogMessage<< "GeneralCoarsenMatrix "<< std::endl; - GridBase *grid = FineGrid(); - - RealD tproj=0.0; - RealD teigen=0.0; - RealD tmat=0.0; - RealD tphase=0.0; - RealD tphaseBZ=0.0; - RealD tinv=0.0; - - ///////////////////////////////////////////////////////////// - // Orthogonalise the subblocks over the basis - ///////////////////////////////////////////////////////////// - CoarseScalar InnerProd(CoarseGrid()); - blockOrthogonalise(InnerProd,V.subspace); - blockOrthogonalise(InnerProd,U.subspace); - - const int npoint = geom.npoint; - - Coordinate clatt = CoarseGrid()->GlobalDimensions(); - int Nd = CoarseGrid()->Nd(); - - /* - * Here, k,l index which possible momentum/shift within the N-points connected by MdagM. - * Matrix index i is mapped to this shift via - * geom.shifts[i] - * - * conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] - * = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > - * = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} - * = M_{kl} A_ji^{b.b+l} - * - * Must assemble and invert matrix M_k,l = e^[i q_k . delta_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} - */ - teigen-=usecond(); - Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint); - Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint); - ComplexD ci(0.0,1.0); - for(int k=0;k phaF(npoint,grid); - std::vector pha(npoint,CoarseGrid()); - - typedef typename CComplex::scalar_type SComplex; - FineComplexField one(grid); one=SComplex(1.0); - FineComplexField zz(grid); zz = Zero(); - tphase=-usecond(); - for(int p=0;p Projector; - Projector.Allocate(nbasis, grid, CoarseGrid()); - Projector.ImportBasis(U.subspace); - - std::vector phaV_batch(npoint, grid); - std::vector MphaV_batch(npoint, grid); - std::vector proj_batch(npoint, CoarseGrid()); - std::vector ComputeProj(npoint, CoarseGrid()); - std::vector FT(npoint, CoarseGrid()); - - // Pre-allocate BLAS_F and BLAS_C to avoid repeated hipMalloc/hipFree of - // ~5.6 GB per blockProject call, which hangs on ROCm for large allocations. - Projector.BLAS_F.resize(Projector.fine_vol * Projector.words * npoint); - Projector.BLAS_C.resize(Projector.coarse_vol * nbasis * npoint); - - for(int i=0;ioSites(); - autoView( A_v , _A[k], AcceleratorWrite); - autoView( FT_v , FT[k], AcceleratorRead); - accelerator_for(sss, osites, nbasis, { -#ifdef GRID_SIMT - int j = acceleratorSIMTlane(nbasis); - A_v[sss](i,j) = FT_v[sss](j); -#else - // CPU build: acceleratorSIMTlane()==0 -- an un-looped SIMT tensor - // index writes ONLY j=0 and silently drops the other nbasis-1 - // columns (caught by Test_schur_dense_coarse import certificate, - // 2026-08-14). Loop explicitly. - for(int j=0;j &out){assert(0);}; -}; - - - -NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h b/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h deleted file mode 100644 index 292f6f035..000000000 --- a/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h +++ /dev/null @@ -1,738 +0,0 @@ -/************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./lib/algorithms/GeneralCoarsenedMatrixMultiRHS.h - - Copyright (C) 2015 - -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); - - -// Fine Object == (per site) type of fine field -// nbasis == number of deflation vectors -template -class MultiGeneralCoarsenedMatrix : public SparseMatrixBase > > { -public: - typedef typename CComplex::scalar_object SComplex; - typedef GeneralCoarsenedMatrix GeneralCoarseOp; - typedef MultiGeneralCoarsenedMatrix 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 Lattice FineComplexField; - typedef CoarseVector Field; - - //////////////////// - // Data members - //////////////////// - GridCartesian * _CoarseGridMulti; - NonLocalStencilGeometry geom; - NonLocalStencilGeometry geom_srhs; - PaddedCell Cell; - GeneralLocalStencil Stencil; - - deviceVector BLAS_B; - deviceVector BLAS_C; - std::vector > BLAS_A; - - std::vector > BLAS_AP; - std::vector > BLAS_BP; - deviceVector BLAS_CP; - - /////////////////////// - // Interface - /////////////////////// - GridBase * Grid(void) { return _CoarseGridMulti; }; // this is all the linalg routines need to know - GridCartesian * CoarseGrid(void) { return _CoarseGridMulti; }; // this is all the linalg routines need to know - - ////////////////////////////////////////////////////////////////////////// - // Bilingual accessors, matching GeneralCoarsenedMatrix. 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]); }; - - // 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]); - } - void CopyMatrix (GeneralCoarseOp &_Op) - { - for(int p=0;plSites(); - int32_t unpadded_sites = CoarseGridMulti->lSites(); - - int32_t nrhs = CoarseGridMulti->FullDimensions()[0]; // # RHS - int32_t orhs = nrhs/CComplex::Nsimd(); - - padded_sites = padded_sites/nrhs; - unpadded_sites = unpadded_sites/nrhs; - - ///////////////////////////////////////////////// - // Device data vector storage - ///////////////////////////////////////////////// - BLAS_A.resize(geom.npoint); - for(int p=0;p lSite - 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 > &linop, - Aggregation & Subspace, - GridBase *CoarseGrid) - { -#if 0 - std::cout << GridLogMessage<< "GeneralCoarsenMatrixMrhs "<< std::endl; - - GridBase *grid = Subspace.FineGrid; - - ///////////////////////////////////////////////////////////// - // Orthogonalise the subblocks over the basis - ///////////////////////////////////////////////////////////// - CoarseScalar InnerProd(CoarseGrid); - blockOrthogonalise(InnerProd,Subspace.subspace); - - const int npoint = geom_srhs.npoint; - - Coordinate clatt = CoarseGrid->GlobalDimensions(); - int Nd = CoarseGrid->Nd(); - /* - * Here, k,l index which possible momentum/shift within the N-points connected by MdagM. - * Matrix index i is mapped to this shift via - * geom.shifts[i] - * - * conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] - * = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > - * = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} - * = M_{kl} A_ji^{b.b+l} - * - * Must assemble and invert matrix M_k,l = e^[i q_k . delta_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} - */ - Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint); - Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint); - ComplexD ci(0.0,1.0); - for(int k=0;k phaF(npoint,grid); - std::vector pha(npoint,CoarseGrid); - - CoarseVector coarseInner(CoarseGrid); - - typedef typename CComplex::scalar_type SComplex; - FineComplexField one(grid); one=SComplex(1.0); - FineComplexField zz(grid); zz = Zero(); - for(int p=0;p _A; - _A.resize(geom_srhs.npoint,CoarseGrid); - - std::vector ComputeProj(npoint,CoarseGrid); - CoarseVector FT(CoarseGrid); - for(int i=0;ioSites(); - autoView( A_v , _A[k], AcceleratorWrite); - autoView( FT_v , FT, AcceleratorRead); - accelerator_for(sss, osites, 1, { - for(int j=0;j > Projector; - Projector.Allocate(nbasis,grid,CoarseGrid); - Projector.ImportBasis(Subspace.subspace); - - const int npoint = geom_srhs.npoint; - - Coordinate clatt = CoarseGrid->GlobalDimensions(); - int Nd = CoarseGrid->Nd(); - /* - * Here, k,l index which possible momentum/shift within the N-points connected by MdagM. - * Matrix index i is mapped to this shift via - * geom.shifts[i] - * - * conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block] - * = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > - * = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} - * = M_{kl} A_ji^{b.b+l} - * - * Must assemble and invert matrix M_k,l = e^[i q_k . delta_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} - */ - Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint); - Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint); - ComplexD ci(0.0,1.0); - for(int k=0;k phaF(npoint,grid); - std::vector pha(npoint,CoarseGrid); - - CoarseVector coarseInner(CoarseGrid); - - tphase=-usecond(); - typedef typename CComplex::scalar_type SComplex; - FineComplexField one(grid); one=SComplex(1.0); - FineComplexField zz(grid); zz = Zero(); - for(int p=0;p _A; - _A.resize(geom_srhs.npoint,CoarseGrid); - - // Count use small chunks than npoint == 81 and save memory - int batch = 9; - std::vector _MphaV(batch,grid); - std::vector TmpProj(batch,CoarseGrid); - - std::vector ComputeProj(npoint,CoarseGrid); - CoarseVector FT(CoarseGrid); - for(int i=0;ioSites(); - autoView( A_v , _A[k], AcceleratorWrite); - autoView( FT_v , FT, AcceleratorRead); - accelerator_for(sss, osites, 1, { - for(int j=0;jM(in,out); - } - void M (const CoarseVector &in, CoarseVector &out) - { - // std::cout << GridLogMessage << "New Mrhs coarse"< 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(); - GridtoBLAS(pin,BLAS_B); - t_GtoB+=usecond(); - - GridBLAS BLAS; - - t_mult=-usecond(); - for(int p=0;p &out){assert(0);}; -}; - -NAMESPACE_END(Grid); diff --git a/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h b/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h index 734ba2e35..0af4b0ca7 100644 --- a/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h +++ b/Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h @@ -37,7 +37,9 @@ template class MultiGeneralCoarsenedOperatorV2 : public SparseMatrixBase > > { public: typedef typename CComplex::scalar_object SComplex; - typedef GeneralCoarsenedMatrix GeneralCoarseOp; + // 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; @@ -82,9 +84,80 @@ public: deviceVector BLAS_C; std::vector > BLAS_A; - std::vector > BLAS_AP; - std::vector > BLAS_BP; - deviceVector BLAS_CP; + 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 @@ -115,6 +188,15 @@ public: // 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]); }; @@ -130,30 +212,7 @@ public: GRID_ASSERT(A.size()==geom_srhs.npoint); BLAStoGrid(A[p],BLAS_A[p]); } - void CopyMatrix (GeneralCoarseOp &_Op) - { - 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() @@ -209,6 +282,11 @@ public: 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 lSite, D dim nbr = nbr*nrhs; // D -> D+1, rhs innermost GRID_ASSERT(nbrGlobalDimensions(); + 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); @@ -424,13 +558,21 @@ public: ComplexD phase(0.0,0.0); for(int mu=0;mu svd(Mkl); + RealD cond = svd.singularValues()(0)/svd.singularValues()(npoint-1); + std::cout << GridLogMessage << "CoarsenOperator: probe momenta "<GlobalDimensions(); int Nd = CoarseGrid->Nd(); + Coordinate K; CoarsenMomenta(CoarseGrid,K); ComplexD ci(0.0,1.0); - typedef typename CComplex::scalar_type SComplex; - FineComplexField one(grid); one=SComplex(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); @@ -492,9 +637,9 @@ public: for(int mu=0;mu > &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; @@ -585,7 +745,7 @@ public: BlockComplexField InnerProd(&BlockGrid); blockOrthogonalise(InnerProd,Subspace); - MultiRHSBlockProject > Projector; + // the caller owns it; import the basis we have just orthonormalised Projector.Allocate(nbasis,grid,CoarseGrid); Projector.ImportBasis(Subspace); @@ -674,6 +834,16 @@ public: 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; @@ -690,7 +860,7 @@ public: BlockComplexField InnerProd(&BlockGrid); blockOrthogonalise(InnerProd,Subspace); - MultiRHSBlockProject > Projector; + // the caller owns it; import the basis we have just orthonormalised Projector.Allocate(nbasis,grid,CoarseGrid); Projector.ImportBasis(Subspace); @@ -787,16 +957,15 @@ public: GRID_TRACE("CoarseV2Mult"); t_tot=-usecond(); - CoarseVector tin=in; t_exch=-usecond(); - // lambda scope so the roctx range covers exactly the exchange; the - // PaddedCellFwd/BwdMPI markers inside it then nest properly. + // The exchange takes its input by const reference, so it reads the + // caller's field directly. lambda scope so the roctx range covers + // exactly the exchange; the PaddedCellFwd/BwdMPI markers inside it then + // nest properly. CoarseVector pin = [&](){ GRID_TRACE("CoarseV2Exchange"); - return CellMulti->ExchangePeriodic(tin); }(); //padded input + return CellMulti->ExchangePeriodic(in); }(); //padded input t_exch+=usecond(); - CoarseVector pout(pin.Grid()); - int npoint = geom.npoint; typedef calcMatrix* Aview; typedef LatticeView Vview; @@ -825,22 +994,64 @@ public: t_mult=-usecond(); { GRID_TRACE("CoarseV2StencilGEMM"); - 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 #pragma once #include +#include NAMESPACE_BEGIN(Grid); -////////////////////////////////////////////////////////////////////// -// mrhs LinearFunction interface: vector-of-fields in, vector out. -// The outer level carries mrhs as std::vector; below it mrhs -// is PACKED into a single D+1 field (rhs = dim 0) and the coarse -// classes are plain LinearFunctions on that. -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; ////////////////////////////////////////////////////////////////////// // Single-polynomial mrhs PGCR: one step length and one set of @@ -54,18 +44,21 @@ public: RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; int ZeroGuess = 0; int FirstCycle = 0; std::string name = "Level 1"; + // Trace range names carry the instance name, so a profile separates the + // solvers that otherwise nest as one shared range. + std::string trace_op = "MrhsPGCR::vOp", trace_orthog = "MrhsPGCR orthog"; 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; } + 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("MrhsPGCR::vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[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; @@ -108,7 +101,7 @@ public: vOp(p[peri_kp],q[peri_kp]); int northog=((kp)>(mmax-1))?(mmax-1):(kp); { - GRID_TRACE("MrhsPGCR orthog"); + GRID_TRACE(trace_orthog.c_str()); // Classical Gram-Schmidt: all coefficients against the UN-updated new q // (independent, batchable), then apply. Complex coefficient: the // operator is non-Hermitian, real() alone left q's non-orthogonal. @@ -203,19 +196,23 @@ public: // is wired by the composer to PVdagMLinearOperator::SloppyComms (a // no-op by default), replacing the file-scope global the example used. ////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { +// Projector_t defaults to a transfer operator whose STORE matches FineField, +// which is the case whenever the coarse sector and the fine level share a +// precision; pass it explicitly when they differ. +template > +class MrhsTwoLevelMG : public MrhsPreconditioner { public: typedef MrhsCoarseVector CoarseVector; LinearOperatorBase &_FineOperator; FineSmoother &_PostSmoother; - MultiRHSBlockProject &_Projector; + Projector_t &_Projector; // store precision need not match FineField LinearFunction &_CoarseSolve; GridBase *_CoarseGrid, *_CoarseGridMrhs; std::function SetSloppy = [](int){}; int SloppyComms = 0; // value passed to SetSloppy on entry MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, - MultiRHSBlockProject &Projector, LinearFunction &CoarseSolve, + Projector_t &Projector, LinearFunction &CoarseSolve, GridBase *CoarseGrid, GridBase *CoarseGridMrhs) : _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve), _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){} @@ -253,4 +250,52 @@ public: } }; +////////////////////////////////////////////////////////////////////// +// The fp64/fp32 seam of the solve chain. The outer Krylov hands fp64 +// residuals to its preconditioner; this adapter converts them to fp32, +// runs an fp32 preconditioner (the whole V-cycle), and converts the +// correction back. Two precisionChange per outer step per rhs, on +// workspaces built once. The outer operator never sees fp32. +// +// The fp32 scratch is allocated ON FIRST USE, not in the constructor: the +// solver builds this seam whichever precision is selected, and at Nrhs 12 +// on a 48^3x96 Ls=24 rank the scratch is 2 GB of device memory. A seam +// that is never called must cost nothing. Same idiom as the GCR history +// vectors: re-made only if the grid or the rhs count changes. +////////////////////////////////////////////////////////////////////// +template +class MrhsMixedPrecPreconditioner : public MrhsPreconditioner { +public: + MrhsPreconditioner &_Inner; + precisionChangeWorkspace _ws_d2f; // out fp32, in fp64 + precisionChangeWorkspace _ws_f2d; // out fp64, in fp32 + GridBase *_gridF; + std::vector _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){ + GRID_TRACE("MGPrecisionSeam"); + int nrhs=in.size(); Scratch(nrhs); + for(int r=0;r &x, std::vector &src){ + GRID_TRACE("MGPrecisionSeamVstart"); + int nrhs=src.size(); Scratch(nrhs); + for(int r=0;r #include #include #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 +#include +#include #include #include #include // PVdagMOperators.h / MrhsMultiGrid.h / PVdagMMultiGrid.h / -// DenseCoarseMatrix.h are NOT in this umbrella: consumers of the PVdagM -// chain include PVdagMMultiGrid.h explicitly (it pulls the dense stack -// and BLAS). +// DenseCoarseMatrix.h / MultiGridIO.h / HDCGMultiGrid.h are NOT in this +// umbrella: consumers include PVdagMMultiGrid.h or HDCGMultiGrid.h +// explicitly (they pull the dense stack, BLAS, and the scidac I/O, which +// is declared after Algorithms.h in Grid.h). Keeping MrhsMultiGrid.h out +// of the Algorithms.h chain also keeps its class names away from the +// pre-2026 drivers that define their own copies. diff --git a/Grid/algorithms/multigrid/PVdagMMultiGrid.h b/Grid/algorithms/multigrid/PVdagMMultiGrid.h index 30d4bb737..1e646044d 100644 --- a/Grid/algorithms/multigrid/PVdagMMultiGrid.h +++ b/Grid/algorithms/multigrid/PVdagMMultiGrid.h @@ -20,6 +20,7 @@ Author: Peter Boyle #pragma once #include +#include #include #include @@ -48,36 +49,13 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// -// Subspace I/O: bare-vector scidac records. Loaded vectors are RAW -- -// deliberately NOT re-orthogonalised: CoarsenOperator block -// orthonormalises in place, and projecting a block-orthonormal vector -// onto its own block-orthonormalised aggregation gives e_k with the -// near-null content silently gone. GramGuard below catches that. +// Subspace I/O is in MultiGridIO.h (shared with the HDCG chain). +// Loaded vectors are RAW -- deliberately NOT re-orthogonalised: +// CoarsenOperator block orthonormalises in place, and projecting a +// block-orthonormal vector onto its own block-orthonormalised +// aggregation gives e_k with the near-null content silently gone. +// GramGuard below catches that. ////////////////////////////////////////////////////////////////////// -template -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 -} ////////////////////////////////////////////////////////////////////// // || - I||_F over a set of coarse vectors. Small means the raw @@ -146,14 +124,14 @@ public: clatt.resize(4); cclatt.resize(4); for(int d=0;d<4;d++){ - GRID_ASSERT( fdims[d+1] % P.Block[d] == 0 ); - clatt[d] = fdims[d+1] / P.Block[d]; + GRID_ASSERT( fdims[d+1] % P.Block1[d] == 0 ); + clatt[d] = fdims[d+1] / P.Block1[d]; } for(int d=0;d<4;d++){ GRID_ASSERT( clatt[d] % P.Block2[d] == 0 ); cclatt[d] = clatt[d] / P.Block2[d]; } - std::cout << GridLogMessage << "MGCoarseGrids: Block " << P.Block << " coarse lattice " << clatt << std::endl; + std::cout << GridLogMessage << "MGCoarseGrids: Block1 " << P.Block1 << " coarse lattice " << clatt << std::endl; std::cout << GridLogMessage << "MGCoarseGrids: Block2 " << P.Block2 << " coarse-coarse lattice " << cclatt << std::endl; Coordinate c5latt({1,clatt[0],clatt[1],clatt[2],clatt[3]}); @@ -181,6 +159,41 @@ public: } }; +////////////////////////////////////////////////////////////////////// +// The fp32 fine grids: same lattice and decomposition as the fp64 fine +// grid, vComplexF SIMD layout. For the fp32 fine level inside the +// preconditioner (the fp32 fermion operators and the fp32 transfer +// operator borrow these). Owned here; declare BEFORE the coarsening. +////////////////////////////////////////////////////////////////////// +class MGFineGridsF { +public: + GridCartesian *UGridF; + GridRedBlackCartesian *UrbGridF; + GridCartesian *FGridF; + GridRedBlackCartesian *FrbGridF; + + MGFineGridsF(GridCartesian *FGrid) + { + Coordinate fdims = FGrid->FullDimensions(); // {Ls, x,y,z,t} + Coordinate fmpi = FGrid->_processors; + GRID_ASSERT( fdims.size() == 5 ); + int Ls = fdims[0]; + Coordinate latt4({fdims[1],fdims[2],fdims[3],fdims[4]}); + Coordinate mpi4 ({fmpi[1], fmpi[2], fmpi[3], fmpi[4]}); + UGridF = SpaceTimeGrid::makeFourDimGrid(latt4,GridDefaultSimd(Nd,vComplexF::Nsimd()),mpi4); + UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF); + FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF); + FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF); + } + ~MGFineGridsF() + { + delete FrbGridF; + delete FGridF; + delete UrbGridF; + delete UGridF; + } +}; + ////////////////////////////////////////////////////////////////////// // ALL the coarsening information for one gauge configuration. // @@ -202,19 +215,41 @@ public: typedef MultiGeneralCoarsenedOperatorV2 CoarseOperator; typedef typename CoarseOperator::CoarseVector CoarseVector; typedef typename CoarseVector::vector_object CoarseSiteObj; - typedef iScalar CComplex2; // coarsening deepens the nest by one iScalar - typedef MultiGeneralCoarsenedOperatorV2 CoarseCoarseOperator; + // 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 typename CoarseCoarseOperator::CoarseVector CoarseCoarseVector; - typedef DenseCoarseMatrix DenseBottom; + typedef DenseCoarseMatrix DenseBottom; + // The fp32 fine field, derived from the fp64 one: the fp32 fine level of + // the preconditioner projects through the SAME basis in fp32 layout. + typedef typename GridTypeMapper::SinglePrecision FobjF; + typedef Lattice FineFieldF; + + ////////////////////////////////////////////////////////////////////// + // ONE fine transfer operator. Its STORE follows the coarse sector -- + // the coarse space is what it feeds -- while its import and export + // accept EITHER fine precision, because they are already a layout + // transformation and a scalar conversion inside one costs nothing. + // So an fp64 outer level and an fp32 V-cycle share this single object, + // and no second basis store, handover or retained basis is needed. + ////////////////////////////////////////////////////////////////////// + static const bool CoarseIsSingle = + ( sizeof(typename GridTypeMapper::scalar_type) == sizeof(ComplexF) ); + typedef typename std::conditional::type ProjectorField; + typedef MultiRHSBlockProject ProjectorL1_t; + typedef MultiRHSBlockProject ProjectorL2_t; MGCoarseGrids &Grids; // borrowed + MGFineGridsF &GridsF; // borrowed MGSetupParams Params; NextToNearestStencilGeometry5D geom; NextToNearestStencilGeometry5D geom2; CoarseOperator CoarseOpPV; CoarseCoarseOperator CoarseOpL2; - MultiRHSBlockProject MrhsProjector; - MultiRHSBlockProject MrhsProjectorL2; + ProjectorL1_t MrhsProjectorL1; // store follows the coarse sector + ProjectorL2_t MrhsProjectorL2; DenseBottom *DenseCC; std::vector rawNull; // RAW fine near-null basis std::vector rawPsi; // RAW coarse near-null basis @@ -222,8 +257,9 @@ public: GridCartesian *CCMrhs; int nrhs; - PVdagMMultiGridCoarsening(MGCoarseGrids &_Grids, const MGSetupParams &P) + PVdagMMultiGridCoarsening(MGCoarseGrids &_Grids, MGFineGridsF &_GridsF, const MGSetupParams &P) : Grids(_Grids), + GridsF(_GridsF), Params(P), geom (_Grids.Coarse5d), geom2(_Grids.CoarseCoarse5d), @@ -289,7 +325,7 @@ public: // halos follow the FineSloppyComms policy, restored EXACT on exit). //////////////////////////////////////////////////////////////////// template - void Coarsen(PVdagMOp &FineOp) + void Coarsen(PVdagMOp &FineOp, int retain_basis=0) { GRID_ASSERT( rawNull.size() == nbasis ); @@ -301,18 +337,22 @@ public: std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L1 CoarsenOperator, batch " << Grids.batch << std::endl; + // CoarsenOperator orthonormalises sub in place and imports it into OUR + // transfer operator, which it then uses. One basis store, imported once. + // The projector's grid carries the STORE's SIMD layout, which is the fp32 + // fine grid when the coarse sector is fp32; fp64 vectors convert on import. + MrhsProjectorL1.Allocate(nbasis, + CoarseIsSingle ? (GridBase *)GridsF.FGridF : (GridBase *)Grids.FGrid, + Grids.Coarse5d); FineOp.SloppyComms(Params.FineSloppyComms); - CoarseOpPV.CoarsenOperator(FineOp,sub,Grids.Coarse5d,Grids.batch); + CoarseOpPV.CoarsenOperator(FineOp,sub,Grids.Coarse5d,Grids.batch,MrhsProjectorL1); FineOp.SloppyComms(0); - - // Transfer operators from the orthonormalised basis; then the - // Galerkin images of the RAW basis define the L2 null space. - MrhsProjector.Allocate(nbasis,Grids.FGrid,Grids.Coarse5d); - MrhsProjector.ImportBasis(sub); + // The Lattice copy of the basis has served its purpose: the transfer + // operator holds it in BLAS layout from here on. sub.clear(); sub.shrink_to_fit(); std::vector psi(nbasis,Grids.Coarse5d); - MrhsProjector.blockProject(rawNull,psi); + MrhsProjectorL1.blockProject(rawNull,psi); GramGuard("psi_coarse",psi,Grids.Coarse5d); rawPsi.clear(); @@ -324,16 +364,25 @@ public: NonHermitianLinearOperator LinOpCoarse(CoarseOpPV); std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L2 CoarsenOperator, batch " << Grids.batch << std::endl; - CoarseOpL2.CoarsenOperator(LinOpCoarse,Grids.CoarseBatch,psi,Grids.CoarseCoarse5d); - MrhsProjectorL2.Allocate(nbasis,Grids.Coarse5d,Grids.CoarseCoarse5d); - MrhsProjectorL2.ImportBasis(psi); // now block orthonormal + CoarseOpL2.CoarsenOperator(LinOpCoarse,Grids.CoarseBatch,psi,Grids.CoarseCoarse5d,MrhsProjectorL2); { std::vector psi_cc(nbasis,Grids.CoarseCoarse5d); MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in GramGuard("psi_cc",psi_cc,Grids.CoarseCoarse5d); } + // The RAW basis has done its work (the L2 null space above was its last + // use) and is the single largest resident block of the setup, so it is + // freed. retain_basis keeps it for a caller that coarsens again from the + // same basis; Setup.RetainSubspace keeps it for the object's whole life, + // which is what a fixed-basis rebuild on a changed gauge field (HMC) + // needs. A valence solve never rebuilds. + if ( !retain_basis && !Params.RetainSubspace ) DiscardBasis(); + MrhsProjectorL1.ReleaseScratch(); + MrhsProjectorL2.ReleaseScratch(); + nrhs = -1; // both operators left at the batch grids + ReportDeviceFootprint("after Coarsen"); } //////////////////////////////////////////////////////////////////// @@ -373,6 +422,7 @@ public: CoarseOpL2.ReleaseGrid(); // let go before the grid dies delete CoarseCoarseOne; nrhs = -1; + ReportDeviceFootprint("after BuildDenseBottom"); } //////////////////////////////////////////////////////////////////// @@ -398,6 +448,66 @@ public: std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: operators at Nrhs " << nr << std::endl; } + //////////////////////////////////////////////////////////////////// + // Galerkin certificate for the L1 coarsening: ||A_c x - P^dag A P x|| + // over ||A_c x|| on a random coarse vector, through the SAME transfer + // operators the solve uses. The coarse operator is compared with its + // own definition, so this reads the ROUNDING level of the coarsening + // (fp64 coarse ~1e-13, fp32 coarse ~1e-6) independently of how the + // solve converges. Drives the fine operator with EXACT halos; at more + // than one rank the sloppy-halo coarsening error is included in the + // reading. Nrhs 1; the caller's Nrhs is restored by its own SetNrhs. + //////////////////////////////////////////////////////////////////// + template + RealD CertifyCoarsening(PVdagMOp &FineOp) + { + SetNrhs(1); + CoarseVector xc(CMrhs), Acx(CMrhs), PtAPx(CMrhs); + GridParallelRNG cRNG(CMrhs); cRNG.SeedFixedIntegers(std::vector({21,22,23,24})); + random(cRNG,xc); + + CoarseOpPV.M(xc,Acx); // A_c x + + std::vector Px(1,Grids.FGrid), APx(1,Grids.FGrid); + MrhsProjectorL1.blockPromote(Px,xc); // P x + FineOp.SloppyComms(0); + FineOp.Op(Px[0],APx[0]); // A P x + MrhsProjectorL1.blockProject(APx,PtAPx); // P^dag A P x + + CoarseVector d(CMrhs); d = Acx - PtAPx; + RealD rel = std::sqrt(norm2(d)/norm2(Acx)); + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening: L1 GALERKIN CERTIFICATE " + << "||A_c x - P^dag A P x||/||A_c x|| = " << rel << std::endl; + return rel; + } + + + //////////////////////////////////////////////////////////////////// + // The device memory this coarsening holds that the memory manager + // CANNOT evict: BLAS stores, which are plain device allocations. + // Lattice fields are omitted on purpose -- they are evictable, so they + // cost eviction traffic, not failure. What is reported is what has to + // fit alongside the manager's cache, the comms buffers and the shm + // segment, which is the budget an out-of-memory run actually breaks. + //////////////////////////////////////////////////////////////////// + void ReportDeviceFootprint(const std::string &when) + { + uint64_t p1 = MrhsProjectorL1.DeviceBytes(); + uint64_t p2 = MrhsProjectorL2.DeviceBytes(); + uint64_t o1 = CoarseOpPV.DeviceBytes(); + uint64_t o2 = CoarseOpL2.DeviceBytes(); + uint64_t dn = DenseCC ? DenseCC->DeviceBytes() : 0; + uint64_t tot = p1+p2+o1+o2+dn; + const double G = 1024.*1024.*1024.; + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening (" << when << "): resident device memory/rank " + << tot/G << " GB = transfer L1 " << p1/G << " + L2 " << p2/G + << " + operator L1 " << o1/G << " + L2 " << o2/G + << " + dense " << dn/G << " GB" << std::endl; + std::cout << GridLogMessage << "PVdagMMultiGridCoarsening (" << when << "): not evictable; must fit alongside the " + << MemoryManager::DeviceMaxBytes/G + << " GB manager cache, the comms buffers and the shm segment" << std::endl; + } + //////////////////////////////////////////////////////////////////// // Free the retained RAW bases (valence use: coarsening is final for // this configuration and the fine basis is ~GB-scale host memory). @@ -420,16 +530,18 @@ public: // (its applications define what "converged" means), asserted by the // FINAL true-residual report in Solve. ////////////////////////////////////////////////////////////////////// -template +template class PVdagMMultiGridSolver { public: typedef typename Coarsening::FineField FineField; + typedef typename Coarsening::FineFieldF FineFieldF; typedef typename Coarsening::CoarseOperator CoarseOperator; typedef typename Coarsening::CoarseVector CoarseVector; typedef typename Coarsening::CoarseCoarseOperator CoarseCoarseOperator; typedef typename Coarsening::CoarseCoarseVector CoarseCoarseVector; typedef typename Coarsening::DenseBottom DenseBottom; typedef PrecGeneralisedConjugateResidualNonHermitian FineSmoother_t; + typedef PrecGeneralisedConjugateResidualNonHermitian FineSmootherF_t; typedef PrecGeneralisedConjugateResidualNonHermitian CoarseKrylov_t; Coarsening &C; @@ -447,11 +559,20 @@ public: CoarseKrylov_t CoarseSmootherGCR; MrhsCoarseThreeLevelPrec L2to3Precon; CoarseKrylov_t L2PGCR; + // The fp64 fine level of the preconditioner FineSmoother_t SmootherGCR; - MrhsTwoLevelMG ThreeLevelPrecon; + MrhsTwoLevelMG ThreeLevelPrecon; + // The fp32 fine level of the preconditioner, behind the fp64/fp32 seam. + // Both share the coarse chain (L2PGCR); the outer Krylov picks one. + PVdagMLinearOperator PVdagMF; + ShiftedPVdagMLinearOperator ShiftedPVdagMF; + TrivialPrecon simple_fine_f; + FineSmootherF_t SmootherGCRF; + MrhsTwoLevelMG ThreeLevelPreconF; + MrhsMixedPrecPreconditioner PreconSeam; MrhsPGCRNonHermitian L1PGCR; - PVdagMMultiGridSolver(Matrix &Ddwf, Matrix &Dpv, + PVdagMMultiGridSolver(Matrix &Ddwf, Matrix &Dpv, MatrixF &DdwfF, MatrixF &DpvF, Coarsening &_C, const PVdagMMultiGridParams &P, int nr) : C(_C), Params(P), nrhs(nr), _regrid((C.SetNrhs(nr),0)), // members below capture C.CMrhs/C.CCMrhs @@ -468,11 +589,20 @@ public: L2PGCR(P.CoarseSolver.Tol,P.CoarseSolver.Order/16,LinOpC,L2to3Precon,P.CoarseSolver.Mmax,16), // Fine smoother: one restart of nstep GCR steps, tolerance 0 = fixed work. SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,P.FineSmoother.Mmax,P.FineSmoother.Nstep), - ThreeLevelPrecon(PVdagM,SmootherGCR,C.MrhsProjector,L2PGCR,C.Grids.Coarse5d,C.CMrhs), - L1PGCR(P.Outer.Tol,P.Outer.MaxIterations,PVdagM,ThreeLevelPrecon,P.Outer.Mmax,P.Outer.Nstep) + ThreeLevelPrecon(PVdagM,SmootherGCR,C.MrhsProjectorL1,L2PGCR,C.Grids.Coarse5d,C.CMrhs), + PVdagMF(DdwfF,DpvF), + ShiftedPVdagMF(P.FineSmoother.Shift,DdwfF,DpvF), + SmootherGCRF(0.0,1,ShiftedPVdagMF,simple_fine_f,P.FineSmoother.Mmax,P.FineSmoother.Nstep), + ThreeLevelPreconF(PVdagMF,SmootherGCRF,C.MrhsProjectorL1,L2PGCR,C.Grids.Coarse5d,C.CMrhs), + PreconSeam(ThreeLevelPreconF,Ddwf.FermionGrid(),DdwfF.FermionGrid(),nr), + L1PGCR(P.Outer.Tol,P.Outer.MaxIterations,PVdagM, + (P.Setup.FinePrecision==MGPrecision::fp32) + ? static_cast&>(PreconSeam) + : static_cast&>(ThreeLevelPrecon), + P.Outer.Mmax,P.Outer.Nstep) { GRID_ASSERT( C.DenseCC != nullptr ); - + CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1); @@ -480,23 +610,39 @@ public: L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1); - + SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1); + SmootherGCRF.Level(1); + SmootherGCRF.Name("Fsmoother"); + SmootherGCRF.SetZeroGuess(1); + L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1); - + // The V-cycle is preconditioner: sloppy halos inside, exact restored on exit. ThreeLevelPrecon.SetSloppy = [this](int s){ PVdagM.SloppyComms(s); }; ThreeLevelPrecon.SloppyComms = Params.Setup.FineSloppyComms; + ThreeLevelPreconF.SetSloppy = [this](int s){ PVdagMF.SloppyComms(s); }; + ThreeLevelPreconF.SloppyComms = Params.Setup.FineSloppyComms; PVdagM.SloppyComms(0); - std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs " << nr - << ", fine halo policy: preconditioner+coarsening " - << (Params.Setup.FineSloppyComms ? "SLOPPY (fp32 wire)" : "exact") - << ", outer Krylov EXACT" << std::endl; + PVdagMF.SloppyComms(0); + // Three stages, three arithmetics, and a wire format that follows the + // ARITHMETIC of the stage it serves: fp32 on an fp64 operator, bf16 on + // an fp32 one. The coarsening always applies the fp64 operator, so its + // wire is fp32 whatever the preconditioner runs in. + const bool sloppy = Params.Setup.FineSloppyComms; + const bool pfp32 = (Params.Setup.FinePrecision==MGPrecision::fp32); + std::cout << GridLogMessage << "PVdagMMultiGridSolver: Nrhs " << nr << std::endl; + std::cout << GridLogMessage << " outer Krylov : fp64, halos EXACT (the true residual is measured here)" << std::endl; + std::cout << GridLogMessage << " preconditioner fine: " << (pfp32 ? "fp32 (seam at the outer Krylov)" : "fp64") + << ", halos " << (sloppy ? (pfp32 ? "SLOPPY bf16 wire" : "SLOPPY fp32 wire") : "exact") << std::endl; + std::cout << GridLogMessage << " coarsening (done) : fp64 operator, halos " + << (sloppy ? "SLOPPY fp32 wire" : "exact") + << " -- the coarse operator does not follow FinePrecision" << std::endl; } void Solve(std::vector &src, std::vector &sol) diff --git a/Grid/algorithms/multigrid/PVdagMMultiGridParams.h b/Grid/algorithms/multigrid/PVdagMMultiGridParams.h index adfabd925..39e44d0c5 100644 --- a/Grid/algorithms/multigrid/PVdagMMultiGridParams.h +++ b/Grid/algorithms/multigrid/PVdagMMultiGridParams.h @@ -42,6 +42,12 @@ NAMESPACE_BEGIN(Grid); // not consumer knobs. ////////////////////////////////////////////////////////////////////////////////////// +// Arithmetic precision of the fine level INSIDE the preconditioner -- the +// smoother and the V-cycle's own fine residuals. The outer Krylov is +// always fp64: the fp64/fp32 seam sits at its preconditioner call. +// Orthogonal to FineSloppyComms, which is the halo WIRE format. +GRID_SERIALIZABLE_ENUM(MGPrecision, undef, fp64, 1, fp32, 2); + // The smoother is the adaptive shifted PGCR -- the one correct route for // this non-Hermitian chain (stationary replay and Chebyshev were explored // and did not win here; Chebyshev remains the HERMITIAN chain's smoother). @@ -79,14 +85,17 @@ struct MGDenseParams : Serializable { struct MGSetupParams : Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(MGSetupParams, - std::vector, Block, // fine -> coarse blocking + std::vector, Block1, // fine -> coarse blocking std::vector, Block2, // coarse -> coarse-coarse blocking int, CoarsenBatch, std::string, SubspaceFile, // scidac; empty = create from noise, no I/O - int, FineSloppyComms);// fp32 wire INSIDE the preconditioner only + int, FineSloppyComms, // reduced-precision halo WIRE inside the preconditioner only: fp32 on the fp64 operator, bf16 on the fp32 operator + MGPrecision, FinePrecision, // fine-level ARITHMETIC inside the preconditioner + int, RetainSubspace); // keep the RAW basis after setup (nbasis fine vectors); needed ONLY for a fixed-basis rebuild on a changed gauge field (HMC). A valence solve never rebuilds, so 0 frees it. MGSetupParams() - : Block({2,2,3,3}), Block2({4,4,2,4}), CoarsenBatch(9), - SubspaceFile(""), FineSloppyComms(1) {}; + : Block1({2,2,3,3}), Block2({4,4,2,4}), CoarsenBatch(9), + SubspaceFile(""), FineSloppyComms(1), FinePrecision(MGPrecision::fp64), + RetainSubspace(0) {}; }; struct PVdagMMultiGridParams : Serializable { @@ -104,9 +113,10 @@ struct PVdagMMultiGridParams : Serializable { inline void CheckValidity(const PVdagMMultiGridParams &P) { - GRID_ASSERT( P.Setup.Block.size() == 4 ); + GRID_ASSERT( P.Setup.Block1.size() == 4 ); GRID_ASSERT( P.Setup.Block2.size() == 4 ); GRID_ASSERT( P.Setup.CoarsenBatch >= 1 ); + GRID_ASSERT( P.Setup.FinePrecision != MGPrecision::undef ); GRID_ASSERT( P.FineSmoother.Nstep > 0 ); GRID_ASSERT( P.CoarseSmoother.Nstep > 0 ); GRID_ASSERT( P.CoarseSolver.Tol > 0.0 ); diff --git a/Grid/algorithms/multigrid/Smoothers.h b/Grid/algorithms/multigrid/Smoothers.h index 1fce1cdb7..6883905a8 100644 --- a/Grid/algorithms/multigrid/Smoothers.h +++ b/Grid/algorithms/multigrid/Smoothers.h @@ -36,6 +36,8 @@ NAMESPACE_BEGIN(Grid); // smoother operator such as shifted PVdagM. // ChebyshevInverter one Chebyshev-corrected step with residual print. // MirsSmoother shifted-MdagM CG, HDCG arXiv:1402.2585. +// CGSmoother fixed-iteration CG on an already-shifted +// Hermitian operator: the mrhs HDCG smoother. // GCRReplaySmoother replays a GCR's recorded step lengths a_k and // orthogonalisation coefficients b_kj with NO // inner products: one matvec per step, zero @@ -191,6 +193,30 @@ public: } }; +////////////////////////////////////////////////////////////////////////////// +// Fixed-iteration CG on an (already shifted) Hermitian operator, tolerance +// zero so the work is fixed. NOT a stationary preconditioner: the CG +// coefficients depend on the input, so the map is nonlinear. Use under a +// flexible outer (fPcg); BlockCGrQ needs the Chebyshev smoother. +////////////////////////////////////////////////////////////////////////////// +template class CGSmoother : public LinearFunction +{ +public: + using LinearFunction::operator(); + LinearOperatorBase &_Op; + int iters; + CGSmoother(int _iters, LinearOperatorBase &Op) : _Op(Op), iters(_iters) + { + std::cout << GridLogMessage << "CGSmoother order " << iters << std::endl; + } + void operator() (const Field &in, Field &out) + { + ConjugateGradient CG(0.0,iters,false); // never converges: by design + out = Zero(); + CG(_Op,in,out); + } +}; + ////////////////////////////////////////////////////////////////////////////// // Replay of a recorded GCR with a trivial preconditioner: // p_0 = r_0 ; x_{k+1} = x_k + a_k p_k ; r_{k+1} = r_k - a_k A p_k ; diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index b58037b19..3de398e88 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -210,6 +210,14 @@ accelerator_inline void convertType(vComplexD2 & out, const vComplexF & in) { precisionChange(out,in); } +// Precision change in the lex (unvectorised) chart: one lane, so a plain +// scalar conversion. Reached by the block operations on an fp32 coarse +// space, whose inner products accumulate in double and convert back. +accelerator_inline void convertType(sComplexF & out, const sComplexD & in) { out.v = ComplexF(in.v); } +accelerator_inline void convertType(sComplexD & out, const sComplexF & in) { out.v = ComplexD(in.v); } +accelerator_inline void convertType(sRealF & out, const sRealD & in) { out.v = RealF(in.v); } +accelerator_inline void convertType(sRealD & out, const sRealF & in) { out.v = RealD(in.v); } + template accelerator_inline void convertType(iScalar & out, const iScalar & in) { convertType(out._internal,in._internal); diff --git a/Grid/tensors/Tensor_traits.h b/Grid/tensors/Tensor_traits.h index 8ac8f7e00..f4f8048df 100644 --- a/Grid/tensors/Tensor_traits.h +++ b/Grid/tensors/Tensor_traits.h @@ -110,6 +110,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef RealD DoublePrecision; typedef RealD DoublePrecision2; + typedef RealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; @@ -124,6 +125,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef RealD DoublePrecision; typedef RealD DoublePrecision2; + typedef RealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; @@ -138,6 +140,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef ComplexD DoublePrecision; typedef ComplexD DoublePrecision2; + typedef ComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; @@ -152,6 +155,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef ComplexD DoublePrecision; typedef ComplexD DoublePrecision2; + typedef ComplexF SinglePrecision; }; #if defined(GRID_CUDA) || defined(GRID_HIP) @@ -168,6 +172,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef scalar_typeD DoublePrecision; typedef scalar_typeD DoublePrecision2; + typedef std::complex SinglePrecision; }; template<> struct GridTypeMapper > : public GridTypeMapper_Base { typedef std::complex scalar_type; @@ -182,6 +187,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef scalar_typeD DoublePrecision; typedef scalar_typeD DoublePrecision2; + typedef std::complex SinglePrecision; }; #endif @@ -198,6 +204,7 @@ NAMESPACE_BEGIN(Grid); typedef Integer Integerified; typedef void DoublePrecision; typedef void DoublePrecision2; + typedef void SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { @@ -213,6 +220,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vRealD DoublePrecision; typedef vRealD2 DoublePrecision2; + typedef vRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; @@ -227,6 +235,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vRealD DoublePrecision; typedef vRealD DoublePrecision2; + typedef vRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; @@ -241,6 +250,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vRealD2 DoublePrecision; typedef vRealD2 DoublePrecision2; + typedef vRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { // Fixme this is incomplete until Grid supports fp16 or bfp16 arithmetic types @@ -256,6 +266,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vRealD DoublePrecision; typedef vRealD DoublePrecision2; + typedef vRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { // Fixme this is incomplete until Grid supports fp16 or bfp16 arithmetic types @@ -271,6 +282,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vComplexD DoublePrecision; typedef vComplexD DoublePrecision2; + typedef vComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; @@ -285,6 +297,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vComplexD DoublePrecision; typedef vComplexD2 DoublePrecision2; + typedef vComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; @@ -299,6 +312,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vComplexD DoublePrecision; typedef vComplexD DoublePrecision2; + typedef vComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; @@ -313,6 +327,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef vComplexD2 DoublePrecision; typedef vComplexD2 DoublePrecision2; + typedef vComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; @@ -327,6 +342,7 @@ NAMESPACE_BEGIN(Grid); typedef vInteger Integerified; typedef void DoublePrecision; typedef void DoublePrecision2; + typedef void SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealF scalar_type; @@ -341,6 +357,7 @@ NAMESPACE_BEGIN(Grid); typedef sInteger Integerified; typedef sRealD DoublePrecision; typedef sRealD DoublePrecision2; + typedef sRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; @@ -355,6 +372,7 @@ NAMESPACE_BEGIN(Grid); typedef sInteger Integerified; typedef sRealD DoublePrecision; typedef sRealD DoublePrecision2; + typedef sRealF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; @@ -369,6 +387,7 @@ NAMESPACE_BEGIN(Grid); typedef sInteger Integerified; typedef sComplexD DoublePrecision; typedef sComplexD DoublePrecision2; + typedef sComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; @@ -383,6 +402,7 @@ NAMESPACE_BEGIN(Grid); typedef sInteger Integerified; typedef sComplexD DoublePrecision; typedef sComplexD DoublePrecision2; + typedef sComplexF SinglePrecision; }; template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; @@ -397,6 +417,7 @@ NAMESPACE_BEGIN(Grid); typedef sInteger Integerified; typedef void DoublePrecision; typedef void DoublePrecision2; + typedef void SinglePrecision; }; #define GridTypeMapper_RepeatedTypes \ @@ -417,6 +438,7 @@ NAMESPACE_BEGIN(Grid); using Integerified = iScalar; using DoublePrecision = iScalar; using DoublePrecision2= iScalar; + using SinglePrecision = iScalar; static constexpr int Rank = BaseTraits::Rank + 1; static constexpr std::size_t count = BaseTraits::count; static constexpr int Dimension(int dim) { @@ -433,6 +455,7 @@ NAMESPACE_BEGIN(Grid); using Integerified = iVector; using DoublePrecision = iVector; using DoublePrecision2= iVector; + using SinglePrecision = iVector; static constexpr int Rank = BaseTraits::Rank + 1; static constexpr std::size_t count = BaseTraits::count * N; static constexpr int Dimension(int dim) { @@ -449,6 +472,7 @@ NAMESPACE_BEGIN(Grid); using Integerified = iMatrix; using DoublePrecision = iMatrix; using DoublePrecision2= iMatrix; + using SinglePrecision = iMatrix; static constexpr int Rank = BaseTraits::Rank + 2; static constexpr std::size_t count = BaseTraits::count * N * N; static constexpr int Dimension(int dim) { diff --git a/benchmarks/Benchmark_usqcd.cc b/benchmarks/Benchmark_usqcd.cc index 34a0d6ecd..c739ac108 100644 --- a/benchmarks/Benchmark_usqcd.cc +++ b/benchmarks/Benchmark_usqcd.cc @@ -276,14 +276,31 @@ public: GridBLAS blas; + // The batched GEMMs of a coarse multigrid level are memory bound, not + // flop bound: each call streams the matrices once and does O(1) flops per + // element. So report the effective bandwidth beside the flop rate -- + // measured against the device's peak it says whether a shape is running + // well, which the flop rate alone does not. beta=1, so C is read as well + // as written. + auto Report = [&](int M,int N,int K,int BATCH) + { + double gf = blas.benchmark(M,N,K,BATCH); + double flops = 8.0*M*N*K*BATCH; + double bytes = 1.0*sizeof(CComplex)*(1.0*M*K + 1.0*K*N + 2.0*M*N)*BATCH; + double gbs = gf*bytes/flops; + fprintf(FP,"%d, %d, %d, %d, %f, %f\n", M, N, K, BATCH, gf, gbs); + std::cout<(M,N,K,BATCH); - - fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p); - - std::cout<(M,N,K,BATCH); - - fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p); - std::cout<(M,N,K,BATCH); - - fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p); - std::cout< using namespace std; using namespace Grid; -// Routes Op/AdjOp -> HermOp so that CoarsenOperator and CreateSubspace -// both see the HPD operator M†M rather than bare M. -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase &wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out) { wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { - wrapped.HermOp(in, out); - ComplexD dot = innerProduct(in, out); - n1 = real(dot); - n2 = norm2(out); - } -}; -// Fixed-iteration CG smoother: runs exactly `iters` steps of CG on the -// shifted operator. tolerance=0 so CG never exits early. -template -class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator &_SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) - : _SmootherOperator(SmootherOperator), iters(_iters) - { - std::cout << GridLogMessage << " CGSmoother order " << iters << std::endl; - } - void operator()(const Field &in, Field &out) - { - ConjugateGradient CG(0.0, iters, false); - out = Zero(); - CG(_SmootherOperator, in, out); - } -}; int main (int argc, char ** argv) { diff --git a/examples/Example_mdagm_cg.cc b/examples/Example_mdagm_cg.cc index ef8001609..1a67a9ffb 100644 --- a/examples/Example_mdagm_cg.cc +++ b/examples/Example_mdagm_cg.cc @@ -57,46 +57,7 @@ Author: Peter Boyle using namespace std; using namespace Grid; -// Wraps any LinearOperatorBase so that Op = AdjOp = HermOp. -// Required when coarsening an HPD operator whose Op != HermOp -// (e.g. MdagMLinearOperator where Op=M, HermOp=M†M). -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase &wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {} - void Op (const Field &in, Field &out) { wrapped.HermOp(in, out); } - void HermOp (const Field &in, Field &out) { wrapped.HermOp(in, out); } - void AdjOp (const Field &in, Field &out) { wrapped.HermOp(in, 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 HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { - wrapped.HermOp(in, out); - ComplexD dot = innerProduct(in, out); - n1 = real(dot); n2 = norm2(out); - } -}; -// Fixed-iteration CG as a smoother (LinearFunction). -// Used as the IR-shifted smoother: solves (M†M + lo*I) x = b approximately. -template -class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - LinearOperatorBase &_op; - int iters; - CGSmoother(int _iters, LinearOperatorBase &op) : _op(op), iters(_iters) { - std::cout << GridLogMessage << "CGSmoother order " << iters << std::endl; - } - void operator()(const Field &in, Field &out) { - ConjugateGradient CG(0.0, iters, false); - out = Zero(); - CG(_op, in, out); - } -}; // Two-level V-cycle preconditioner (LinearFunction). template diff --git a/examples/Example_pvdagm_mrhs.cc b/examples/Example_pvdagm_mrhs.cc index cf22afa9c..c0d14fbcb 100644 --- a/examples/Example_pvdagm_mrhs.cc +++ b/examples/Example_pvdagm_mrhs.cc @@ -58,6 +58,7 @@ Author: Peter Boyle #include #include #include +#include using namespace std; using namespace Grid; @@ -183,194 +184,6 @@ public: } }; -////////////////////////////////////////////////////////////////////// -// Minimal multi-RHS function interface (preconditioner slot) -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; - -////////////////////////////////////////////////////////////////////// -// Single-polynomial multi-RHS PGCR (non-Hermitian). -// -// Verbatim adaptation of PrecGeneralisedConjugateResidualNonHermitian -// to std::vector: every innerProduct / norm2 is SUMMED over the -// RHS index, so one alpha/beta per step is shared by all RHS -- the -// single GCR on the enlarged block-diagonal system. -////////////////////////////////////////////////////////////////////// -template -class MrhsPGCRNonHermitian { -public: - RealD Tolerance; - Integer MaxIterations; - int mmax; - int nstep; - int steps; - int level; - int ZeroGuess = 0; int FirstCycle = 0; // caller contract: zero guess => first-cycle r0 = src - std::string name = "Level 1"; - LinearOperatorBase &Linop; - MrhsLinearFunction &Preconditioner; - - void Level(int lv) { name = "Level " + std::to_string(lv); level=lv; }; - void Name(std::string n) { name = n; }; - 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; } - - /////////////////////////////////////////////////////////////// - // vector-of-fields linear algebra, reductions summed over rhs - /////////////////////////////////////////////////////////////// - static RealD vnorm2(std::vector &x){ - RealD s=0.0; for(auto &f : x) s+=norm2(f); return s; - } - static ComplexD vinnerProduct(std::vector &x, std::vector &y){ - ComplexD s(0.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){ - 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; - - std::vector r(nrhs,grid); - - GridStopWatch SolverTimer; - SolverTimer.Start(); - - steps=0; - FirstCycle=1; - for(int k=0;k &src, std::vector &psi, RealD rsq){ - - RealD cp; - ComplexD a, b, rq; - RealD zAAz; - - int nrhs = src.size(); - GridBase *grid = src[0].Grid(); - - std::vector r (nrhs,grid); - std::vector z (nrhs,grid); - std::vector Az(nrhs,grid); - - //////////////////////////////// - // history for flexible orthog: [mmax][nrhs] - //////////////////////////////// - std::vector< std::vector > q(mmax, std::vector(nrhs,grid)); - std::vector< std::vector > p(mmax, std::vector(nrhs,grid)); - std::vector qq(mmax); - - std::cout<(mmax-1))?(mmax-1):(kp); - for(int back=0;back=0); - b = -real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; - vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); - vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); - } - qq[peri_kp]=vnorm2(q[peri_kp]); - } - GRID_ASSERT(0); // never reached - return cp; - } -}; - ////////////////////////////////////////////////////////////////////// // Trivial multi-RHS preconditioner ////////////////////////////////////////////////////////////////////// @@ -382,98 +195,6 @@ public: } }; -////////////////////////////////////////////////////////////////////// -// MultiRHS two-level V-cycle. -// -// Mirrors MGPreconditioner in Example_pvdagm: -// out = in (trivial pre) [per rhs] -// r1 = in - A out [per rhs] -// batched blockProject -> pack -> ONE mrhs coarse PGCR -> unpack -// -> batched blockPromote; out += correction -// r2 = in - A out [per rhs] -// per-RHS fine post-smoother; out += smooth(r2) -////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { -public: - typedef MrhsCoarseVector CoarseVector; // same lattice type on Coarse5d and CoarseMrhs - - LinearOperatorBase &_FineOperator; - FineSmoother &_PostSmoother; // single-RHS smoother, looped - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseSolve; // PGCR on the 6D mrhs field - GridBase *_CoarseGrid; // Coarse5d (single rhs) - GridBase *_CoarseGridMrhs; // 6D - - MrhsTwoLevelMG(LinearOperatorBase &FineOp, - FineSmoother &Post, - MultiRHSBlockProject &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){ - int nrhs = in.size(); - GridBase *fgrid = in[0].Grid(); - double t; - - std::vector vec1(nrhs,fgrid); - std::vector vec2(nrhs,fgrid); - - // Trivial pre-smoother: out = in (as in Example_pvdagm with simple_fine) - for(int r=0;rcoarse, pack rhs into 6D field - std::vector Csrc_split(nrhs,_CoarseGrid); - std::vector Csol_split(nrhs,_CoarseGrid); - CoarseVector CsrcMrhs(_CoarseGridMrhs); - CoarseVector CsolMrhs(_CoarseGridMrhs); - - t=-usecond(); - _Projector.blockProject(vec1,Csrc_split); - for(int r=0;rfine, add correction - t=-usecond(); - for(int r=0;r #include #include #include +#include using namespace std; using namespace Grid; @@ -167,206 +168,6 @@ public: void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } }; -////////////////////////////////////////////////////////////////////// -// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from Example_pvdagm_mrhs.cc): -// reductions summed over rhs -> one alpha/beta per step for the enlarged system. -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; - -template -class MrhsPGCRNonHermitian { -public: - RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; - int ZeroGuess = 0; int FirstCycle = 0; // caller contract: zero guess => first-cycle r0 = src - std::string name = "Level 1"; - LinearOperatorBase &Linop; - MrhsLinearFunction &Preconditioner; - void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } - void Name(std::string n){ name = n; } - 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){ 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; - std::vector r(nrhs,grid); - GridStopWatch T; T.Start(); steps=0; FirstCycle=1; - for(int k=0;k &src,std::vector &psi,RealD rsq){ - RealD cp; ComplexD a,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid(); - std::vector r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid); - std::vector< std::vector > q(mmax,std::vector(nrhs,grid)); - std::vector< std::vector > p(mmax,std::vector(nrhs,grid)); - std::vector qq(mmax); - std::cout<(mmax-1))?(mmax-1):(kp); - for(int back=0;back=0); - b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; - vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } - qq[peri_kp]=vnorm2(q[peri_kp]); - } - GRID_ASSERT(0); return cp; - } -}; - -////////////////////////////////////////////////////////////////////// -// L2->L3 mrhs V-cycle: a LinearFunction on the 6D mrhs COARSE field. -// Mirrors Example_pvdagm_mrhs.cc's MrhsTwoLevelMG one level down, and the -// single-RHS MGPreconditioner of Example_pvdagm_3level_SVDdefl.cc: -// out = in (trivial pre) -// r = in - A_coarse out -// restrict (unpack 6D coarse -> blockProject -> pack 6D coarse-coarse) -// ONE coarse-coarse solve (L3, GEMM) -// prolong (unpack -> blockPromote -> pack); out += correction -// r = in - A_coarse out -// coarse smoother (shifted 6D coarse op); out += smooth(r) -////////////////////////////////////////////////////////////////////// -template -class MrhsCoarseThreeLevelPrec : public LinearFunction { -public: - LinearOperatorBase &_CoarseOp; // mrhs coarse op (6D) - LinearFunction &_CoarseSmoother; // shifted 6D coarse smoother - MultiRHSBlockProject &_Projector; // L2->L3 (vector-based) - LinearFunction &_CoarseCoarseSolve; // L3 solve (6D cc) - GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs; - int _nrhs; - - MrhsCoarseThreeLevelPrec(LinearOperatorBase &CoarseOp, - LinearFunction &CoarseSmoother, - MultiRHSBlockProject &Projector, - LinearFunction &CoarseCoarseSolve, - GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs) - : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector), - _CoarseCoarseSolve(CoarseCoarseSolve), - _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {} - - using LinearFunction::operator(); - virtual void operator()(const CoarseField &in, CoarseField &out) { - int nrhs=_nrhs; double t; - CoarseField vec1(in.Grid()); - CoarseField vec2(in.Grid()); - - // trivial pre-smoother - out = in; - - // residual (6D coarse) - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - - // restrict: unpack 6D coarse -> vector -> blockProject -> vector -> pack 6D cc - std::vector csplit(nrhs,_Coarse5d); - std::vector ccsplit(nrhs,_CoarseCoarse5d); - CoarseCoarseField CCsrc(_CoarseCoarseMrhs); - CoarseCoarseField CCsol(_CoarseCoarseMrhs); - - t=-usecond(); - for(int r=0;rL3 restrict took "< blockPromote -> pack 6D coarse; add correction - t=-usecond(); - for(int r=0;rL3 prolong took "<L2 mrhs V-cycle (verbatim from Example_pvdagm_mrhs.cc): -// per-rhs fine smoother + batched restriction + ONE coarse solve + batched prolong. -// The coarse solve passed in is now itself three-level (preconditioned by L2->L3). -////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { -public: - typedef MrhsCoarseVector CoarseVector; - LinearOperatorBase &_FineOperator; - FineSmoother &_PostSmoother; - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseSolve; - GridBase *_CoarseGrid, *_CoarseGridMrhs; - MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, - MultiRHSBlockProject &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){ - int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t; - std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); - for(int r=0;r Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); - CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); - t=-usecond(); - _Projector.blockProject(vec1,Csrc_split); - for(int r=0;r #include #include #include +#include #include #include @@ -173,281 +174,6 @@ public: void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); } }; -////////////////////////////////////////////////////////////////////// -// Dense CC solve on the PACKED 6D mrhs coarse-coarse field: drop-in -// for the L3 PGCR, delegating to the library DenseCoarseMatrix. -////////////////////////////////////////////////////////////////////// -template -class MrhsDenseCCSolve : public LinearFunction { -public: - DenseType &_Dense; - GridBase *_CoarseCoarse5d; - int _nrhs; - MrhsDenseCCSolve(DenseType &D, GridBase *cc5d, int nrhs) - : _Dense(D), _CoarseCoarse5d(cc5d), _nrhs(nrhs) {} - using LinearFunction::operator(); - virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){ - if ( getenv("DENSE_CC_CHECK") ) { - // Audit path: per-rhs 5D unpack so ApplyBatch can run the _Op defect - // check per rhs. ~50ms/call of slice/split overhead -- audit only. - CoarseCoarseField tmp(in.Grid()); - tmp = in; - std::vector split_in (_nrhs,_CoarseCoarse5d); - std::vector split_out(_nrhs,_CoarseCoarse5d); - for(int r=0;r<_nrhs;r++) ExtractSliceFast(split_in[r], tmp, r, 0); - _Dense.ApplyBatch(split_in, split_out); - for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0); - } else { - GRID_TRACE("MrhsDensCCSolve::ApplyBatch6D"); - _Dense.ApplyBatch6D(in, out, _nrhs); - } - } -}; - -////////////////////////////////////////////////////////////////////// -// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from the -// frozen Example_pvdagm_mrhs_3level_dense.cc) -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; - -template -class MrhsPGCRNonHermitian { -public: - RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; - int ZeroGuess = 0; int FirstCycle = 0; // caller contract: zero guess => first-cycle r0 = src - std::string name = "Level 1"; - LinearOperatorBase &Linop; - MrhsLinearFunction &Preconditioner; - void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } - void Name(std::string n){ name = n; } - 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){ GRID_TRACE("GCR-vnorm2"); RealD s=0; for(auto &f:x) s+=norm2(f); return s; } - static ComplexD vinnerProduct(std::vector &x,std::vector &y){ GRID_TRACE("GCR-vinnerProduct"); 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){ GRID_TRACE("GCR-vaxpy"); 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("GCR-vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); } - void operator()(std::vector &src,std::vector &psi){ - GRID_TRACE((name+"MrhsPGCRNonHermitian").c_str()); - 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; - for(int k=0;k &src,std::vector &psi,RealD rsq){ - RealD cp; ComplexD a,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid(); - std::vector r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid); - std::vector< std::vector > q(mmax,std::vector(nrhs,grid)); - std::vector< std::vector > p(mmax,std::vector(nrhs,grid)); - std::vector qq(mmax); - std::cout<(mmax-1))?(mmax-1):(kp); - for(int back=0;back=0); - b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; - vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } - qq[peri_kp]=vnorm2(q[peri_kp]); - } - // q[peri_kp]=Az; p[peri_kp]=z; - // int northog=((kp)>(mmax-1))?(mmax-1):(kp); - // for(int back=0;back=0); - // b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; - // vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } - // qq[peri_kp]=vnorm2(q[peri_kp]); - } - GRID_ASSERT(0); return cp; - } -}; - -////////////////////////////////////////////////////////////////////// -// L2->L3 mrhs V-cycle: LinearFunction on the 6D mrhs COARSE field. -// The coarse-coarse solve slot takes EITHER the dense mrhs solve -// (DENSE_CC=1) or the L3 PGCR (DENSE_CC=0). -////////////////////////////////////////////////////////////////////// -template -class MrhsCoarseThreeLevelPrec : public LinearFunction { -public: - LinearOperatorBase &_CoarseOp; // mrhs coarse op (6D) - LinearFunction &_CoarseSmoother; // shifted 6D coarse smoother - MultiRHSBlockProject &_Projector; // L2->L3 (vector-based) - LinearFunction &_CoarseCoarseSolve; // L3 solve (6D cc) - GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs; - int _nrhs; - - MrhsCoarseThreeLevelPrec(LinearOperatorBase &CoarseOp, - LinearFunction &CoarseSmoother, - MultiRHSBlockProject &Projector, - LinearFunction &CoarseCoarseSolve, - GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs) - : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector), - _CoarseCoarseSolve(CoarseCoarseSolve), - _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {} - - using LinearFunction::operator(); - virtual void operator()(const CoarseField &in, CoarseField &out) { - int nrhs=_nrhs; double t; - CoarseField vec1(in.Grid()); - CoarseField vec2(in.Grid()); - - // trivial pre-smoother - out = in; - - // residual (6D coarse) - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - - // restrict: unpack 6D coarse -> vector -> blockProject -> vector -> pack 6D cc - std::vector csplit(nrhs,_Coarse5d); - std::vector ccsplit(nrhs,_CoarseCoarse5d); - CoarseCoarseField CCsrc(_CoarseCoarseMrhs); - CoarseCoarseField CCsol(_CoarseCoarseMrhs); - - t=-usecond(); - { - GRID_TRACE("L2L3-Vcycle - Extract/blockProject/Insert"); - for(int r=0;rL3 restrict took "< blockPromote -> pack 6D coarse; add correction - t=-usecond(); - { - GRID_TRACE("L2L3-Vcycle - Ext/blockPromote/Ins"); - for(int r=0;rL3 prolong took "<L2 mrhs V-cycle (verbatim from the frozen example) -////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { -public: - typedef MrhsCoarseVector CoarseVector; - LinearOperatorBase &_FineOperator; - FineSmoother &_PostSmoother; - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseSolve; - GridBase *_CoarseGrid, *_CoarseGridMrhs; - MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, - MultiRHSBlockProject &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){ - int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t; - - std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); - { - GRID_TRACE("L1L2-Vcycle - Resid"); - for(int r=0;r Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); - CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); - t=-usecond(); - { - GRID_TRACE("L1L2-Vcycle - blockProject"); - _Projector.blockProject(vec1,Csrc_split); - for(int r=0;rImport(LittleDiracOpL2); - MrhsDenseCC.reset(new MrhsDenseCCSolve(*DenseCC, CoarseCoarse5d, nrhs)); + MrhsDenseCC.reset(new MrhsDenseCCSolve(*DenseCC, nrhs)); } ////////////////////////////////////////////////////////////////////// diff --git a/examples/Example_pvdagm_mrhs_3level_dense.cc b/examples/Example_pvdagm_mrhs_3level_dense.cc index b70416852..320e63644 100644 --- a/examples/Example_pvdagm_mrhs_3level_dense.cc +++ b/examples/Example_pvdagm_mrhs_3level_dense.cc @@ -60,6 +60,7 @@ Author: Peter Boyle #include #include #include +#include #include #include @@ -830,226 +831,6 @@ public: } }; -////////////////////////////////////////////////////////////////////// -// Dense CC solve on the PACKED 6D mrhs coarse-coarse field: unpack per -// rhs, ONE batched dense apply, repack. Drop-in for the L3 PGCR. -////////////////////////////////////////////////////////////////////// -template -class MrhsDenseCCSolve : public LinearFunction { -public: - DistributedDenseInverse &_Dense; - GridBase *_CoarseCoarse5d; - int _nrhs; - MrhsDenseCCSolve(DistributedDenseInverse &D, GridBase *cc5d, int nrhs) - : _Dense(D), _CoarseCoarse5d(cc5d), _nrhs(nrhs) {} - using LinearFunction::operator(); - virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){ - if ( getenv("DENSE_CC_CHECK") ) { - // Audit path: per-rhs 5D unpack so ApplyBatch can run the _Op defect - // check per rhs. ~50ms/call of slice/split overhead -- audit only. - CoarseCoarseField tmp(in.Grid()); - tmp = in; - std::vector split_in (_nrhs,_CoarseCoarse5d); - std::vector split_out(_nrhs,_CoarseCoarse5d); - for(int r=0;r<_nrhs;r++) ExtractSliceFast(split_in[r], tmp, r, 0); - _Dense.ApplyBatch(split_in, split_out); - for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0); - } else { - _Dense.ApplyBatch6D(in, out, _nrhs); - } - } -}; - -////////////////////////////////////////////////////////////////////// -// mrhs interfaces + single-polynomial mrhs PGCR (verbatim from Example_pvdagm_mrhs.cc) -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; - -template -class MrhsPGCRNonHermitian { -public: - RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; - int ZeroGuess = 0; int FirstCycle = 0; // caller contract: zero guess => first-cycle r0 = src - std::string name = "Level 1"; - LinearOperatorBase &Linop; - MrhsLinearFunction &Preconditioner; - void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } - void Name(std::string n){ name = n; } - 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){ 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; - std::vector r(nrhs,grid); - GridStopWatch T; T.Start(); steps=0; FirstCycle=1; - for(int k=0;k &src,std::vector &psi,RealD rsq){ - RealD cp; ComplexD a,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid(); - std::vector r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid); - std::vector< std::vector > q(mmax,std::vector(nrhs,grid)); - std::vector< std::vector > p(mmax,std::vector(nrhs,grid)); - std::vector qq(mmax); - std::cout<(mmax-1))?(mmax-1):(kp); - for(int back=0;back=0); - b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; - vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } - qq[peri_kp]=vnorm2(q[peri_kp]); - } - GRID_ASSERT(0); return cp; - } -}; - -////////////////////////////////////////////////////////////////////// -// L2->L3 mrhs V-cycle: LinearFunction on the 6D mrhs COARSE field. -// The coarse-coarse solve slot now takes EITHER the dense mrhs solve -// (DENSE_CC=1) or the L3 PGCR (DENSE_CC=0). -////////////////////////////////////////////////////////////////////// -template -class MrhsCoarseThreeLevelPrec : public LinearFunction { -public: - LinearOperatorBase &_CoarseOp; // mrhs coarse op (6D) - LinearFunction &_CoarseSmoother; // shifted 6D coarse smoother - MultiRHSBlockProject &_Projector; // L2->L3 (vector-based) - LinearFunction &_CoarseCoarseSolve; // L3 solve (6D cc) - GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs; - int _nrhs; - - MrhsCoarseThreeLevelPrec(LinearOperatorBase &CoarseOp, - LinearFunction &CoarseSmoother, - MultiRHSBlockProject &Projector, - LinearFunction &CoarseCoarseSolve, - GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs) - : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector), - _CoarseCoarseSolve(CoarseCoarseSolve), - _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {} - - using LinearFunction::operator(); - virtual void operator()(const CoarseField &in, CoarseField &out) { - int nrhs=_nrhs; double t; - CoarseField vec1(in.Grid()); - CoarseField vec2(in.Grid()); - - // trivial pre-smoother - out = in; - - // residual (6D coarse) - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - - // restrict: unpack 6D coarse -> vector -> blockProject -> vector -> pack 6D cc - std::vector csplit(nrhs,_Coarse5d); - std::vector ccsplit(nrhs,_CoarseCoarse5d); - CoarseCoarseField CCsrc(_CoarseCoarseMrhs); - CoarseCoarseField CCsol(_CoarseCoarseMrhs); - - t=-usecond(); - for(int r=0;rL3 restrict took "< blockPromote -> pack 6D coarse; add correction - t=-usecond(); - for(int r=0;rL3 prolong took "<L2 mrhs V-cycle (verbatim from Example_pvdagm_mrhs.cc) -////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { -public: - typedef MrhsCoarseVector CoarseVector; - LinearOperatorBase &_FineOperator; - FineSmoother &_PostSmoother; - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseSolve; - GridBase *_CoarseGrid, *_CoarseGridMrhs; - MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, - MultiRHSBlockProject &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){ - int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t; - std::vector vec1(nrhs,fgrid),vec2(nrhs,fgrid); - for(int r=0;r Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); - CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); - t=-usecond(); - _Projector.blockProject(vec1,Csrc_split); - for(int r=0;r> DenseCC; - std::unique_ptr> MrhsDenseCC; + std::unique_ptr,CoarseCoarseVector>> MrhsDenseCC; if (UseDenseCC) { std::cout << GridLogMessage << "**********************************************" << std::endl; std::cout << GridLogMessage << " Dense CC inverse setup (mrhs bottom)" << std::endl; std::cout << GridLogMessage << "**********************************************" << std::endl; DenseCC.reset(new DistributedDenseInverse(LinOpCC5d, CoarseCoarse5d, nbasis)); - MrhsDenseCC.reset(new MrhsDenseCCSolve(*DenseCC, CoarseCoarse5d, nrhs)); + MrhsDenseCC.reset(new MrhsDenseCCSolve,CoarseCoarseVector>(*DenseCC, nrhs)); } ////////////////////////////////////////////////////////////////////// diff --git a/examples/Example_pvdagm_multigrid.cc b/examples/Example_pvdagm_multigrid.cc index 06bf00cb8..a52b64c0a 100644 --- a/examples/Example_pvdagm_multigrid.cc +++ b/examples/Example_pvdagm_multigrid.cc @@ -33,7 +33,11 @@ Author: Peter Boyle // parameter. The effective parameters are always printed, so the log // describes its own run. // -// Compile-time: -DNBASIS=8 cuts the basis down for laptop runs. +// Compile-time: NBASIS defaults to 60, the production basis; a subspace file +// may hold more vectors, the load reads only the first NBASIS. -DNBASIS=8 +// cuts it down for laptop runs. +// examples/Makefile.am builds the fp32-coarse and fp32-dense-inversion +// variants below as their own binaries. // #include #include @@ -48,6 +52,16 @@ using namespace Grid; #define NBASIS 60 #endif +// Precision of the coarse + coarse-coarse sector is a compile-time +// instantiation: -DCOARSE_SINGLE builds the fp32 coarse space (the dense +// bottom's apply slab is fp32 either way; its inversion is a configure +// option). Both levels carry the same site type iVector. +#ifdef COARSE_SINGLE +typedef sTComplexF CoarseScalar_t; +#else +typedef sTComplexD CoarseScalar_t; +#endif + struct PVdagMDriverParams : Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(PVdagMDriverParams, int, Ls, @@ -123,22 +137,46 @@ int main (int argc, char ** argv) ////////////////////////////////////////////////////////////////////// // Grids -> coarsening -> dense bottom. Scope order is lifetime order. + // + // The fine operator applied during setup is always fp64, on the fp64 + // basis. Everything the coarsening produces -- the Galerkin matrix + // elements and the transfer operator's store -- takes the coarse precision + // (CoarseScalar_t, a compile-time choice). + // + // There is exactly ONE fine transfer operator. Its STORE follows the + // coarse sector, since the coarse space is what it feeds, while its import + // and export accept either fine precision: they are already a layout + // transformation, and a scalar conversion inside one is free. So the fp64 + // setup and an fp32 V-cycle share the same object and the same basis store. + // FinePrecision selects only the fine operator and smoother; the outer + // Krylov and its true-residual check stay fp64 throughout. ////////////////////////////////////////////////////////////////////// - typedef PVdagMMultiGridCoarsening Coarsening_t; + typedef PVdagMMultiGridCoarsening Coarsening_t; + std::cout << GridLogMessage << "Coarse sector precision (compiled): " + << (sizeof(typename GridTypeMapper::scalar_type)==sizeof(ComplexF) ? "fp32" : "fp64") + << ", nbasis " << NBASIS << std::endl; MGCoarseGrids CGrids(FGrid, P.MultiGrid.Setup); - Coarsening_t Coarsening(CGrids, P.MultiGrid.Setup); + MGFineGridsF FGridsF(FGrid); + + LatticeGaugeFieldF UmuF(FGridsF.UGridF); + precisionChange(UmuF,Umu); + MobiusFermionF DdwfF(UmuF,*FGridsF.FGridF,*FGridsF.FrbGridF,*FGridsF.UGridF,*FGridsF.UrbGridF,P.Mass,P.M5,P.MobiusB,P.MobiusC); + MobiusFermionF DpvF (UmuF,*FGridsF.FGridF,*FGridsF.FrbGridF,*FGridsF.UGridF,*FGridsF.UrbGridF,1.0, P.M5,P.MobiusB,P.MobiusC); + + Coarsening_t Coarsening(CGrids, FGridsF, P.MultiGrid.Setup); Coarsening.GetSubspace(RNG5, PVdagM); Coarsening.Coarsen(PVdagM); Coarsening.BuildDenseBottom(); + Coarsening.CertifyCoarsening(PVdagM); ////////////////////////////////////////////////////////////////////// // Solves: mrhs then (optionally) single RHS through the SAME objects. ////////////////////////////////////////////////////////////////////// auto RunSolve = [&](int nr) { - PVdagMMultiGridSolver Solver(Ddwf,Dpv,Coarsening,P.MultiGrid,nr); + PVdagMMultiGridSolver Solver(Ddwf,Dpv,DdwfF,DpvF,Coarsening,P.MultiGrid,nr); std::vector src(nr,FGrid), sol(nr,FGrid); for(int r=0;r #include #include #include +#include #include #include @@ -396,220 +397,6 @@ void PowerIteration(const std::string &name, LinearOperatorBase &Op, Grid << std::endl; } -////////////////////////////////////////////////////////////////////// -// Dense L3 solve on the packed D+1 coarse-coarse field -////////////////////////////////////////////////////////////////////// -template -class MrhsDenseCCSolve : public LinearFunction { -public: - DenseType &_Dense; - int _nrhs; - MrhsDenseCCSolve(DenseType &D, int nrhs) : _Dense(D), _nrhs(nrhs) {} - using LinearFunction::operator(); - virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){ - _Dense.ApplyBatch6D(in, out, _nrhs); - } -}; - -////////////////////////////////////////////////////////////////////// -// mrhs interfaces + single-polynomial mrhs PGCR -////////////////////////////////////////////////////////////////////// -template -class MrhsLinearFunction { -public: - virtual void operator()(std::vector &in, std::vector &out) = 0; -}; - -template -class MrhsPGCRNonHermitian { -public: - RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level; - int ZeroGuess = 0; int FirstCycle = 0; - std::string name = "Level 1"; - LinearOperatorBase &Linop; - MrhsLinearFunction &Preconditioner; - std::function OnStep; // called with the outer step count after every step (smoother switching) - void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } - void Name(std::string n){ name = n; } - 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("MrhsPGCR::vOp"); 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; - std::vector r(nrhs,grid); - GridStopWatch T; T.Start(); steps=0; FirstCycle=1; - for(int k=0;k &src,std::vector &psi,RealD rsq){ - RealD cp; ComplexD a,b,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); - MemoryManager::Snapshot(name+" orthog begin step "+std::to_string(steps)); - { - GRID_TRACE("MrhsPGCR orthog"); - // Classical Gram-Schmidt: all coefficients against the UN-updated new q - // (independent, batchable), then apply. Complex coefficient: the - // operator is non-Hermitian, real() alone left q's non-orthogonal. - // Batched per rhs (one fused kernel + one reduction each), the shared - // coefficient summed over rhs on the host, ONE GlobalSumVector. - 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]; } - 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;backL3 mrhs V-cycle on the D+1 coarse field -////////////////////////////////////////////////////////////////////// -template -class MrhsCoarseThreeLevelPrec : public LinearFunction { -public: - LinearOperatorBase &_CoarseOp; - LinearFunction &_CoarseSmoother; - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseCoarseSolve; - GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs; - int _nrhs; - MrhsCoarseThreeLevelPrec(LinearOperatorBase &CoarseOp, - LinearFunction &CoarseSmoother, - MultiRHSBlockProject &Projector, - LinearFunction &CoarseCoarseSolve, - GridBase *Coarse5d, GridBase *CoarseCoarse5d, GridBase *CoarseCoarseMrhs, int nrhs) - : _CoarseOp(CoarseOp), _CoarseSmoother(CoarseSmoother), _Projector(Projector), - _CoarseCoarseSolve(CoarseCoarseSolve), - _Coarse5d(Coarse5d), _CoarseCoarse5d(CoarseCoarse5d), _CoarseCoarseMrhs(CoarseCoarseMrhs), _nrhs(nrhs) {} - using LinearFunction::operator(); - virtual void operator()(const CoarseField &in, CoarseField &out) { - int nrhs=_nrhs; - 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); - - // CCsol=Zero(); // Is this necessary? - _CoarseCoarseSolve(CCsrc,CCsol); - - _Projector.blockPromote(vec1,CCsol); - add(out,out,vec1); - - _CoarseOp.Op(out,vec1); sub(vec1,in,vec1); - // vec2=Zero(); // Zero guess - _CoarseSmoother(vec1,vec2); - add(out,out,vec2); - } -}; - -////////////////////////////////////////////////////////////////////// -// L1->L2 mrhs V-cycle -////////////////////////////////////////////////////////////////////// -template -class MrhsTwoLevelMG : public MrhsLinearFunction { -public: - typedef MrhsCoarseVector CoarseVector; - LinearOperatorBase &_FineOperator; - FineSmoother &_PostSmoother; - MultiRHSBlockProject &_Projector; - LinearFunction &_CoarseSolve; - GridBase *_CoarseGrid, *_CoarseGridMrhs; - MrhsTwoLevelMG(LinearOperatorBase &FineOp, FineSmoother &Post, - MultiRHSBlockProject &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){ - // The whole V-cycle is preconditioner: its fine residuals and the - // smoother run with sloppy halos; the caller (the outer Krylov) gets - // the exact operator back on exit. - GRID_TRACE("MGVcycle"); - SetFineSloppy(FineSloppyComms); - 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"); - _Projector.blockProject(vec1,CsrcMrhs); - } - CsolMrhs=Zero(); - { GRID_TRACE("MGCoarseSolve"); - _CoarseSolve(CsrcMrhs,CsolMrhs); - } - { GRID_TRACE("MGPromote"); - _Projector.blockPromote(vec1,CsolMrhs); - for(int r=0;r 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); diff --git a/examples/Makefile.am b/examples/Makefile.am index 31cbdc860..4a50ac030 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -2,5 +2,23 @@ SUBDIRS = . include Make.inc +# Precision variants of the two multigrid drivers. The coarse-sector and +# dense-inversion precisions are compile-time instantiations, so each is its +# own binary built from the same source; the fine-level precision is a +# run-time parameter and needs no variant. Per-target flags give each +# variant its own object file. +bin_PROGRAMS += Example_pvdagm_multigrid_fp32coarse \ + Example_pvdagm_multigrid_fp32coarse_fp32dense \ + Example_hdcg_multigrid_fp32coarse +Example_pvdagm_multigrid_fp32coarse_SOURCES = Example_pvdagm_multigrid.cc +Example_pvdagm_multigrid_fp32coarse_CPPFLAGS = -DCOARSE_SINGLE +Example_pvdagm_multigrid_fp32coarse_LDADD = $(top_builddir)/Grid/libGrid.a +Example_pvdagm_multigrid_fp32coarse_fp32dense_SOURCES = Example_pvdagm_multigrid.cc +Example_pvdagm_multigrid_fp32coarse_fp32dense_CPPFLAGS = -DCOARSE_SINGLE -DGRID_DENSE_INVERSE_SINGLE +Example_pvdagm_multigrid_fp32coarse_fp32dense_LDADD = $(top_builddir)/Grid/libGrid.a + +Example_hdcg_multigrid_fp32coarse_SOURCES = Example_hdcg_multigrid.cc +Example_hdcg_multigrid_fp32coarse_CPPFLAGS = -DCOARSE_SINGLE +Example_hdcg_multigrid_fp32coarse_LDADD = $(top_builddir)/Grid/libGrid.a diff --git a/systems/Frontier/pvdagm_multigrid.job b/systems/Frontier/pvdagm_multigrid.job index 92e88d1ae..017a7809c 100644 --- a/systems/Frontier/pvdagm_multigrid.job +++ b/systems/Frontier/pvdagm_multigrid.job @@ -21,10 +21,14 @@ # only code path (2D block-cyclic dense inverse, ring-allgather apply, # inverseLU big leaves) and are gone from the environment. # -# The binary must be built with -DNBASIS=64 to match the subspace file. +# The binary is a plain `make` target with the default basis, 60. The +# subspace load reads the first 60 vectors of the file. +# selects the arithmetic of the fine level +# INSIDE the preconditioner (fp64 here; the outer Krylov is fp64 and exact +# either way). The precision sweep is in pvdagm_mixed_precision.job. ############################################################################## -root=$HOME/PVdagM/Grid/systems/Frontier +root=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/MGrewrite/Grid/systems/Frontier source $root/sourceme-rocm7.2.sh # Paths substituted into the XML below. @@ -48,6 +52,8 @@ exec numactl -m \$NUMA -N \$NUMA \$* EOF chmod +x ./select_gpu +# sourceme sets FI_HMEM_ROCR_USE_DMABUF=0 -- the standing avoidance of the CXI +# NO_TRANSLATION fault at NRHS>=12 (libfabric #12775). export OMP_NUM_THREADS=7 ulimit -c 0 # no 22 GB GPU core dumps @@ -89,11 +95,13 @@ cat << EOF > params.xml 1 - 2233 + 2233 4424 9 $SUBSPACE 1 + fp64 + 0 0.164 diff --git a/tests/debug/Test_general_coarse.cc b/tests/debug/Test_general_coarse.cc index b47693568..2d0ca32cb 100644 --- a/tests/debug/Test_general_coarse.cc +++ b/tests/debug/Test_general_coarse.cc @@ -40,29 +40,6 @@ using namespace Grid; // a second definition here is a duplicate-symbol link error the moment the // archive member is pulled in. -/////////////////////// -// Tells little dirac op to use MdagM as the .Op() -/////////////////////// -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - 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){ - wrapped.HermOp(in,out); - } - void AdjOp (const Field &in, Field &out){ - wrapped.HermOp(in,out); - } - void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } - void HermOp(const Field &in, Field &out){ - wrapped.HermOp(in,out); - } -}; int main (int argc, char ** argv) diff --git a/tests/debug/Test_general_coarse_hdcg.cc b/tests/debug/Test_general_coarse_hdcg.cc index 9c6be62eb..73c3b2721 100644 --- a/tests/debug/Test_general_coarse_hdcg.cc +++ b/tests/debug/Test_general_coarse_hdcg.cc @@ -33,45 +33,6 @@ Author: Peter Boyle using namespace std; using namespace Grid; -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; - -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; - int main (int argc, char ** argv) { @@ -270,15 +231,9 @@ int main (int argc, char ** argv) ShiftedHermOpLinearOperator ShiftedFineHermOp(HermOpEO,MirsShift); CGSmoother CGsmooth(ord,ShiftedFineHermOp) ; - TwoLevelADEF2mrhs - HDCGmrhs(1.0e-8, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 500, FineHermOp, ADEF2, FrbGrid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid); diff --git a/tests/debug/Test_general_coarse_hdcg_phys.cc b/tests/debug/Test_general_coarse_hdcg_phys.cc index aff41901a..6e2145c50 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys.cc @@ -113,21 +113,6 @@ void LoadBasis(aggregation &Agg, std::string file) } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; /* void operator() (const Field &in, Field &out) { @@ -137,25 +122,6 @@ public: } }; */ -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - CG(_SmootherOperator,in,out); - } -}; int main (int argc, char ** argv) diff --git a/tests/debug/Test_general_coarse_hdcg_phys48.cc b/tests/debug/Test_general_coarse_hdcg_phys48.cc index 5383ff60b..910e404f6 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys48.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys48.cc @@ -155,44 +155,7 @@ void LoadEigenvectors(std::vector &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; int main (int argc, char ** argv) @@ -408,15 +371,9 @@ int main (int argc, char ** argv) MrhsGuesser.ImportEigenBasis(evec,eval); CGSmoother CGsmooth(Refineord,ShiftedFineHermOp) ; - TwoLevelADEF2mrhs - HDCGmrhsRefine(RefineTol, 500, - RefineFineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhsRefine_ADEF2(RefineFineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhsRefine(RefineTol, 500, RefineFineHermOp, HDCGmrhsRefine_ADEF2, MrhsProjector.fine_grid); // Reload the first pass aggregates, because we orthogonalised them LoadBasis(Aggregates,subspace_file); @@ -474,15 +431,9 @@ int main (int argc, char ** argv) MrhsProjector.Allocate(nbasis,FrbGrid,Coarse5d); MrhsProjector.ImportBasis(Aggregates.subspace); - TwoLevelADEF2mrhs - HDCGmrhs(1.0e-8, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhs_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 500, FineHermOp, HDCGmrhs_ADEF2, MrhsProjector.fine_grid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid); diff --git a/tests/debug/Test_general_coarse_hdcg_phys48_blockcg.cc b/tests/debug/Test_general_coarse_hdcg_phys48_blockcg.cc index 36bc6d13b..e9cde1399 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys48_blockcg.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys48_blockcg.cc @@ -156,21 +156,6 @@ void LoadEigenvectors(std::vector &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; template class FixedCGPolynomial : public LinearFunction { @@ -267,28 +252,6 @@ public: } }; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; @@ -615,15 +578,9 @@ int main (int argc, char ** argv) std::cout << "**************************************"< - HDCGmrhs(1.0e-8, 300, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhs_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 300, FineHermOp, HDCGmrhs_ADEF2, MrhsProjector.fine_grid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid); @@ -667,15 +624,9 @@ int main (int argc, char ** argv) for(auto tol : tols) { - TwoLevelADEF2mrhs - HDCGmrhsSloppy(tol, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhsSloppy_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhsSloppy(tol, 500, FineHermOp, HDCGmrhsSloppy_ADEF2, MrhsProjector.fine_grid); // Solve again to 10^-5 for(int r=0;r &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; template class FixedCGPolynomial : public LinearFunction { @@ -267,28 +252,6 @@ public: } }; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; @@ -614,15 +577,9 @@ int main (int argc, char ** argv) std::cout << "**************************************"< - HDCGmrhs(1.0e-8, 300, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhs_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 300, FineHermOp, HDCGmrhs_ADEF2, MrhsProjector.fine_grid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid); @@ -647,15 +604,9 @@ int main (int argc, char ** argv) for(auto tol : tols) { - TwoLevelADEF2mrhs - HDCGmrhsSloppy(tol, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhsSloppy_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhsSloppy(tol, 500, FineHermOp, HDCGmrhsSloppy_ADEF2, MrhsProjector.fine_grid); // Solve again to 10^-5 for(int r=0;r &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; int main (int argc, char ** argv) diff --git a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc index fb806bf2c..e57fcc050 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc @@ -100,44 +100,7 @@ void LoadEigenvectors(std::vector &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; int main (int argc, char ** argv) @@ -348,15 +311,9 @@ int main (int argc, char ** argv) MrhsProjector.Allocate(nbasis,FrbGrid,Coarse5d); MrhsProjector.ImportBasis(Aggregates.subspace); - TwoLevelADEF2mrhs - HDCGmrhs(1.0e-8, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhs_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 500, FineHermOp, HDCGmrhs_ADEF2, MrhsProjector.fine_grid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid); diff --git a/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc b/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc index abfb0afe7..f525d3fda 100644 --- a/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc +++ b/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc @@ -100,44 +100,7 @@ void LoadEigenvectors(std::vector &eval, #endif } -// Want Op in CoarsenOp to call MatPcDagMatPc -template -class HermOpAdaptor : public LinearOperatorBase -{ - LinearOperatorBase & wrapped; -public: - HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; - void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); } - void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); } - void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,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 HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } -}; -template class CGSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - int iters; - CGSmoother(int _iters, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - iters(_iters) - { - std::cout << GridLogMessage<<" Mirs smoother order "< CG(0.0,iters,false); // non-converge is just fine in a smoother - - out=Zero(); - - CG(_SmootherOperator,in,out); - } -}; int main (int argc, char ** argv) @@ -356,15 +319,9 @@ int main (int argc, char ** argv) CGSmoother CGsmooth(ord,ShiftedFineHermOp) ; - TwoLevelADEF2mrhs - HDCGmrhs(1.0e-8, 500, - FineHermOp, - CGsmooth, - HPDSolveMrhs, // Used in M1 - HPDSolveMrhs, // Used in Vstart - MrhsProjector, - MrhsGuesser, - CoarseMrhs); + MrhsADEF2Preconditioner + HDCGmrhs_ADEF2(FineHermOp, CGsmooth, HPDSolveMrhs, HPDSolveMrhs, MrhsProjector, MrhsGuesser, CoarseMrhs); + TwoLevelCGmrhs HDCGmrhs(1.0e-8, 500, FineHermOp, HDCGmrhs_ADEF2, MrhsProjector.fine_grid); std::vector src_mrhs(nrhs,FrbGrid); std::vector res_mrhs(nrhs,FrbGrid);