portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-08-22 18:39:36 +01:00
Code Releases Activity

Ready for faster 3 level solve with dense coarse and BLAS on SRHS !

Browse Source 
Exciting
This commit is contained in:
Peter Boyle
2026-08-20 11:44:26 -04:00
parent 2fde4f3e50
commit 97c9178785
8 changed files with 573 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CLAUDE.md
+2
View File
@@ -174,6 +174,8 @@ Every programme is wrapped in `Grid_init(&argc, &argv)` / `Grid_finalize()` (`Gr
- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`. - Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`.
- **Tensor indices are positional, not labelled.** The `Grid/tensors/` arithmetic recurses structurally over the `iScalar`/`iVector`/`iMatrix` nest: each level defines only the {scalar,vector,matrix}² products at its own level, with element types resolved by automatic type deduction, so every colour/spin/lorentz combination composes from ~200 lines (versus the pre-C++11 QDP++/PETE approach of machine-generating every case). An index's meaning derives entirely from its nesting depth counted from the outside; `iScalar` is the identity/broadcast case at every level. Never insert or remove a nesting level casually — the multiplication tables contract by position. - **Tensor indices are positional, not labelled.** The `Grid/tensors/` arithmetic recurses structurally over the `iScalar`/`iVector`/`iMatrix` nest: each level defines only the {scalar,vector,matrix}² products at its own level, with element types resolved by automatic type deduction, so every colour/spin/lorentz combination composes from ~200 lines (versus the pre-C++11 QDP++/PETE approach of machine-generating every case). An index's meaning derives entirely from its nesting depth counted from the outside; `iScalar` is the identity/broadcast case at every level. Never insert or remove a nesting level casually — the multiplication tables contract by position.
- **Multigrid coarsening deepens the tensor nest by one level.** A coarse site vector is `iVector<CComplex,nbasis>`, and `innerProduct` on it returns `iScalar<CComplex>` — one level deeper than the fine block scalar. So the block-inner-product scalar type gains one `iScalar` wrapper per MG level (fine: `vTComplex`; level 2: `iScalar<vTComplex>`; see `examples/Example_pvdagm_3level.cc`). When calling `blockInnerProduct`/`blockZAXPY`/`blockOrthogonalise` on coarse fields, the coarse scalar type must match `decltype(innerProduct(siteVector(),siteVector()))` exactly; a wrong depth fails to compile (no viable `operator=` deep in the instantiation chain) rather than mis-contracting. - **Multigrid coarsening deepens the tensor nest by one level.** A coarse site vector is `iVector<CComplex,nbasis>`, and `innerProduct` on it returns `iScalar<CComplex>` — one level deeper than the fine block scalar. So the block-inner-product scalar type gains one `iScalar` wrapper per MG level (fine: `vTComplex`; level 2: `iScalar<vTComplex>`; see `examples/Example_pvdagm_3level.cc`). When calling `blockInnerProduct`/`blockZAXPY`/`blockOrthogonalise` on coarse fields, the coarse scalar type must match `decltype(innerProduct(siteVector(),siteVector()))` exactly; a wrong depth fails to compile (no viable `operator=` deep in the instantiation chain) rather than mis-contracting.
- **Grids are borrowed, never owned.** `conformable` is pointer identity, so every object that interoperates must hold the *same* `GridCartesian *`; a class that minted its own grid internally could never conform with anything else. Ownership is therefore not available, and lifetime is managed by scope discipline instead of reference counting: whoever creates a grid retains it beyond every object it handed a reference to. Anything *derived* from a grid inherits this — `~PaddedCell` dereferences its `unpadded_grid`, so a `PaddedCell` cannot even be **destroyed** after its parent grid, only used. Where a consumer must let go early, it offers an explicit hand-back (`MultiGeneralCoarsenedOperatorV2::ReleaseGrid()`) to be called *before* the grid is destroyed.
- The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`. - The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`.
- Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI. - Use `GRID_ASSERT(cond)` (defined in `Grid/GridStd.h`), not bare `assert` — it prints a Grid-formatted message and aborts cleanly under MPI.
- Logging is stream-based, not macro-based: `std::cout << GridLogMessage << ... << std::endl;`. Channels declared in `Grid/log/Log.h` include `GridLogError`, `GridLogWarning`, `GridLogDebug`, `GridLogPerformance`, `GridLogIterative`, `GridLogSolver`, `GridLogHMC`, `GridLogComms`, `GridLogMemory`, `GridLogDslash`, `GridLogIRL`, `GridLogMG`. A subset is switched on at runtime with e.g. `--log Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours` (names given without the `GridLog` prefix). - Logging is stream-based, not macro-based: `std::cout << GridLogMessage << ... << std::endl;`. Channels declared in `Grid/log/Log.h` include `GridLogError`, `GridLogWarning`, `GridLogDebug`, `GridLogPerformance`, `GridLogIterative`, `GridLogSolver`, `GridLogHMC`, `GridLogComms`, `GridLogMemory`, `GridLogDslash`, `GridLogIRL`, `GridLogMG`. A subset is switched on at runtime with e.g. `--log Error,Warning,Message,Performance,Iterative,Integrator,Debug,Colours` (names given without the `GridLog` prefix).
Grid/algorithms/multigrid/DenseCoarseMatrix.h
+115 -45
View File
@@ -81,19 +81,25 @@ NAMESPACE_BEGIN(Grid);
// Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD // Tensor-depth agnostic: site scalar objects treated as contiguous ComplexD
// (iScalar wrappers add no data), so any MG level's coarse operator imports. // (iScalar wrappers add no data), so any MG level's coarse operator imports.
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
template<class Fobj,class CComplex,int nbasis> //
class DenseCoarseMatrix : public LinearFunction<typename GeneralCoarsenedMatrix<Fobj,CComplex,nbasis>::CoarseVector> { // Depends on a coarse operator only to extract its matrix elements; thereafter
// it is given a coarse vector and applies the inverse. Import() is a template
// member so any of the coarse classes will do, and the type of a
// DenseCoarseMatrix does not record which one built it.
//
template<class CComplex,int nbasis>
class DenseCoarseMatrix : public LinearFunction<Lattice<iVector<CComplex,nbasis> > > {
public: public:
typedef GeneralCoarsenedMatrix<Fobj,CComplex,nbasis> GeneralCoarseOp; typedef iVector<CComplex,nbasis > siteVector;
typedef typename GeneralCoarseOp::CoarseVector Field; typedef Lattice<siteVector> CoarseVector;
typedef typename GeneralCoarseOp::CoarseMatrix CoarseMatrix; typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
typedef CoarseVector Field;
using LinearFunction<Field>::operator(); using LinearFunction<Field>::operator();
typedef typename Field::vector_object vobj; typedef typename Field::vector_object vobj;
typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_object sobj;
typedef typename CoarseMatrix::vector_object Mvobj; typedef typename CoarseMatrix::vector_object Mvobj;
typedef typename Mvobj::scalar_object Msobj; typedef typename Mvobj::scalar_object Msobj;
GeneralCoarseOp &_Op; // the coarse operator: stencil source + certificate oracle
GridBase *grid; GridBase *grid;
int nd; int nd;
int64_t N; // dense rank = gSites * nbasis int64_t N; // dense rank = gSites * nbasis
@@ -121,12 +127,11 @@ public:
int devSum; int devSum;
double schurAuditRel; // DENSE_SCHUR=2: rel slab diff single-vs-schur (-1 = not run) double schurAuditRel; // DENSE_SCHUR=2: rel slab diff single-vs-schur (-1 = not run)
DenseCoarseMatrix(GeneralCoarseOp &Op, GridBase *g) DenseCoarseMatrix(GridBase *g)
: _Op(Op), grid(g) : grid(g)
{ {
GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(ComplexD) ); GRID_ASSERT( sizeof(sobj) == nbasis*sizeof(ComplexD) );
GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(ComplexD) ); GRID_ASSERT( sizeof(Msobj) == nbasis*nbasis*sizeof(ComplexD) );
GRID_ASSERT( grid == Op.Grid() );
nd = grid->_ndimension; nd = grid->_ndimension;
N = grid->gSites() * nbasis; N = grid->gSites() * nbasis;
lsites = grid->lSites(); lsites = grid->lSites();
@@ -158,7 +163,16 @@ public:
} }
slab.resize((uint64_t)nrows * N); slab.resize((uint64_t)nrows * N);
}
////////////////////////////////////////////////////////////////////
// The only place a coarse operator is needed: pull its elements, invert,
// and make the slab resident. Any class exposing Geometry() and
// ExtractMatrix(p,A) will do -- single RHS or either multiRHS.
////////////////////////////////////////////////////////////////////
template<class CoarseOp>
void Import(CoarseOp &Op)
{
double t0 = usecond(); double t0 = usecond();
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
// 0. Slab cache: SLAB_FILE=<stem> -> per-rank raw file <stem>.<rank>. // 0. Slab cache: SLAB_FILE=<stem> -> per-rank raw file <stem>.<rank>.
@@ -187,9 +201,9 @@ public:
} }
} }
if (!loaded) { if (!loaded) {
ImportDense(); // slab <- my rows of A (LOCAL, no comms) ImportDense(Op); // slab <- my rows of A (LOCAL, no comms)
ImportCertificate(); // dense apply == Op.M on a non-constant vector ImportCertificate(Op); // dense apply == Op.M, before inversion
InvertDense(); // slab <- my rows of A^{-1} InvertDense(Op); // slab <- my rows of A^{-1}
double t1 = usecond(); double t1 = usecond();
if (sfile) { if (sfile) {
FILE *f = fopen(slabfile.c_str(),"wb"); FILE *f = fopen(slabfile.c_str(),"wb");
@@ -251,7 +265,7 @@ public:
double ta = usecond(); double ta = usecond();
(*this)(x, y); (*this)(x, y);
double tb = usecond(); double tb = usecond();
_Op.M(y, z); ApplyOracle(Op, y, z);
z = z - x; z = z - x;
RealD rel = std::sqrt(norm2(z)/norm2(x)); RealD rel = std::sqrt(norm2(z)/norm2(x));
std::cout << GridLogMessage << "DenseCoarseMatrix: VERIFY ||A Ainv x - x||/||x|| = " std::cout << GridLogMessage << "DenseCoarseMatrix: VERIFY ||A Ainv x - x||/||x|| = "
@@ -262,11 +276,52 @@ public:
<< (usecond()-t0)/1.0e6 << " s" << std::endl; << (usecond()-t0)/1.0e6 << " s" << std::endl;
} }
////////////////////////////////////////////////////////////////////
// Apply the source operator to a D dimensional field, whichever kind it is.
//
// A multiRHS operator lives on the D+1 grid, so drive it with several right
// hand sides at once: slice r carries (r+1)*in, and linearity says the
// results must scale likewise. One apply, and unlike a single rhs check it
// also catches rhs mixing. Cheap check, not a production path.
////////////////////////////////////////////////////////////////////
template<class CoarseOp>
void ApplyOracle(CoarseOp &Op,const Field &in, Field &out)
{
if ( Op.Grid() == grid ) { Op.M(in,out); return; }
GridBase *mgrid = Op.Grid();
GRID_ASSERT(mgrid->_ndimension == nd+1);
int nr = mgrid->_fdimensions[0];
Field min(mgrid), mout(mgrid);
for(int r=0;r<nr;r++){
Field scaled(grid);
scaled = ComplexD(r+1.0,0.0)*in;
InsertSliceFast(scaled,min,r,0);
}
Op.M(min,mout);
ExtractSliceFast(out,mout,0,0);
for(int r=1;r<nr;r++){
Field sr(grid),d(grid);
ExtractSliceFast(sr,mout,r,0);
d = sr - ComplexD(r+1.0,0.0)*out;
RealD rel = std::sqrt(norm2(d)/norm2(sr));
if ( rel >= 1.0e-6 ) {
std::cout << GridLogMessage << "DenseCoarseMatrix: oracle rhs "<<r
<<" inconsistent with rhs 0, rel "<<rel<<std::endl;
}
GRID_ASSERT( rel < 1.0e-6 );
}
}
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
// 1. Direct stencil -> dense import of MY ROWS of A (no comms): // 1. Direct stencil -> dense import of MY ROWS of A (no comms):
// Dense[(s,a),(wrap(s+shift_p),b)] += A[p][s]_{a,b} // Dense[(s,a),(wrap(s+shift_p),b)] += A[p][s]_{a,b}
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
void ImportDense(void) template<class CoarseOp>
void ImportDense(CoarseOp &Op)
{ {
double t = -usecond(); double t = -usecond();
Coordinate gdims = grid->GlobalDimensions(); Coordinate gdims = grid->GlobalDimensions();
@@ -276,12 +331,12 @@ public:
uint64_t nelem = (uint64_t)nrows * N; uint64_t nelem = (uint64_t)nrows * N;
thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); }); thread_for(i, nelem, { slab[i] = ComplexF(0.0,0.0); });
for(int p=0; p<_Op.geom.npoint; p++){ for(int p=0; p<Op.Geometry().npoint; p++){
Coordinate shift = _Op.geom.shifts[p]; Coordinate shift = Op.Geometry().shifts[p];
// _A[p] is PADDED after ExchangeCoarseLinks (end of CoarsenOperator): // _A[p] is PADDED after ExchangeCoarseLinks (end of CoarsenOperator):
// extract the unpadded field before peeking with unpadded coordinates // extract the unpadded field before peeking with unpadded coordinates
// (exactly as MultiGeneralCoarsenedMatrix::CopyMatrix does). // (exactly as MultiGeneralCoarsenedMatrix::CopyMatrix does).
CoarseMatrix Aun = _Op.Cell.Extract(_Op._A[p]); CoarseMatrix Aun(grid); Op.ExtractMatrix(p,Aun);
if ( getenv("DENSE_IMPORT_DEBUG") ) { if ( getenv("DENSE_IMPORT_DEBUG") ) {
// Peek-path vs field-norm audit: sum |peekLocalSite|^2 must match norm2 // Peek-path vs field-norm audit: sum |peekLocalSite|^2 must match norm2
double pk = 0.0; double pk = 0.0;
@@ -295,7 +350,7 @@ public:
RealD gpk = pk; RealD gpk = pk;
grid->GlobalSumVector(&gpk, 1); grid->GlobalSumVector(&gpk, 1);
std::cout << GridLogMessage << "DenseCoarseMatrix: DEBUG p=" << p std::cout << GridLogMessage << "DenseCoarseMatrix: DEBUG p=" << p
<< " norm2(_A[p]) " << norm2(_Op._A[p]) << " norm2(_A[p]) " << norm2(Aun)
<< " norm2(Extract) " << norm2(Aun) << " norm2(Extract) " << norm2(Aun)
<< " sum|peek|^2 " << gpk << std::endl; << " sum|peek|^2 " << gpk << std::endl;
} }
@@ -346,7 +401,7 @@ public:
grid->GlobalSumVector(&gz, 1); grid->GlobalSumVector(&gz, 1);
std::cout << GridLogMessage << "DenseCoarseMatrix: stencil->dense import took " std::cout << GridLogMessage << "DenseCoarseMatrix: stencil->dense import took "
<< t/1.0e6 << " s (" << _Op.geom.npoint << " points, local, no comms)" << t/1.0e6 << " s (" << Op.Geometry().npoint << " points, local, no comms)"
<< " zero rows " << (int64_t)gz << "/" << N << std::endl; << " zero rows " << (int64_t)gz << "/" << N << std::endl;
// Debug: coordinate pattern of live sites (mechanism fingerprint) // Debug: coordinate pattern of live sites (mechanism fingerprint)
@@ -375,7 +430,8 @@ public:
// 2. IMPORT CERTIFICATE: dense rows vs Op.M on a NON-CONSTANT vector. // 2. IMPORT CERTIFICATE: dense rows vs Op.M on a NON-CONSTANT vector.
// (Constant x has x[s+d]==x[s-d]: blind to a shift-sign error.) // (Constant x has x[s+d]==x[s-d]: blind to a shift-sign error.)
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
void ImportCertificate(void) template<class CoarseOp>
void ImportCertificate(CoarseOp &Op)
{ {
Field x(grid); Field Ax(grid); Field Dx(grid); Field x(grid); Field Ax(grid); Field Dx(grid);
for(int ss=0; ss<lsites; ss++){ for(int ss=0; ss<lsites; ss++){
@@ -406,7 +462,7 @@ public:
for(int b=0; b<nbasis; b++) ((ComplexD *)&s)[b] = yh[ss*nbasis+b]; for(int b=0; b<nbasis; b++) ((ComplexD *)&s)[b] = yh[ss*nbasis+b];
pokeLocalSite(s, Dx, myLcoor[ss]); pokeLocalSite(s, Dx, myLcoor[ss]);
} }
_Op.M(x, Ax); ApplyOracle(Op, x, Ax);
Field d(grid); d = Dx - Ax; Field d(grid); d = Dx - Ax;
RealD rel = std::sqrt(norm2(d)/norm2(Ax)); RealD rel = std::sqrt(norm2(d)/norm2(Ax));
std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE ||Dense x - A x||/||A x|| = " std::cout << GridLogMessage << "DenseCoarseMatrix: IMPORT CERTIFICATE ||Dense x - A x||/||A x|| = "
@@ -428,7 +484,8 @@ public:
// slab difference; keep the Schur result // slab difference; keep the Schur result
// (so VERIFY certifies the new path). // (so VERIFY certifies the new path).
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
void InvertDense(void) template<class CoarseOp>
void InvertDense(CoarseOp &Op)
{ {
char *sc = getenv("DENSE_SCHUR"); char *sc = getenv("DENSE_SCHUR");
int mode = sc ? atoi(sc) : 0; int mode = sc ? atoi(sc) : 0;
@@ -440,7 +497,7 @@ public:
} }
if ( mode == 1 ) if ( mode == 1 )
{ {
InvertDenseSchur(); InvertDenseSchur(Op);
return; return;
} }
GRID_ASSERT( mode == 2 ); GRID_ASSERT( mode == 2 );
@@ -448,7 +505,7 @@ public:
InvertDenseSingle(); InvertDenseSingle();
std::vector<ComplexF> ref(slab); // Ainv, single path std::vector<ComplexF> ref(slab); // Ainv, single path
slab = Aimp; slab = Aimp;
InvertDenseSchur(); // slab = Ainv, Schur path InvertDenseSchur(Op); // slab = Ainv, Schur path
// NaN-PROOF comparison: max() masks NaN, so count non-finite // NaN-PROOF comparison: max() masks NaN, so count non-finite
// entries in each result explicitly. // entries in each result explicitly.
@@ -713,17 +770,18 @@ public:
// different precision order). NaN-proof: non-finite entries are // different precision order). NaN-proof: non-finite entries are
// counted explicitly since max() silently masks NaN. // counted explicitly since max() silently masks NaN.
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
void ImportDenseFP64(BlockRows &S, std::vector<int64_t> &g2rm) template<class CoarseOp>
void ImportDenseFP64(CoarseOp &Op, BlockRows &S, std::vector<int64_t> &g2rm)
{ {
Coordinate gdims = grid->GlobalDimensions(); Coordinate gdims = grid->GlobalDimensions();
int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1; int sign = getenv("DENSE_IMPORT_SIGN") ? atoi(getenv("DENSE_IMPORT_SIGN")) : 1;
GRID_ASSERT( sign==1 || sign==-1 ); GRID_ASSERT( sign==1 || sign==-1 );
std::vector<ComplexD> h((uint64_t)nrows*N, ComplexD(0.0,0.0)); std::vector<ComplexD> h((uint64_t)nrows*N, ComplexD(0.0,0.0));
for(int p=0; p<_Op.geom.npoint; p++) for(int p=0; p<Op.Geometry().npoint; p++)
{ {
Coordinate shift = _Op.geom.shifts[p]; Coordinate shift = Op.Geometry().shifts[p];
CoarseMatrix Aun = _Op.Cell.Extract(_Op._A[p]); CoarseMatrix Aun(grid); Op.ExtractMatrix(p,Aun);
autoView(Av, Aun, CpuRead); autoView(Av, Aun, CpuRead);
thread_for(ss, lsites, { thread_for(ss, lsites, {
Coordinate ncoor(nd); Coordinate ncoor(nd);
@@ -785,7 +843,8 @@ public:
// slab. Everything downstream (device residency, split-K apply, // slab. Everything downstream (device residency, split-K apply,
// VERIFY, SLAB_FILE) is untouched. // VERIFY, SLAB_FILE) is untouched.
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
void InvertDenseSchur(void) template<class CoarseOp>
void InvertDenseSchur(CoarseOp &Op)
{ {
double t1 = usecond(); double t1 = usecond();
int P = grid->ProcessorCount(); int P = grid->ProcessorCount();
@@ -812,7 +871,7 @@ public:
} }
BlockRows S; BlockRows S;
ImportDenseFP64(S, g2rm); ImportDenseFP64(Op, S, g2rm);
int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES")) int64_t panelBytes = getenv("DENSE_PANEL_BYTES") ? atol(getenv("DENSE_PANEL_BYTES"))
: (int64_t)1024*1024*1024; : (int64_t)1024*1024*1024;
@@ -921,13 +980,21 @@ public:
((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]); ((ComplexD *)&s)[b] = ComplexD(hY[ss*nbasis + b]);
pokeLocalSite(s, psi, myLcoor[ss]); pokeLocalSite(s, psi, myLcoor[ss]);
} }
if ( getenv("DENSE_CC_CHECK") ) { }
Field tmp(grid);
_Op.M(psi, tmp); ////////////////////////////////////////////////////////////////////
tmp = tmp - src; // Defect of an applied inverse. Was a DENSE_CC_CHECK block inside
std::cout << GridLogMessage << "DenseCoarseMatrix: apply defect ||A x - b||/||b|| = " // operator(), but that is virtual and cannot take an operator, so the
<< std::sqrt(norm2(tmp)/norm2(src)) << std::endl; // caller now asks for it explicitly.
} ////////////////////////////////////////////////////////////////////
template<class CoarseOp>
void CheckApply(CoarseOp &Op,const Field &src,const Field &psi)
{
Field tmp(grid);
Op.M(psi, tmp);
tmp = tmp - src;
std::cout << GridLogMessage << "DenseCoarseMatrix: apply defect ||A x - b||/||b|| = "
<< std::sqrt(norm2(tmp)/norm2(src)) << std::endl;
} }
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
@@ -960,14 +1027,17 @@ public:
double t1 = usecond(); double t1 = usecond();
std::cout << GridLogMessage << "DenseCoarseMatrix: batched apply " << nr << " rhs took " std::cout << GridLogMessage << "DenseCoarseMatrix: batched apply " << nr << " rhs took "
<< (t1-t0)/1000.0 << " ms (" << (t1-t0)/1000.0/nr << " ms/rhs)" << std::endl; << (t1-t0)/1000.0 << " ms (" << (t1-t0)/1000.0/nr << " ms/rhs)" << std::endl;
if ( getenv("DENSE_CC_CHECK") ) { }
Field tmp(grid);
for(int rr=0; rr<nr; rr++){ template<class CoarseOp>
_Op.M(psi[rr], tmp); void CheckApplyBatch(CoarseOp &Op,std::vector<Field> &src,std::vector<Field> &psi,int nr)
tmp = tmp - src[rr]; {
std::cout << GridLogMessage << "DenseCoarseMatrix: batch defect["<<rr<<"] = " Field tmp(grid);
<< std::sqrt(norm2(tmp)/norm2(src[rr])) << std::endl; for(int rr=0; rr<nr; rr++){
} Op.M(psi[rr], tmp);
tmp = tmp - src[rr];
std::cout << GridLogMessage << "DenseCoarseMatrix: batch defect["<<rr<<"] = "
<< std::sqrt(norm2(tmp)/norm2(src[rr])) << std::endl;
} }
} }
Grid/algorithms/multigrid/GeneralCoarsenedMatrix.h
+11
View File
@@ -74,6 +74,17 @@ public:
/////////////////////// ///////////////////////
// Interface // 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 * 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 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 GridCartesian * CoarseGrid(void) { return _CoarseGrid; }; // this is all the linalg routines need to know
Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h
+12 -3
View File
@@ -77,13 +77,22 @@ public:
GridBase * Grid(void) { return _CoarseGridMulti; }; // this is all the linalg routines need to know 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 GridCartesian * CoarseGrid(void) { return _CoarseGridMulti; }; // this is all the linalg routines need to know
// Can be used to do I/O on the operator matrices externally //////////////////////////////////////////////////////////////////////////
void SetMatrix (int p,CoarseMatrix & A) // 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<CoarseMatrix> & A)
{ {
GRID_ASSERT(A.size()==geom_srhs.npoint); GRID_ASSERT(A.size()==geom_srhs.npoint);
GridtoBLAS(A[p],BLAS_A[p]); GridtoBLAS(A[p],BLAS_A[p]);
} }
void GetMatrix (int p,CoarseMatrix & A) void GetMatrix (int p,std::vector<CoarseMatrix> & A)
{ {
GRID_ASSERT(A.size()==geom_srhs.npoint); GRID_ASSERT(A.size()==geom_srhs.npoint);
BLAStoGrid(A[p],BLAS_A[p]); BLAStoGrid(A[p],BLAS_A[p]);
Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h
+12 -3
View File
@@ -110,13 +110,22 @@ public:
} }
} }
// Can be used to do I/O on the operator matrices externally //////////////////////////////////////////////////////////////////////////
void SetMatrix (int p,CoarseMatrix & A) // Bilingual accessors, matching GeneralCoarsenedMatrix. Note Grid() is the
// D+1 multiRHS grid here, so a consumer wanting the space the elements live
// on must ask for CoarseGridD().
//////////////////////////////////////////////////////////////////////////
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<CoarseMatrix> & A)
{ {
GRID_ASSERT(A.size()==geom_srhs.npoint); GRID_ASSERT(A.size()==geom_srhs.npoint);
GridtoBLAS(A[p],BLAS_A[p]); GridtoBLAS(A[p],BLAS_A[p]);
} }
void GetMatrix (int p,CoarseMatrix & A) void GetMatrix (int p,std::vector<CoarseMatrix> & A)
{ {
GRID_ASSERT(A.size()==geom_srhs.npoint); GRID_ASSERT(A.size()==geom_srhs.npoint);
BLAStoGrid(A[p],BLAS_A[p]); BLAStoGrid(A[p],BLAS_A[p]);
examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc
+3 -2
View File
@@ -467,7 +467,7 @@ int main (int argc, char ** argv)
typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis> SubspaceL2; typedef Aggregation<CoarseSiteObj,vTTComplex,nbasis> SubspaceL2;
// The library dense bottom over the L2 coarse operator // The library dense bottom over the L2 coarse operator
typedef DenseCoarseMatrix<CoarseSiteObj,vTTComplex,nbasis> DenseCC_t; typedef DenseCoarseMatrix<vTTComplex,nbasis> DenseCC_t;
PVdagM_t PVdagM(Ddwf,Dpv); PVdagM_t PVdagM(Ddwf,Dpv);
ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv); ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
@@ -573,7 +573,8 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "**********************************************" << std::endl; std::cout << GridLogMessage << "**********************************************" << std::endl;
std::cout << GridLogMessage << " Dense CC inverse setup (library DenseCoarseMatrix)" << std::endl; std::cout << GridLogMessage << " Dense CC inverse setup (library DenseCoarseMatrix)" << std::endl;
std::cout << GridLogMessage << "**********************************************" << std::endl; std::cout << GridLogMessage << "**********************************************" << std::endl;
DenseCC.reset(new DenseCC_t(LittleDiracOpL2, CoarseCoarse5d)); DenseCC.reset(new DenseCC_t(CoarseCoarse5d));
DenseCC->Import(LittleDiracOpL2);
MrhsDenseCC.reset(new MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector>(*DenseCC, CoarseCoarse5d, nrhs)); MrhsDenseCC.reset(new MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector>(*DenseCC, CoarseCoarse5d, nrhs));
} }
examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc
+387 -24
View File
@@ -29,8 +29,8 @@ Author: Peter Boyle <pboyle@bnl.gov>
// //
// PVdagM three level multigrid on the V2 coarse operator. // PVdagM three level multigrid on the V2 coarse operator.
// //
// STAGE ONE: grids, types, subspace, and the L1 coarsening only. The L2 // STAGES ONE AND TWO: grids, types, subspace, and the L1 and L2 coarsenings.
// chain, the dense bottom and the solves are not here yet. // The dense bottom and the solves are not here yet.
// //
// Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc: // Differences from Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc:
// //
@@ -53,14 +53,15 @@ Author: Peter Boyle <pboyle@bnl.gov>
// aggregation gives e_k, and the near null content is silently gone. The // aggregation gives e_k, and the near null content is silently gone. The
// ||<psi|psi> - I||_F guard below is what catches that. // ||<psi|psi> - I||_F guard below is what catches that.
// //
// Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK COARSEN_BATCH // Env: LATT LS MASS NBASIS(compile time) NRHS BLOCK BLOCK2 COARSEN_BATCH
// HOT_START CONFIG SUBSPACE_FILE V1_CHECK // HOT_START CONFIG SUBSPACE_FILE V1_CHECK MRHS_COARSEN
// //
#include <Grid/Grid.h> #include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h> #include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h> #include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h> #include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
#include <memory> #include <memory>
@@ -78,12 +79,32 @@ int Ls = 24;
int CoarsenBatch = 9; int CoarsenBatch = 9;
std::vector<int> lat_size({48,48,48,96}); std::vector<int> lat_size({48,48,48,96});
// Solver tuning, values as in the V1 example
RealD FineSmootherShift = 0.1;
int FineSmootherOrder = 16;
RealD CoarseSmootherShift = 0.1;
int CoarseSmootherNstep = 4;
RealD CoarseSolverTol = 0.03;
int CoarseSolverOrder = 200;
RealD OuterTol = 1.0e-8;
int OuterMmax = 8;
int OuterNstep = 8;
void ParseEnvironment(void) void ParseEnvironment(void)
{ {
if(getenv("MASS")) mass = atof(getenv("MASS")); if(getenv("MASS")) mass = atof(getenv("MASS"));
if(getenv("NRHS")) Nrhs = atoi(getenv("NRHS")); if(getenv("NRHS")) Nrhs = atoi(getenv("NRHS"));
if(getenv("LS")) Ls = atoi(getenv("LS")); if(getenv("LS")) Ls = atoi(getenv("LS"));
if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH")); if(getenv("COARSEN_BATCH")) CoarsenBatch= atoi(getenv("COARSEN_BATCH"));
if(getenv("FineSmootherShift")) FineSmootherShift = atof(getenv("FineSmootherShift"));
if(getenv("FineSmootherOrder")) FineSmootherOrder = atoi(getenv("FineSmootherOrder"));
if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol"));
if(getenv("CoarseSolverOrder")) CoarseSolverOrder = atoi(getenv("CoarseSolverOrder"));
if(getenv("OuterTol")) OuterTol = atof(getenv("OuterTol"));
if(getenv("OuterMmax")) OuterMmax = atoi(getenv("OuterMmax"));
if(getenv("OuterNstep")) OuterNstep = atoi(getenv("OuterNstep"));
if(getenv("LATT")){ if(getenv("LATT")){
Coordinate l; Coordinate l;
GridCmdOptionIntVector(std::string(getenv("LATT")),l); GridCmdOptionIntVector(std::string(getenv("LATT")),l);
@@ -139,9 +160,8 @@ public:
}; };
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
// ||<v|v> - I||_F over a set of coarse vectors. ~0.23 means the raw near // ||<v|v> - I||_F over a set of coarse vectors. Small means the raw near null
// null content survived the projection; ~sqrt(N_sites) means block // content survived the projection; see GramGuard for where a leak lands.
// orthonormal vectors leaked in and every image collapsed to e_k.
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
template<class CoarseField> template<class CoarseField>
RealD GramDefect(std::vector<CoarseField> &v) RealD GramDefect(std::vector<CoarseField> &v)
@@ -157,6 +177,228 @@ RealD GramDefect(std::vector<CoarseField> &v)
return std::sqrt(s2); return std::sqrt(s2);
} }
// On a leak every image collapses to the block unit e_k, the Gram becomes
// N*I, and the defect lands at (N-1)*sqrt(nbasis) -- orders above the ~0.2
// of a content preserving projection. Trip well below that so a mis-set
// threshold costs a log line rather than the run.
template<class CoarseField>
void GramGuard(const std::string &name,std::vector<CoarseField> &v,GridBase *grid)
{
RealD defect = GramDefect(v);
RealD N = (RealD)grid->gSites();
RealD leak = (N-1.0)*std::sqrt((RealD)v.size());
RealD trip = std::sqrt(N);
std::cout << GridLogMessage << "GUARD: ||<"<<name<<"|"<<name<<"> - I||_F = " << defect
<< " (e_k leak would be " << leak << ", trip at " << trip << ")" << std::endl;
GRID_ASSERT( defect < trip );
}
//////////////////////////////////////////////////////////////////////
// Shifted variants for the smoothers
//////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat; Matrix &_PV;
public:
RealD shift;
ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){ Field tmp(in.Grid()); _Mat.M(in,tmp); _PV.Mdag(tmp,out); out = out + shift*in; }
void AdjOp (const Field &in, Field &out){ Field tmp(in.Grid()); _PV.M(tmp,out); _Mat.Mdag(in,tmp); out = out + shift*in; }
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){ Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
};
template<class Field>
class ShiftedLinearOperator : public LinearOperatorBase<Field> {
LinearOperatorBase<Field> &_Op; RealD shift;
public:
ShiftedLinearOperator(RealD _shift, LinearOperatorBase<Field> &Op) : _Op(Op), shift(_shift) {}
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out) { assert(0); }
void Op (const Field &in, Field &out) { _Op.Op(in,out); out = out + shift*in; }
void AdjOp (const Field &in, Field &out) { _Op.AdjOp(in,out); out = out + shift*in; }
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp (const Field &in, Field &out) { Field tmp(in.Grid()); Op(in,tmp); AdjOp(tmp,out); }
};
//////////////////////////////////////////////////////////////////////
// Dense L3 solve on the packed D+1 coarse-coarse field
//////////////////////////////////////////////////////////////////////
template<class DenseType, class CoarseCoarseField>
class MrhsDenseCCSolve : public LinearFunction<CoarseCoarseField> {
public:
DenseType &_Dense;
int _nrhs;
MrhsDenseCCSolve(DenseType &D, int nrhs) : _Dense(D), _nrhs(nrhs) {}
using LinearFunction<CoarseCoarseField>::operator();
virtual void operator()(const CoarseCoarseField &in, CoarseCoarseField &out){
_Dense.ApplyBatch6D(in, out, _nrhs);
}
};
//////////////////////////////////////////////////////////////////////
// mrhs interfaces + single-polynomial mrhs PGCR
//////////////////////////////////////////////////////////////////////
template<class Field>
class MrhsLinearFunction {
public:
virtual void operator()(std::vector<Field> &in, std::vector<Field> &out) = 0;
};
template<class Field>
class MrhsPGCRNonHermitian {
public:
RealD Tolerance; Integer MaxIterations; int mmax,nstep,steps,level;
int ZeroGuess = 0; int FirstCycle = 0;
std::string name = "Level 1";
LinearOperatorBase<Field> &Linop;
MrhsLinearFunction<Field> &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<Field> &_Linop,MrhsLinearFunction<Field> &Prec,int _mmax,int _nstep)
: Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; }
static RealD vnorm2(std::vector<Field> &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; }
static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &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<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); }
void vOp(std::vector<Field> &in,std::vector<Field> &out){ for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
void operator()(std::vector<Field> &src,std::vector<Field> &psi){
RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq;
std::vector<Field> r(nrhs,grid);
GridStopWatch T; T.Start(); steps=0; FirstCycle=1;
for(int k=0;k<MaxIterations;k++){
cp=GCRnStep(src,psi,rsq);
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR("<<mmax<<","<<nstep<<") "<<steps<<" steps cp = "<<cp<<" target "<<rsq<<std::endl;
if(cp<rsq){
T.Stop(); vOp(psi,r); for(int rr=0;rr<nrhs;rr++) axpy(r[rr],-1.0,src[rr],r[rr]);
RealD tr=vnorm2(r);
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR: Converged on iteration "<<steps
<<" computed residual "<<std::sqrt(cp/ssq)<<" true residual "<<std::sqrt(tr/ssq)<<" target "<<Tolerance<<std::endl;
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR Time elapsed: Total "<<T.Elapsed()<<std::endl;
return;
}
}
std::cout<<GridLogMessage<<"MrhsPGCR: did not converge"<<std::endl;
}
RealD GCRnStep(std::vector<Field> &src,std::vector<Field> &psi,RealD rsq){
RealD cp; ComplexD a,b,rq; RealD zAAz; int nrhs=src.size(); GridBase *grid=src[0].Grid();
std::vector<Field> r(nrhs,grid),z(nrhs,grid),Az(nrhs,grid);
std::vector< std::vector<Field> > q(mmax,std::vector<Field>(nrhs,grid));
std::vector< std::vector<Field> > p(mmax,std::vector<Field>(nrhs,grid));
std::vector<RealD> qq(mmax);
if (ZeroGuess && FirstCycle) { for(int rr=0;rr<nrhs;rr++){ psi[rr]=Zero(); r[rr]=src[rr]; } }
else { vOp(psi,Az); for(int rr=0;rr<nrhs;rr++) r[rr]=src[rr]-Az[rr]; }
FirstCycle=0;
Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
p[0]=z; q[0]=Az; qq[0]=zAAz; cp=vnorm2(r);
for(int k=0;k<nstep;k++){
steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"] resid "<<cp<<" target "<<rsq<<std::endl;
if((k==nstep-1)||(cp<rsq)) return cp;
Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
q[peri_kp]=Az; p[peri_kp]=z;
int northog=((kp)>(mmax-1))?(mmax-1):(kp);
for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-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 on the D+1 coarse field
//////////////////////////////////////////////////////////////////////
template<class CoarseField, class CoarseCoarseField>
class MrhsCoarseThreeLevelPrec : public LinearFunction<CoarseField> {
public:
LinearOperatorBase<CoarseField> &_CoarseOp;
LinearFunction<CoarseField> &_CoarseSmoother;
MultiRHSBlockProject<CoarseField> &_Projector;
LinearFunction<CoarseCoarseField> &_CoarseCoarseSolve;
GridBase *_Coarse5d, *_CoarseCoarse5d, *_CoarseCoarseMrhs;
int _nrhs;
MrhsCoarseThreeLevelPrec(LinearOperatorBase<CoarseField> &CoarseOp,
LinearFunction<CoarseField> &CoarseSmoother,
MultiRHSBlockProject<CoarseField> &Projector,
LinearFunction<CoarseCoarseField> &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<CoarseField>::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();
_CoarseCoarseSolve(CCsrc,CCsol);
_Projector.blockPromote(vec1,CCsol);
add(out,out,vec1);
_CoarseOp.Op(out,vec1); sub(vec1,in,vec1);
vec2=Zero();
_CoarseSmoother(vec1,vec2);
add(out,out,vec2);
}
};
//////////////////////////////////////////////////////////////////////
// L1->L2 mrhs V-cycle
//////////////////////////////////////////////////////////////////////
template<class FineField, class MrhsCoarseVector, class FineSmoother>
class MrhsTwoLevelMG : public MrhsLinearFunction<FineField> {
public:
typedef MrhsCoarseVector CoarseVector;
LinearOperatorBase<FineField> &_FineOperator;
FineSmoother &_PostSmoother;
MultiRHSBlockProject<FineField> &_Projector;
LinearFunction<CoarseVector> &_CoarseSolve;
GridBase *_CoarseGrid, *_CoarseGridMrhs;
MrhsTwoLevelMG(LinearOperatorBase<FineField> &FineOp, FineSmoother &Post,
MultiRHSBlockProject<FineField> &Projector, LinearFunction<CoarseVector> &CoarseSolve,
GridBase *CoarseGrid, GridBase *CoarseGridMrhs)
: _FineOperator(FineOp),_PostSmoother(Post),_Projector(Projector),_CoarseSolve(CoarseSolve),
_CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){}
virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){
int nrhs=in.size(); GridBase *fgrid=in[0].Grid();
std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid);
for(int r=0;r<nrhs;r++) out[r]=in[r];
for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
// fine vector -> D+1 coarse, via the mixed blockProject
CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs);
_Projector.blockProject(vec1,CsrcMrhs);
CsolMrhs=Zero();
_CoarseSolve(CsrcMrhs,CsolMrhs);
_Projector.blockPromote(vec1,CsolMrhs);
for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); }
}
};
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
@@ -233,7 +475,8 @@ int main (int argc, char ** argv)
MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c); MobiusFermionD Dpv (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0, M5,b,c);
typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t; typedef PVdagMLinearOperator<MobiusFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<MobiusFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv); PVdagM_t PVdagM(Ddwf,Dpv);
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
@@ -310,13 +553,7 @@ int main (int argc, char ** argv)
MrhsProjector.blockProject(rawNull,psi_coarse); // RAW vectors in MrhsProjector.blockProject(rawNull,psi_coarse); // RAW vectors in
rawNull.clear(); rawNull.shrink_to_fit(); rawNull.clear(); rawNull.shrink_to_fit();
{ GramGuard("psi_coarse",psi_coarse,Coarse5d);
RealD defect = GramDefect(psi_coarse);
RealD leak = std::sqrt((double)Coarse5d->gSites());
std::cout << GridLogMessage << "GUARD: ||<psi_coarse|psi_coarse> - I||_F = " << defect
<< " (~0.23 good; ~sqrt(N_coarse)=" << leak << " = e_k leak)" << std::endl;
GRID_ASSERT( defect < leak );
}
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
// Optional cross check of the coarse matrix elements against the V1 // Optional cross check of the coarse matrix elements against the V1
@@ -366,9 +603,9 @@ int main (int argc, char ** argv)
ComplexD *w2=(ComplexD *)&h2[0]; ComplexD *w2=(ComplexD *)&h2[0];
int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD); int64_t words = sites*sizeof(calcMatrix)/sizeof(ComplexD);
for(int64_t i=0;i<words;i++){ for(int64_t i=0;i<words;i++){
ComplexD d=w1[i]-w2[i]; ComplexD d=w1[i]-w2[i];
num += real(d)*real(d)+imag(d)*imag(d); num += real(d)*real(d)+imag(d)*imag(d);
den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]); den += real(w1[i])*real(w1[i])+imag(w1[i])*imag(w1[i]);
} }
} }
std::cout << GridLogMessage << "V1_CHECK: |A_V1|^2 = " << den << std::endl; std::cout << GridLogMessage << "V1_CHECK: |A_V1|^2 = " << den << std::endl;
@@ -441,11 +678,7 @@ int main (int argc, char ** argv)
std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d); std::vector<CoarseCoarseVector> psi_cc(nbasis,CoarseCoarse5d);
MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in MrhsProjectorL2.blockProject(rawPsi,psi_cc); // RAW vectors in
RealD defect = GramDefect(psi_cc); GramGuard("psi_cc",psi_cc,CoarseCoarse5d);
RealD leak = std::sqrt((double)CoarseCoarse5d->gSites());
std::cout << GridLogMessage << "GUARD: ||<psi_cc|psi_cc> - I||_F = " << defect
<< " (~0.23 good; ~sqrt(N_cc)=" << leak << " = e_k leak)" << std::endl;
GRID_ASSERT( defect < leak );
} }
rawPsi.clear(); rawPsi.shrink_to_fit(); rawPsi.clear(); rawPsi.shrink_to_fit();
@@ -470,7 +703,137 @@ int main (int argc, char ** argv)
GRID_ASSERT( norm2(ccout) > 0.0 ); GRID_ASSERT( norm2(ccout) > 0.0 );
} }
std::cout << GridLogMessage << "*** stage two complete: L1 and L2 coarse operators built ***" << std::endl; //////////////////////////////////////////////////////////////////////
// STAGE THREE (part one): the dense bottom on L2.
//
// DenseCoarseMatrix is bilingual: it takes the elements through
// Geometry()/ExtractMatrix(), so the V2 operator serves directly. It does
// detect that a multiRHS op cannot apply on the D dimensional grid and
// skips its own certificate and VERIFY, so the equivalent check is done
// here instead, driving the L2 operator at Nrhs 1 through a slice.
//////////////////////////////////////////////////////////////////////
typedef DenseCoarseMatrix<CComplexS2,nbasis> DenseCC_t;
std::unique_ptr<DenseCC_t> DenseCC;
if ( getenv("DENSE_CC")==nullptr || atoi(getenv("DENSE_CC")) ) {
std::cout << GridLogMessage << "*** L3 dense bottom: import from the V2 L2 operator ***" << std::endl;
DenseCC.reset(new DenseCC_t(CoarseCoarse5d));
DenseCC->Import(CoarseOpL2);
////////////////////////////////////////////////////////////////////
// ||A Ainv x - x|| / ||x||, the check Import could not run itself
////////////////////////////////////////////////////////////////////
Coordinate cc1latt({1,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CoarseCoarseOne = new GridCartesian(cc1latt,cmsimd,cmmpi);
CoarseOpL2.SetGrid(CoarseCoarseOne);
CoarseCoarseVector x(CoarseCoarse5d),y(CoarseCoarse5d),z(CoarseCoarse5d);
GridParallelRNG dRNG(CoarseCoarse5d); dRNG.SeedFixedIntegers({11,12,13,14});
random(dRNG,x);
(*DenseCC)(x,y); // y = Ainv x
CoarseCoarseVector y1(CoarseCoarseOne),z1(CoarseCoarseOne);
InsertSliceFast(y,y1,0,0);
CoarseOpL2.M(y1,z1); // z = A y
ExtractSliceFast(z,z1,0,0);
z = z - x;
RealD rel = std::sqrt(norm2(z)/norm2(x));
std::cout << GridLogMessage << "L3 dense: ||A Ainv x - x||/||x|| = " << rel << std::endl;
GRID_ASSERT( rel < 1.0e-2 );
CoarseOpL2.SetGrid(CoarseCoarseMrhs);
delete CoarseCoarseOne;
}
//////////////////////////////////////////////////////////////////////
// STAGE THREE (part two): the solves.
//
// Both operators are driven from the SAME objects at whatever Nrhs is
// asked for -- the matrix elements were built once and survive SetGrid --
// so single RHS and multiRHS are the same code path with a different grid.
//////////////////////////////////////////////////////////////////////
GRID_ASSERT(DenseCC != nullptr); // the PGCR bottom is not ported yet
typedef PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> FineSmoother_t;
ShiftedPVdagM_t ShiftedPVdagM(FineSmootherShift,Ddwf,Dpv);
TrivialPrecon<LatticeFermionD> simple_fine;
TrivialPrecon<CoarseVector> simpleC;
auto RunSolve = [&](int nr)
{
std::cout << GridLogMessage << "**********************************************" << std::endl;
std::cout << GridLogMessage << " V2 THREE-level solve, Nrhs = " << nr << std::endl;
std::cout << GridLogMessage << "**********************************************" << std::endl;
Coordinate cml({nr,1,clatt[0],clatt[1],clatt[2],clatt[3]});
Coordinate ccml({nr,1,cclatt[0],cclatt[1],cclatt[2],cclatt[3]});
GridCartesian *CMrhs = new GridCartesian(cml, cmsimd,cmmpi);
GridCartesian *CCMrhs = new GridCartesian(ccml,cmsimd,cmmpi);
CoarseOpPV.SetGrid(CMrhs);
CoarseOpL2.SetGrid(CCMrhs);
NonHermitianLinearOperator<CoarseOperator,CoarseVector> LinOpC (CoarseOpPV);
NonHermitianLinearOperator<CoarseCoarseOperator,CoarseCoarseVector> LinOpCC(CoarseOpL2);
MrhsDenseCCSolve<DenseCC_t,CoarseCoarseVector> ccSolve(*DenseCC,nr);
ShiftedLinearOperator<CoarseVector> ShiftedC(CoarseSmootherShift, LinOpC);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherNstep,CoarseSmootherNstep);
CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
L2to3Precon(LinOpC, CoarseSmootherGCR, MrhsProjectorL2, ccSolve,
Coarse5d, CoarseCoarse5d, CCMrhs, nr);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
L2PGCR(CoarseSolverTol, CoarseSolverOrder/16, LinOpC, L2to3Precon, 16, 16);
L2PGCR.Level(2); L2PGCR.Name("Couter"); L2PGCR.SetZeroGuess(1);
FineSmoother_t SmootherGCR(0.0,1,ShiftedPVdagM,simple_fine,FineSmootherOrder,FineSmootherOrder);
SmootherGCR.Level(1); SmootherGCR.Name("Fsmoother"); SmootherGCR.SetZeroGuess(1);
MrhsTwoLevelMG<LatticeFermionD,CoarseVector,FineSmoother_t>
ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
MrhsPGCRNonHermitian<LatticeFermionD>
L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1);
std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
GridStopWatch w; w.Start();
L1PGCR(src,sol);
w.Stop();
std::cout << GridLogMessage << "V2 3-level solve Nrhs "<<nr<<" total " << w.Elapsed()
<< " (per RHS: " << w.useconds()/1.0e6/nr << " s)" << std::endl;
{ LatticeFermionD Ax(FGrid); RealD worst=0.0;
for(int r=0;r<nr;r++){ PVdagM.Op(sol[r],Ax); Ax=Ax-src[r];
RealD rn=std::sqrt(norm2(Ax)/norm2(src[r]));
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": rhs["<<r<<"] true residual = " << rn << std::endl;
worst=std::max(worst,rn); }
std::cout << GridLogMessage << "FINAL Nrhs "<<nr<<": worst-case residual = " << worst << std::endl;
}
// The operators borrow these grids and build a PaddedCell on them, so
// they must let go before the grids are destroyed.
CoarseOpPV.ReleaseGrid();
CoarseOpL2.ReleaseGrid();
delete CMrhs; delete CCMrhs;
};
RunSolve(nrhs);
if ( getenv("SOLVE_SRHS")==nullptr || atoi(getenv("SOLVE_SRHS")) ) RunSolve(1);
std::cout << GridLogMessage << "*** stage three complete: solves done ***" << std::endl;
Grid_finalize(); Grid_finalize();
} }
tests/debug/Test_coarse_v2_coarsen.cc
+31
View File
@@ -223,6 +223,37 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "|A_V1 - A_V2srhs|^2 / |A_V1|^2 = " << nums/den << std::endl; std::cout << GridLogMessage << "|A_V1 - A_V2srhs|^2 / |A_V1|^2 = " << nums/den << std::endl;
GRID_ASSERT( nums/den < 1.0e-20 ); GRID_ASSERT( nums/den < 1.0e-20 );
////////////////////////////////////////////////
// GetMatrix/SetMatrix round trip through the BLAS layout array. This is
// how a bilingual DenseCoarseMatrix will retrieve the elements, and it
// needs no SetGrid since the matrix elements are Nrhs independent.
////////////////////////////////////////////////
{
typedef MrhsV2::CoarseMatrix CoarseMatrixS;
std::vector<CoarseMatrixS> Aget(npoint,CoarseS);
for(int p=0;p<npoint;p++) OpV2.GetMatrix(p,Aget);
MrhsV2 OpV2rt(geomS,CoarseS);
for(int p=0;p<npoint;p++) OpV2rt.SetMatrix(p,Aget);
std::vector<std::vector<calcMatrix> > A4;
ReadMatrix(OpV2rt,npoint,A4);
RealD numrt=0.0;
for(int p=0;p<npoint;p++){
GRID_ASSERT(A2[p].size()==A4[p].size());
ComplexD *w2 = (ComplexD *)&A2[p][0];
ComplexD *w4 = (ComplexD *)&A4[p][0];
int64_t words = A2[p].size()*sizeof(calcMatrix)/sizeof(ComplexD);
for(int64_t i=0;i<words;i++){
ComplexD d = w2[i]-w4[i];
numrt += real(d)*real(d)+imag(d)*imag(d);
}
}
std::cout << GridLogMessage << "GetMatrix/SetMatrix round trip |diff|^2 = " << numrt << std::endl;
GRID_ASSERT( numrt == 0.0 );
}
//////////////////////////////////////////////// ////////////////////////////////////////////////
// and V2 applies the matrix it just built // and V2 applies the matrix it just built
//////////////////////////////////////////////// ////////////////////////////////////////////////