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

First test of distributed schur recursive inverse

Browse Source
This commit is contained in:
Peter Boyle
2026-08-14 17:53:35 -04:00
parent 55c064de28
commit c22473f15d
2 changed files with 1052 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/debug/Test_schur_dense_coarse.cc
+289
View File
@@ -0,0 +1,289 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Test_schur_dense_coarse.cc
Copyright (C) 2026
Author: Peter Boyle <pboyle@bnl.gov>
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 */
//
// T6 of the RecursiveSchurInverse regression chain
// (schur_recursive_inverse_plan.txt 4B.5): the DenseCoarseMatrix GLUE,
// on a real (tiny) lattice coarse operator, CPU laptop build.
//
// Builds a genuine GeneralCoarsenedMatrix (DWF MdagM + 0.5 shift for a
// guaranteed-invertible Galerkin coarse op, random aggregation basis,
// nbasis=8, 4^4 x Ls/1 blocking) and constructs DenseCoarseMatrix in
// DENSE_SCHUR=2 AUDIT mode with small DENSE_PANEL_BYTES (multi-panel
// gathers exercised through the glue). The constructor then runs, in
// order, all the certificates this stage exists to check:
// - fresh ImportDense (no SLAB_FILE) + IMPORT CERTIFICATE vs Op.M
// - InvertDenseSingle (the oracle)
// - InvertDenseSchur: self-certifying rank-major map, fp64 diagonal
// import certificate vs the fp32 slab, distributed recursion,
// growth telemetry
// - AUDIT: max|Ainv_schur - Ainv_single| over the full slab
// - VERIFY ||A Ainv x - x||/||x|| through the SCHUR result
// This program adds asserts on the audit number and a random-vector
// round trip.
//
// Uniform local volume 12.12.12.12 (fine), per-dim blocks {4,4,3,3},
// coarse 3.3.4.4/rank, nbasis 4 (N = 576n):
// mpirun -n 1 ./Test_schur_dense_coarse --grid 12.12.12.12 --mpi 1.1.1.1
// mpirun -n 2 ./Test_schur_dense_coarse --grid 12.12.12.24 --mpi 1.1.1.2
// mpirun -n 3 ./Test_schur_dense_coarse --grid 12.12.12.36 --mpi 1.1.1.3
// mpirun -n 4 ./Test_schur_dense_coarse --grid 12.12.12.48 --mpi 1.1.1.4
//
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/multigrid/DenseCoarseMatrix.h>
using namespace std;
using namespace Grid;
///////////////////////////////////////////////////////////////////////
// MdagM + shift: Galerkin projection of a PD operator plus sigma I is
// safely invertible whatever the (random) subspace quality.
///////////////////////////////////////////////////////////////////////
template<class Field>
class ShiftedHermOpAdaptor : public LinearOperatorBase<Field>
{
LinearOperatorBase<Field> &wrapped;
RealD shift;
public:
ShiftedHermOpAdaptor(LinearOperatorBase<Field> &wrapme, RealD s)
: wrapped(wrapme), shift(s) {};
void Op(const Field &in, Field &out)
{
wrapped.HermOp(in, out);
out = out + shift*in;
}
void AdjOp(const Field &in, Field &out)
{
Op(in, out);
}
void HermOp(const Field &in, Field &out)
{
Op(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<Field> &out) { GRID_ASSERT(0); }
void HermOpAndNorm(const Field &in, Field &out, RealD &n1, RealD &n2) { GRID_ASSERT(0); }
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls = 4;
const int nbasis = 4;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
GridDefaultSimd(Nd,vComplex::Nsimd()),
GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// Per-dimension blocking {4,4,3,3}: fine 12.12.12.12 -> coarse
// 3.3.4.4. Two constraints meet here (both MEASURED today):
// - coarse dims of 2 hit the probing pathology (health probe below)
// - GEN-simd lanes {1,1,2,2} must land on even coarse dims, so the
// odd production-like 3s go on the lane-free x,y axes (exactly the
// production [3,6,8,8] trick).
Coordinate blocks({4,4,3,3});
Coordinate clatt = GridDefaultLatt();
for(int d=0; d<clatt.size(); d++)
{
GRID_ASSERT( (clatt[d] % blocks[d]) == 0 );
clatt[d] = clatt[d]/blocks[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt,
GridDefaultSimd(Nd,vComplex::Nsimd()),
GridDefaultMpi());
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({9,10,11,12});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG CRNG(Coarse5d); CRNG.SeedFixedIntegers(cseeds);
LatticeGaugeField Umu(UGrid);
SU<Nc>::HotConfiguration(RNG4,Umu);
RealD mass = 0.1;
RealD M5 = 1.8;
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
ShiftedHermOpAdaptor<LatticeFermionD> HOA(HermDefOp, 0.5);
std::cout << GridLogMessage << "Building random aggregation space, nbasis " << nbasis << std::endl;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse5d,FGrid,0);
Aggregates.CreateSubspaceRandom(RNG5);
std::cout << GridLogMessage << "Coarsening shifted MdagM" << std::endl;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
LittleDiracOp.CoarsenOperator(HOA,Aggregates);
///////////////////////////////////////////////////////////////////////
// Operator health probes (independent of DenseCoarseMatrix).
//
// MEASURED PATHOLOGY, banked 2026-08-14: on coarse dims of 2 (fine
// 8.8.8.8, block 4 -> coarse 2.2.2.2) the coarsened operator is
// rank 16/128 with 112 zero ROWS (output support = 2 of 16 sites)
// and Hermiticity violation 0.17 -- the probing construction breaks
// on the size-2 torus. The import certificate cannot see this
// (dense and M share _A). Out of scope here; coarse dims >= 3.
//
// Cheap any-size probes: output support + Hermiticity via inner
// products on random vectors.
///////////////////////////////////////////////////////////////////////
{
CoarseVector px(Coarse5d);
CoarseVector py(Coarse5d);
CoarseVector Mx(Coarse5d);
CoarseVector My(Coarse5d);
random(CRNG, px);
random(CRNG, py);
LittleDiracOp.M(px, Mx);
LittleDiracOp.M(py, My);
ComplexD ip1 = innerProduct(py, Mx); // <y, M x>
ComplexD ip3 = innerProduct(px, My); // <x, M y>
RealD hermdev = abs(ip1 - conj(ip3)) / std::sqrt(norm2(Mx)*norm2(py));
RealD support = norm2(Mx) / norm2(px);
std::cout << GridLogMessage << "Operator health: ||Mx||^2/||x||^2 = " << support
<< " herm-dev " << hermdev << std::endl;
// Hermitian fine op => exactly Hermitian Galerkin coarse op.
// (A measured herm-dev of 1.3e-4 here was the CPU SIMT-lane
// CoarsenOperator bug -- fixed 2026-08-14, now 4e-15. A loud
// failure here means _A population is broken again.)
GRID_ASSERT( support > 1.0e-3 );
GRID_ASSERT( hermdev < 1.0e-10 );
}
///////////////////////////////////////////////////////////////////////
// Full-matrix conditioning probe at small N: dense columns by
// applying M to unit vectors, fp64 Eigen SVD.
///////////////////////////////////////////////////////////////////////
{
int64_t Nprobe = Coarse5d->gSites() * nbasis;
if ( Nprobe <= 700 )
{
Eigen::MatrixXcd eA(Nprobe, Nprobe);
CoarseVector e(Coarse5d);
CoarseVector Me(Coarse5d);
for(int64_t j=0; j<Nprobe; j++)
{
int64_t gsite = j / nbasis;
int b = j % nbasis;
e = Zero();
Coordinate gcoor(Coarse5d->_ndimension);
Lexicographic::CoorFromIndex(gcoor, gsite, Coarse5d->GlobalDimensions());
typedef typename CoarseVector::vector_object::scalar_object csobj;
csobj s;
s = Zero();
((ComplexD *)&s)[b] = ComplexD(1.0,0.0);
pokeSite(s, e, gcoor);
LittleDiracOp.M(e, Me);
for(int64_t i=0; i<Nprobe; i++)
{
int64_t gsi = i / nbasis;
int bi = i % nbasis;
Coordinate gci(Coarse5d->_ndimension);
Lexicographic::CoorFromIndex(gci, gsi, Coarse5d->GlobalDimensions());
csobj si;
peekSite(si, Me, gci);
eA(i,j) = ((ComplexD *)&si)[bi];
}
}
Eigen::JacobiSVD<Eigen::MatrixXcd> svd(eA);
double smax = svd.singularValues()(0);
double smin = svd.singularValues()(Nprobe-1);
int64_t rank = 0;
for(int64_t i=0; i<Nprobe; i++)
{
if ( svd.singularValues()(i) > 1.0e-10*smax ) rank++;
}
double herm = (eA - eA.adjoint()).cwiseAbs().maxCoeff();
int64_t zrows = 0;
int64_t zcols = 0;
for(int64_t i=0; i<Nprobe; i++)
{
if ( eA.row(i).cwiseAbs().maxCoeff() < 1.0e-12 ) zrows++;
if ( eA.col(i).cwiseAbs().maxCoeff() < 1.0e-12 ) zcols++;
}
std::cout << GridLogMessage << "Operator probe: N=" << Nprobe
<< " sigma_max " << smax
<< " sigma_min " << smin
<< " rank " << rank << "/" << Nprobe
<< " herm-dev " << herm
<< " zero rows/cols " << zrows << "/" << zcols
<< std::endl;
}
}
///////////////////////////////////////////////////////////////////////
// T6: AUDIT mode, fresh import, multi-panel gathers. The constructor
// runs every certificate in the chain (see banner).
///////////////////////////////////////////////////////////////////////
setenv("DENSE_SCHUR","2",1);
setenv("DENSE_PANEL_BYTES","65536",1);
unsetenv("SLAB_FILE");
typedef DenseCoarseMatrix<vSpinColourVector,vTComplex,nbasis> DenseCC;
DenseCC dcm(LittleDiracOp, Coarse5d);
std::cout << GridLogMessage << "T6 audit relative slab difference (schur vs single) = "
<< dcm.schurAuditRel << std::endl;
GRID_ASSERT( dcm.schurAuditRel >= 0.0 ); // audit actually ran
GRID_ASSERT( dcm.schurAuditRel < 1.0e-3 );
///////////////////////////////////////////////////////////////////////
// Random-vector round trip through the SCHUR inverse
///////////////////////////////////////////////////////////////////////
CoarseVector x(Coarse5d);
CoarseVector y(Coarse5d);
CoarseVector z(Coarse5d);
random(CRNG, x);
dcm(x, y);
LittleDiracOp.M(y, z);
z = z - x;
RealD rel = std::sqrt(norm2(z)/norm2(x));
std::cout << GridLogMessage << "T6 round trip ||A Ainv x - x||/||x|| (random x) = "
<< rel << std::endl;
GRID_ASSERT( rel < 1.0e-2 );
std::cout << GridLogMessage << "Test_schur_dense_coarse: T6 ALL PASS" << std::endl;
Grid_finalize();
}
tests/debug/Test_schur_inverse.cc
+763
View File
@@ -0,0 +1,763 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Test_schur_inverse.cc
Copyright (C) 2026
Author: Peter Boyle <pboyle@bnl.gov>
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 */
//
// Staged regression gate for RecursiveSchurInverse (distributed dense
// inversion by recursive Schur complement) -- the laptop-side certificate
// chain of schur_recursive_inverse_plan.txt section 4B.5. Runs on a
// CPU-only build (Eigen BLAS backends) under mpirun:
//
// mpirun -n 1 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.1
// mpirun -n 2 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.2
// mpirun -n 3 ./Test_schur_inverse --grid 8.8.8.12 --mpi 1.1.1.3
// mpirun -n 4 ./Test_schur_inverse --grid 8.8.8.8 --mpi 1.1.1.4
//
// (n=3 exercises uneven row splits throughout.) The lattice exists only
// to furnish the communicator; no field is ever constructed.
//
// PRECISION: the inversion runs ENTIRELY in fp64 (decision 2026-08-14,
// superseding the fp32-merge design); certificates are eps64-scaled.
// The single terminal fp32 rounding belongs to the caller (tested at the
// glue level, Test_schur_dense_coarse).
//
// Stages present (cumulative -- earlier tests are never removed):
// T1a : ownership tables -- CheckRowStart on synthetic uneven partitions,
// MakeRowStart allgather vs closed form on the live communicator.
// T1b : STORAGE-CONVENTION PIN -- column-major + ld + window-offset
// semantics fixed once via identity multiplies through the
// explicit-ld gemmBatched, on INTEGER-VALUED data so all three
// cases below are EXACT (values well within the mantissa):
// (1) alpha=1,beta=0 read from an input column window
// (2) alpha=-1,beta=1 accumulate (the S-formation case)
// (3) write INTO an output column window, neighbours untouched
// No later failure can be a transposition/convention ambiguity.
// T2 : GatherGemm vs naive fp64 oracle (owner sub-ranges, alpha-beta
// cases, tiny+huge panels, half-participation call shape).
// T3 : LeafInvert in-place residual certificate.
// T4 : full recursive Invert vs Eigen fp64 oracle, growth-scaled
// certification, adversarial near-singular-A11 family with
// telemetry-spike assertion.
//
// Hard asserts throughout; thresholds pre-registered in the plan.
//
#include <Grid/Grid.h>
#include <Grid/Grid_Eigen_Dense.h>
#include <Grid/algorithms/multigrid/RecursiveSchurInverse.h>
using namespace std;
using namespace Grid;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
GridCartesian Comm(GridDefaultLatt(),
GridDefaultSimd(Nd,vComplex::Nsimd()),
GridDefaultMpi());
GridBase *grid = &Comm;
////////////////////////////////////////////////////////////////
// T1a : ownership tables
////////////////////////////////////////////////////////////////
{
// Synthetic partitions of N=97 (prime: every P>1 is uneven)
const int64_t N = 97;
for(int P=1; P<=4; P++)
{
std::vector<int64_t> table(P+1);
table[0] = 0;
for(int r=0; r<P; r++)
{
int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
table[r+1] = table[r] + nr;
}
RecursiveSchurInverse::CheckRowStart(table, N);
}
// Live allgather: deliberately uneven local counts, closed-form oracle
int P = grid->ProcessorCount();
int me = grid->ThisRank();
int64_t myNrows = 3 + me;
std::vector<int64_t> table = RecursiveSchurInverse::MakeRowStart(grid, myNrows);
std::vector<int64_t> expect(P+1);
expect[0] = 0;
for(int r=0; r<P; r++)
{
expect[r+1] = expect[r] + (3 + r);
}
GRID_ASSERT( (int)table.size() == P+1 );
for(int r=0; r<=P; r++)
{
GRID_ASSERT( table[r] == expect[r] );
}
// Constructor smoke: derived ownership matches
RecursiveSchurInverse RSI(grid, table[P], table, 1024*1024);
GRID_ASSERT( RSI.P == P );
GRID_ASSERT( RSI.me == me );
GRID_ASSERT( RSI.myRow0 == expect[me] );
GRID_ASSERT( RSI.myNrows == myNrows );
std::cout << GridLogMessage
<< "T1a ownership tables (synthetic P=1..4, live allgather, ctor) PASS"
<< std::endl;
}
////////////////////////////////////////////////////////////////
// T1b : storage-convention pin (every rank, local, exact)
////////////////////////////////////////////////////////////////
{
const int64_t rows = 5;
const int64_t cols = 13;
const int64_t col0 = 6; // input window start
const int64_t w = 4; // window width
// f(i,j): integer-valued, unique per element
auto f = [](int64_t i, int64_t j) -> ComplexD
{
return ComplexD( (RealD)(1 + i + 10*j), (RealD)(i - j) );
};
BlockRows A;
A.Resize(rows, cols);
{
std::vector<ComplexD> Ahost((uint64_t)rows*cols);
for(int64_t j=0; j<cols; j++)
{
for(int64_t i=0; i<rows; i++)
{
Ahost[(uint64_t)(i + j*rows)] = f(i,j);
}
}
acceleratorCopyToDevice(&Ahost[0], &A.data[0], (uint64_t)rows*cols*sizeof(ComplexD));
}
// Identity I_w, column major
deviceVector<ComplexD> Idev((uint64_t)w*w);
{
std::vector<ComplexD> Ihost((uint64_t)w*w, ComplexD(0.0,0.0));
for(int64_t d=0; d<w; d++)
{
Ihost[(uint64_t)(d + d*w)] = ComplexD(1.0,0.0);
}
acceleratorCopyToDevice(&Ihost[0], &Idev[0], (uint64_t)w*w*sizeof(ComplexD));
}
GridBLAS BLAS;
ComplexD one ( 1.0,0.0);
ComplexD minus (-1.0,0.0);
ComplexD zero ( 0.0,0.0);
deviceVector<ComplexD*> Ap(1);
deviceVector<ComplexD*> Bp(1);
deviceVector<ComplexD*> Cp(1);
std::vector<ComplexD*> ptr_h(1);
auto setptr = [&](deviceVector<ComplexD*> &d, ComplexD *p)
{
ptr_h[0] = p;
acceleratorCopyToDevice(&ptr_h[0], &d[0], sizeof(ComplexD*));
};
////////////////////////////////////////////////////////////
// Case 1: C = A(:, col0:col0+w) . I_w (alpha=1, beta=0)
////////////////////////////////////////////////////////////
{
deviceVector<ComplexD> Cdev((uint64_t)rows*w);
setptr(Ap, A.ColumnWindow(col0));
setptr(Bp, &Idev[0]);
setptr(Cp, &Cdev[0]);
BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
(int)rows, (int)w, (int)w,
one, Ap, (int)A.ld,
Bp, (int)w,
zero, Cp, (int)rows);
BLAS.synchronise();
std::vector<ComplexD> Chost((uint64_t)rows*w);
acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD));
for(int64_t j=0; j<w; j++)
{
for(int64_t i=0; i<rows; i++)
{
GRID_ASSERT( Chost[(uint64_t)(i + j*rows)] == f(i, col0+j) );
}
}
}
////////////////////////////////////////////////////////////
// Case 2: C = C0 - A(:, col0:col0+w) . I_w (alpha=-1, beta=1)
// -- the S-formation accumulate; exact on integer data
////////////////////////////////////////////////////////////
{
auto g = [](int64_t i, int64_t j) -> ComplexD
{
return ComplexD( (RealD)(100 + i + j), (RealD)7 );
};
deviceVector<ComplexD> Cdev((uint64_t)rows*w);
{
std::vector<ComplexD> Chost((uint64_t)rows*w);
for(int64_t j=0; j<w; j++)
{
for(int64_t i=0; i<rows; i++)
{
Chost[(uint64_t)(i + j*rows)] = g(i,j);
}
}
acceleratorCopyToDevice(&Chost[0], &Cdev[0], (uint64_t)rows*w*sizeof(ComplexD));
}
setptr(Ap, A.ColumnWindow(col0));
setptr(Bp, &Idev[0]);
setptr(Cp, &Cdev[0]);
BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
(int)rows, (int)w, (int)w,
minus, Ap, (int)A.ld,
Bp, (int)w,
one, Cp, (int)rows);
BLAS.synchronise();
std::vector<ComplexD> Chost((uint64_t)rows*w);
acceleratorCopyFromDevice(&Cdev[0], &Chost[0], (uint64_t)rows*w*sizeof(ComplexD));
for(int64_t j=0; j<w; j++)
{
for(int64_t i=0; i<rows; i++)
{
ComplexD expect = g(i,j) - f(i, col0+j);
GRID_ASSERT( Chost[(uint64_t)(i + j*rows)] == expect );
}
}
}
////////////////////////////////////////////////////////////
// Case 3: write INTO a column window of a wider C;
// columns outside the window must be untouched
////////////////////////////////////////////////////////////
{
const int64_t ccols = 6;
const int64_t cw0 = 2; // output window start
BlockRows C;
C.Resize(rows, ccols);
{
std::vector<ComplexD> Chost((uint64_t)rows*ccols, ComplexD(-999.0, 999.0));
acceleratorCopyToDevice(&Chost[0], &C.data[0], (uint64_t)rows*ccols*sizeof(ComplexD));
}
setptr(Ap, A.ColumnWindow(col0));
setptr(Bp, &Idev[0]);
setptr(Cp, C.ColumnWindow(cw0));
BLAS.gemmBatched(GridBLAS_OP_N, GridBLAS_OP_N,
(int)rows, (int)w, (int)w,
one, Ap, (int)A.ld,
Bp, (int)w,
zero, Cp, (int)C.ld);
BLAS.synchronise();
std::vector<ComplexD> Chost((uint64_t)rows*ccols);
acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)rows*ccols*sizeof(ComplexD));
for(int64_t j=0; j<ccols; j++)
{
for(int64_t i=0; i<rows; i++)
{
ComplexD got = Chost[(uint64_t)(i + j*rows)];
if ( (j >= cw0) && (j < cw0+w) )
{
GRID_ASSERT( got == f(i, col0 + (j-cw0)) );
}
else
{
GRID_ASSERT( got == ComplexD(-999.0, 999.0) );
}
}
}
}
std::cout << GridLogMessage
<< "T1b storage-convention pin (window read / S-accumulate / window write, exact) PASS"
<< std::endl;
}
////////////////////////////////////////////////////////////////
// T2 : GatherGemm vs naive double-precision oracle.
//
// Every rank generates the SAME full N x N random fp64 operands
// from a fixed seed (no comms needed for the oracle), keeps only
// its own rows in BlockRows form, and after each GatherGemm call
// checks its output window element-by-element against a plain
// triple-loop ComplexD accumulation over the same entries.
//
// Sweep: N in {8, 96, 97}; owner ranges full/upper-half/single;
// (alpha,beta) in {(1,0), (-1,1)}; panelBytes tiny (ragged
// many-chunk gathers) and huge (single panel). Sentinel columns
// outside the output window must be untouched. Finally, a
// HALF-PARTICIPATION case rehearses the recursion call pattern:
// lower ranks own B but pass EMPTY A/C (collectives only).
////////////////////////////////////////////////////////////////
{
int P = grid->ProcessorCount();
int me = grid->ThisRank();
std::mt19937 rng(777);
std::uniform_real_distribution<double> dist(-1.0,1.0);
const int64_t nout = 5; // output width
const int64_t colB = 3; // B window offset
const int64_t colC = 2; // C window offset
for(int64_t N : {8L, 96L, 97L})
{
// Ownership: uneven for any P not dividing N
std::vector<int64_t> table(P+1);
table[0] = 0;
for(int r=0; r<P; r++)
{
int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
table[r+1] = table[r] + nr;
}
int64_t r0 = table[me];
int64_t myNr = table[me+1] - table[me];
// Identical full operands on every rank
std::vector<ComplexD> Aglob((uint64_t)N*N);
std::vector<ComplexD> Bglob((uint64_t)N*N);
for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng));
for(uint64_t i=0; i<(uint64_t)N*N; i++) Bglob[i] = ComplexD(dist(rng),dist(rng));
// My rows of a full-matrix operand as a BlockRows
auto fillRows = [&](BlockRows &X, std::vector<ComplexD> &glob,
int64_t row0, int64_t nr)
{
X.Resize(nr, N);
if ( nr == 0 ) return;
std::vector<ComplexD> h((uint64_t)nr*N);
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<nr; i++)
{
h[(uint64_t)(i + j*nr)] = glob[(uint64_t)((row0+i) + j*N)];
}
}
acceleratorCopyToDevice(&h[0], &X.data[0], (uint64_t)nr*N*sizeof(ComplexD));
};
// Owner-range cases: full span, upper half, single interior rank
std::vector<std::pair<int,int> > ranges;
ranges.push_back(std::make_pair(0, P));
if ( P > 1 ) ranges.push_back(std::make_pair(P/2, P));
if ( P > 1 ) ranges.push_back(std::make_pair(1, 2));
for(auto range : ranges)
{
int rB0 = range.first;
int rB1 = range.second;
int64_t ka0 = table[rB0]; // A-column window start = B row span
int64_t k = table[rB1] - table[rB0];
for(int acase=0; acase<2; acase++)
{
ComplexD alpha = ( acase==0 ) ? ComplexD( 1.0,0.0) : ComplexD(-1.0,0.0);
ComplexD beta = ( acase==0 ) ? ComplexD( 0.0,0.0) : ComplexD( 1.0,0.0);
for(int64_t panelBytes : {64L, 1L<<30})
{
RecursiveSchurInverse RSI(grid, N, table, panelBytes);
BlockRows A;
BlockRows B;
BlockRows C;
fillRows(A, Aglob, r0, myNr);
fillRows(B, Bglob, r0, myNr);
// Output: sentinel-filled, window at colC
const ComplexD sentinel(-999.0, 999.0);
const int64_t ccols = colC + nout + 2;
C.Resize(myNr, ccols);
std::vector<ComplexD> C0((uint64_t)myNr*ccols, sentinel);
if ( acase == 1 )
{
// beta=1 needs defined window content: g(i,j), integer-valued
for(int64_t j=0; j<nout; j++)
{
for(int64_t i=0; i<myNr; i++)
{
C0[(uint64_t)(i + (colC+j)*myNr)] = ComplexD((RealD)(50+i+j), (RealD)-3);
}
}
}
if ( myNr > 0 )
{
acceleratorCopyToDevice(&C0[0], &C.data[0], (uint64_t)myNr*ccols*sizeof(ComplexD));
}
RSI.GatherGemm(alpha, A, ka0, k,
rB0, rB1,
B, colB, nout,
beta, C, colC);
std::vector<ComplexD> Chost((uint64_t)myNr*ccols);
if ( myNr > 0 )
{
acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*ccols*sizeof(ComplexD));
}
double tol = 1.0e-14 * (double)k;
for(int64_t j=0; j<ccols; j++)
{
for(int64_t i=0; i<myNr; i++)
{
ComplexD got = Chost[(uint64_t)(i + j*myNr)];
if ( (j >= colC) && (j < colC+nout) )
{
int64_t jj = j - colC;
ComplexD acc(0.0,0.0);
if ( acase == 1 )
{
acc = C0[(uint64_t)(i + j*myNr)];
}
for(int64_t t=0; t<k; t++)
{
acc += alpha
* Aglob[(uint64_t)((r0+i) + (ka0+t)*N)]
* Bglob[(uint64_t)((ka0+t) + (colB+jj)*N)];
}
GRID_ASSERT( abs(got - acc) < tol );
}
else
{
GRID_ASSERT( got == sentinel );
}
}
}
}
}
}
////////////////////////////////////////////////////////////
// Half-participation: owners = [0,ph) hold B; participants
// = [ph,P) hold A/C; owners pass EMPTY A/C and column
// offset 0 (collectives only) -- the recursion call shape.
////////////////////////////////////////////////////////////
if ( P > 1 )
{
int ph = ( P+1 ) / 2;
int64_t ka0 = table[0];
int64_t k = table[ph] - table[0];
int participant = ( me >= ph );
RecursiveSchurInverse RSI(grid, N, table, 64);
BlockRows A;
BlockRows B;
BlockRows C;
fillRows(B, Bglob, r0, myNr);
if ( participant )
{
fillRows(A, Aglob, r0, myNr);
C.Resize(myNr, nout);
}
ComplexD one (1.0,0.0);
ComplexD zero(0.0,0.0);
int64_t cA = participant ? ka0 : 0;
RSI.GatherGemm(one, A, cA, k,
0, ph,
B, colB, nout, // owners deposit from their B window
zero, C, 0);
if ( participant )
{
std::vector<ComplexD> Chost((uint64_t)myNr*nout);
acceleratorCopyFromDevice(&C.data[0], &Chost[0], (uint64_t)myNr*nout*sizeof(ComplexD));
double tol = 1.0e-14 * (double)k;
for(int64_t j=0; j<nout; j++)
{
for(int64_t i=0; i<myNr; i++)
{
ComplexD acc(0.0,0.0);
for(int64_t t=0; t<k; t++)
{
acc += Aglob[(uint64_t)((r0+i) + (ka0+t)*N)]
* Bglob[(uint64_t)((ka0+t) + (colB+j)*N)];
}
GRID_ASSERT( abs(Chost[(uint64_t)(i + j*myNr)] - acc) < tol );
}
}
}
}
}
std::cout << GridLogMessage
<< "T2 GatherGemm vs oracle (N=8/96/97, 3 owner ranges, 2 alpha-beta, tiny+huge panels, half-participation) PASS"
<< std::endl;
}
////////////////////////////////////////////////////////////////
// T3 : LeafInvert -- in-place fp64 inversion of the contiguous
// leaf window. Purely local, every rank runs its own
// uneven-size leaf; residual certificate in ComplexD.
////////////////////////////////////////////////////////////////
{
int P = grid->ProcessorCount();
int me = grid->ThisRank();
int64_t w = 17 + 3*me;
uint64_t len = (uint64_t)w*w;
std::vector<int64_t> table = RecursiveSchurInverse::MakeRowStart(grid, w);
RecursiveSchurInverse RSI(grid, table[P], table, 1<<20);
// A = w I + R : well conditioned
std::mt19937 rng(31 + me);
std::uniform_real_distribution<double> dist(-1.0,1.0);
std::vector<ComplexD> Ahost(len);
for(uint64_t i=0; i<len; i++) Ahost[i] = ComplexD(dist(rng),dist(rng));
for(int64_t d=0; d<w; d++) Ahost[(uint64_t)(d + d*w)] += ComplexD((RealD)w, 0.0);
BlockRows Ar;
Ar.Resize(w, w);
acceleratorCopyToDevice(&Ahost[0], &Ar.data[0], len*sizeof(ComplexD));
RSI.LeafInvert(0, w, Ar);
std::vector<ComplexD> X(len);
acceleratorCopyFromDevice(&Ar.data[0], &X[0], len*sizeof(ComplexD));
double maxdev = 0.0;
for(int64_t j=0; j<w; j++)
{
for(int64_t i=0; i<w; i++)
{
ComplexD acc(0.0,0.0);
for(int64_t t=0; t<w; t++)
{
acc += Ahost[(uint64_t)(i + t*w)] * X[(uint64_t)(t + j*w)];
}
if ( i==j ) acc -= ComplexD(1.0,0.0);
maxdev = std::max(maxdev, abs(acc));
}
}
GRID_ASSERT( maxdev < 1.0e-13 );
std::cout << GridLogMessage
<< "T3 LeafInvert in-place fp64 (residual " << maxdev << ") PASS" << std::endl;
}
////////////////////////////////////////////////////////////////
// T4 : full recursive Invert vs Eigen fp64 oracle.
//
// Every rank builds the SAME N x N fp64 matrix from a fixed seed,
// keeps its rows, inverts through the full SPMD recursion, then
// the test gathers the complete inverse (zero-fill GlobalSum) and
// checks BOTH certificates:
// cert1 = || A X - I ||_max (ComplexD accumulation)
// cert2 = max|X - Xref| / max|Xref| (Xref = Eigen fp64 inverse)
//
// Families (eps64-scaled tolerances; the fp32-era growth data
// rescales by eps64/eps32 ~ 1.9e-9):
// kappa-moderate : A = R + 3 sqrt(N) I
// kappa-large : A = R + 0.3 sqrt(N) I
// adversarial : leading block (rank 0's whole leaf) REPLACED by
// 1e-2 * (R' + 3 sqrt(b) I) inside a well-conditioned
// A -- the growth spike must REGISTER in telemetry
// (asserted > 10 when P > 1); at fp64 the certificate
// barely notices it: that insensitivity IS the point
// of the fp64 conversion.
//
// N=64 runs with panelBytes=128 (ragged many-chunk gathers inside
// the recursion); larger N with 1 MB panels.
////////////////////////////////////////////////////////////////
{
int P = grid->ProcessorCount();
int me = grid->ThisRank();
std::mt19937 rng(2026);
std::uniform_real_distribution<double> dist(-1.0,1.0);
for(int64_t N : {64L, 200L, 513L})
{
std::vector<int64_t> table(P+1);
table[0] = 0;
for(int r=0; r<P; r++)
{
int64_t nr = N/P + ( (r < (int)(N%P)) ? 1 : 0 );
table[r+1] = table[r] + nr;
}
int64_t r0 = table[me];
int64_t myNr = table[me+1] - table[me];
for(int fam=0; fam<3; fam++)
{
const char *famname = (fam==0) ? "kappa-moderate" :
(fam==1) ? "kappa-large" : "adversarial-A11";
double shift = (fam==1) ? 0.3*std::sqrt((double)N) : 3.0*std::sqrt((double)N);
double tol = (fam==0) ? 1.0e-12 :
(fam==1) ? 1.0e-11 : 5.0e-11;
// Identical operand on every rank (all draws rank-independent)
std::vector<ComplexD> Aglob((uint64_t)N*N);
for(uint64_t i=0; i<(uint64_t)N*N; i++) Aglob[i] = ComplexD(dist(rng),dist(rng));
for(int64_t d=0; d<N; d++) Aglob[(uint64_t)(d + d*N)] += ComplexD(shift, 0.0);
if ( fam == 2 )
{
// Leading block = rank 0's whole leaf, scaled down 100x but
// internally well conditioned (shift scales as sqrt(b): a
// FIXED shift makes A11 itself near-singular at large b).
int64_t b = ( P > 1 ) ? table[1] : N/4;
for(int64_t j=0; j<b; j++)
{
for(int64_t i=0; i<b; i++)
{
Aglob[(uint64_t)(i + j*N)] = ComplexD(0.01,0.0)*ComplexD(dist(rng),dist(rng));
}
}
RealD bshift = (RealD)(0.03*std::sqrt((double)b));
for(int64_t d=0; d<b; d++) Aglob[(uint64_t)(d + d*N)] += ComplexD(bshift,0.0);
}
// Eigen fp64 oracle
Eigen::MatrixXcd eA(N,N);
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<N; i++)
{
eA(i,j) = Aglob[(uint64_t)(i + j*N)];
}
}
Eigen::MatrixXcd Xref = eA.inverse();
// Distribute, invert
int64_t panelBytes = ( N == 64 ) ? 128 : (1<<20);
RecursiveSchurInverse RSI(grid, N, table, panelBytes);
BlockRows Arows;
Arows.Resize(myNr, N);
{
std::vector<ComplexD> h((uint64_t)myNr*N);
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<myNr; i++)
{
h[(uint64_t)(i + j*myNr)] = Aglob[(uint64_t)((r0+i) + j*N)];
}
}
acceleratorCopyToDevice(&h[0], &Arows.data[0], (uint64_t)myNr*N*sizeof(ComplexD));
}
RSI.Invert(Arows);
// Gather the full inverse: zero-fill + GlobalSum
std::vector<ComplexD> Xfull((uint64_t)N*N, ComplexD(0.0,0.0));
{
std::vector<ComplexD> h((uint64_t)myNr*N);
acceleratorCopyFromDevice(&Arows.data[0], &h[0], (uint64_t)myNr*N*sizeof(ComplexD));
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<myNr; i++)
{
Xfull[(uint64_t)((r0+i) + j*N)] = h[(uint64_t)(i + j*myNr)];
}
}
}
grid->GlobalSumVector(&Xfull[0], (int)(N*N));
// cert1 = ||A X - I||_max
double cert1 = 0.0;
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<N; i++)
{
ComplexD acc(0.0,0.0);
for(int64_t t=0; t<N; t++)
{
acc += Aglob[(uint64_t)(i + t*N)] * Xfull[(uint64_t)(t + j*N)];
}
if ( i==j ) acc -= ComplexD(1.0,0.0);
cert1 = std::max(cert1, abs(acc));
}
}
// cert2 = max|X - Xref| / max|Xref|
double maxref = 0.0;
double maxdif = 0.0;
for(int64_t j=0; j<N; j++)
{
for(int64_t i=0; i<N; i++)
{
maxref = std::max(maxref, abs(Xref(i,j)));
maxdif = std::max(maxdif, abs(Xfull[(uint64_t)(i + j*N)] - Xref(i,j)));
}
}
double cert2 = maxdif / maxref;
double maxNormB = 0.0;
for(uint64_t i=0; i<RSI.telNormB.size(); i++)
{
maxNormB = std::max(maxNormB, RSI.telNormB[i]);
}
// GROWTH-SCALED certification, eps64 (the fp32-era model with
// eps swapped: cert2 ~ (10-12) ||B||_F sqrt(N) eps; threshold =
// 3x margin, floored at the family tolerance). ||B||_F capped
// per family so growth cannot silently excuse a logic error.
// cert1 remains a loose absolute bound (an O(1) logic error
// gives cert1 ~ 1e2-1e3; fp64 rounding gives ~1e-10).
double eps64 = 2.3e-16;
double tolModel = 30.0 * std::max(1.0, maxNormB) * std::sqrt((double)N) * eps64;
double tolEff = std::max(tol, tolModel);
double capB = (fam==0) ? 100.0 : (fam==1) ? 2000.0 : 10000.0;
std::cout << GridLogMessage
<< "T4 N=" << N << " " << famname
<< " ||AX-I||_max " << cert1
<< " |X-Xref|/|Xref| " << cert2
<< " max||B||_F " << maxNormB
<< " tolEff " << tolEff
<< ( (cert2 < tolEff) && (cert1 < 1.0e-6) ? " PASS" : " FAIL" )
<< std::endl;
GRID_ASSERT( cert2 < tolEff );
GRID_ASSERT( cert1 < 1.0e-6 );
GRID_ASSERT( maxNormB < capB );
if ( (fam == 2) && (P > 1) )
{
GRID_ASSERT( maxNormB > 10.0 ); // the spike must REGISTER
}
}
}
std::cout << GridLogMessage
<< "T4 recursive Invert vs Eigen oracle (N=64/200/513, 3 families) PASS"
<< std::endl;
}
std::cout << GridLogMessage
<< "Test_schur_inverse: ALL STAGES PASS" << std::endl;
Grid_finalize();
}