Improvements for mulktigrid. + Schur setup improvements

Grid/lattice/Lattice_reduction.h

@@ -809,6 +809,21 @@ static void sliceInnerProductMatrix(  Eigen::MatrixXcd &mat, const Lattice<vobj>
 // Same code path for every Lattice<vobj>: fine fermion fields and coarse
 // multi-RHS fields (nrhs folded into the grid) alike.
 /////////////////////////////////////////////////////////////////////////////
+// The batch of views and coefficients is passed to the kernel BY VALUE as
+// lambda-captured kernel arguments (a ViewPack), so there is no per-call
+// deviceVector allocation and no synchronous host->device memcpy of pointer
+// tables: the only sync points are the reductions themselves.
+// LatticeView has no default constructor, so the pack holds raw aligned
+// storage and views are copied in bytewise (as basisRotateJ does through
+// acceleratorPut); the struct is trivially copyable as a kernel argument.
+template<class View,int B>
+struct ViewPack {
+  alignas(View) unsigned char raw[B*sizeof(View)];
+  ComplexD b[B];
+  accelerator_inline const View & v(int j) const { return reinterpret_cast<const View *>(raw)[j]; }
+  void set(int j,const View &view){ memcpy(raw+j*sizeof(View),&view,sizeof(View)); }
+};
+
 template<int B,class vobj>
 void rankInnerProductMultiChunk(ComplexD *out,int m,
 				const std::vector<const Lattice<vobj>*> &left,
@@ -822,14 +837,14 @@ void rankInnerProductMultiChunk(ComplexD *out,int m,
   GridBase *grid = right.Grid();
   const uint64_t sites = grid->oSites();
 
-  hostVector<View>   h_left_v(m);
+  std::vector<View> h_v; h_v.reserve(m);
-  deviceVector<View> d_left_v(m);
+  ViewPack<View,B> pack;
   for(int j=0;j<m;j++){
     conformable(*left[j],right);
-    h_left_v[j] = left[j]->View(AcceleratorRead);
+    h_v.push_back(left[j]->View(AcceleratorRead));
+    pack.set(j,h_v[j]);
   }
-  acceleratorCopyToDevice(&h_left_v[0],&d_left_v[0],m*sizeof(View));
+  for(int j=m;j<B;j++) pack.set(j,h_v[0]);   // valid but unused lanes
-  View *left_vp = &d_left_v[0];
 
   deviceVector<batch_t> partial(sites);
   batch_t *partial_v = &partial[0];
@@ -839,13 +854,13 @@ void rankInnerProductMultiChunk(ComplexD *out,int m,
 	auto r = right_v[ss];
 	batch_t acc;
 	for(int j=0;j<B;j++){
-	  if ( j<m ) acc._internal[j] = innerProductD(left_vp[j][ss],r);
+	  if ( j<m ) acc._internal[j] = innerProductD(pack.v(j)[ss],r);
 	  else       zeroit(acc._internal[j]);
 	}
 	partial_v[ss] = acc;
     });
   }
-  for(int j=0;j<m;j++) h_left_v[j].ViewClose();
+  for(int j=0;j<m;j++) h_v[j].ViewClose();
 
   auto res = sum(partial_v,sites);   // one reduction for the whole batch
   for(int j=0;j<m;j++) out[j] = TensorRemove(res._internal[j]);
@@ -877,48 +892,70 @@ void innerProductMulti(std::vector<ComplexD> &out,
   if ( out.size() ) right.Grid()->GlobalSumVector(&out[0],(int)out.size());
 }
 
-// z = z + sum_j b[j] x[j];  if do_norm, returns global |z|^2 from the same pass.
+// z = z + sum_{j<m} b[j] x[j] for one chunk of at most B vectors; if do_norm,
-template<class vobj>
+// also writes per-site |z|^2 into inner_tmp_v.
-RealD axpyMultiNormImpl(Lattice<vobj> &z,const std::vector<ComplexD> &b,
+template<int B,class vobj>
-			const std::vector<const Lattice<vobj>*> &x,int do_norm)
+void axpyMultiChunk(Lattice<vobj> &z,const ComplexD *b,
+		    const std::vector<const Lattice<vobj>*> &x,int m,
+		    int do_norm,
+		    decltype(innerProduct(vobj(),vobj())) *inner_tmp_v)
 {
   typedef decltype(z.View(AcceleratorRead)) View;
-
+  GRID_ASSERT(m>=1 && m<=B);
-  int m = x.size();
-  GRID_ASSERT((int)b.size()>=m);
   GridBase *grid = z.Grid();
   const uint64_t nsimd = grid->Nsimd();
   const uint64_t sites = grid->oSites();
 
-  hostVector<View>   h_x_v(std::max(m,1));
+  std::vector<View> h_v; h_v.reserve(m);
-  deviceVector<View> d_x_v(std::max(m,1));
+  ViewPack<View,B> pack;
-  hostVector<ComplexD>   h_b(std::max(m,1));
-  deviceVector<ComplexD> d_b(std::max(m,1));
   for(int j=0;j<m;j++){
     conformable(*x[j],z);
     GRID_ASSERT(x[j]!=&z);            // window must not alias the accumulator
-    h_x_v[j] = x[j]->View(AcceleratorRead);
+    h_v.push_back(x[j]->View(AcceleratorRead));
-    h_b[j]   = b[j];
+    pack.set(j,h_v[j]);
+    pack.b[j] = b[j];
   }
-  if ( m ) {
+  for(int j=m;j<B;j++){ pack.set(j,h_v[0]); pack.b[j] = ComplexD(0.0); }
-    acceleratorCopyToDevice(&h_x_v[0],&d_x_v[0],m*sizeof(View));
-    acceleratorCopyToDevice(&h_b[0],&d_b[0],m*sizeof(ComplexD));
-  }
-  View     *x_vp = &d_x_v[0];
-  ComplexD *b_p  = &d_b[0];
 
   autoView(z_v,z,AcceleratorWrite);
-  typedef decltype(innerProduct(z_v[0],z_v[0])) inner_t;
-  deviceVector<inner_t> inner_tmp(do_norm ? sites : 1);
-  inner_t *inner_tmp_v = &inner_tmp[0];
-
   accelerator_for(ss,sites,nsimd,{
       auto acc = coalescedRead(z_v[ss]);
-      for(int j=0;j<m;j++) acc = acc + b_p[j]*coalescedRead(x_vp[j][ss]);
+      for(int j=0;j<B;j++) if ( j<m ) acc = acc + pack.b[j]*coalescedRead(pack.v(j)[ss]);
       coalescedWrite(z_v[ss],acc);
       if ( do_norm ) coalescedWrite(inner_tmp_v[ss],innerProduct(acc,acc));
   });
-  for(int j=0;j<m;j++) h_x_v[j].ViewClose();
+  for(int j=0;j<m;j++) h_v[j].ViewClose();
+}
+
+// z = z + sum_j b[j] x[j];  if do_norm, returns global |z|^2 from the same
+// (last) pass.  Chunks of 16 for windows longer than 16.
+template<class vobj>
+RealD axpyMultiNormImpl(Lattice<vobj> &z,const std::vector<ComplexD> &b,
+			const std::vector<const Lattice<vobj>*> &x,int do_norm)
+{
+  typedef decltype(innerProduct(vobj(),vobj())) inner_t;
+  int m = x.size();
+  GRID_ASSERT((int)b.size()>=m);
+  GridBase *grid = z.Grid();
+  const uint64_t sites = grid->oSites();
+
+  deviceVector<inner_t> inner_tmp(do_norm ? sites : 1);
+  inner_t *inner_tmp_v = &inner_tmp[0];
+
+  if ( m==0 ) {
+    if ( do_norm ) return norm2(z);
+    return 0.0;
+  }
+  for(int j0=0;j0<m;j0+=16){
+    int mm   = std::min(16,m-j0);
+    int last = (j0+mm>=m);
+    std::vector<const Lattice<vobj>*> sub(x.begin()+j0,x.begin()+j0+mm);
+    int dn = do_norm && last;
+    if      ( mm<=2 ) axpyMultiChunk<2> (z,&b[j0],sub,mm,dn,inner_tmp_v);
+    else if ( mm<=4 ) axpyMultiChunk<4> (z,&b[j0],sub,mm,dn,inner_tmp_v);
+    else if ( mm<=8 ) axpyMultiChunk<8> (z,&b[j0],sub,mm,dn,inner_tmp_v);
+    else              axpyMultiChunk<16>(z,&b[j0],sub,mm,dn,inner_tmp_v);
+  }
 
   RealD nrm = 0.0;
   if ( do_norm ) {

systems/Frontier/slate_vs_schur2d.job

@@ -52,15 +52,24 @@ chmod +x ./select_gpu
 root=$HOME/ParallelIO/systems/Frontier
 source $root/sourceme-rocm7.2.sh
 
-export OMP_NUM_THREADS=7
+# ONE host thread per rank, for both codes.  Grid uses no host OpenMP on the
+# GPU build; SLATE with 7 threads issuing device-buffer MPI concurrently
+# deadlocks in getrf's tile broadcast under Cray MPICH (slate_debug.job V1
+# hung, V2 with OMP_NUM_THREADS=1 completed, 2026-08-25).  Each rank has a
+# whole GCD; the parallelism is on the device for both.
+export OMP_NUM_THREADS=1
 export MPICH_GPU_SUPPORT_ENABLED=1
 export MPICH_SMP_SINGLE_COPY_MODE=CMA
 export MPICH_OFI_NIC_POLICY=GPU
-export MPICH_MAX_THREAD_SAFETY=multiple   # SLATE calls MPI from OpenMP tasks
+export MPICH_MAX_THREAD_SAFETY=multiple   # SLATE still requests MPI_THREAD_MULTIPLE
 
 module load libfabric
 
 BIN=$root/tests/debug/Test_schur2d_vs_slate
+# Three legs per stage: Grid 2D Schur, SLATE getrf+getri, SLATE getrf+getrs(I).
+# getri is a host loop in SLATE (minutes at N=138240; measured once already):
+# S2D_SKIP_GETRI=1 drops it so S3 is ~3 min.  Unset to measure it again.
+export S2D_SKIP_GETRI=${S2D_SKIP_GETRI:-1}
 OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000"
 
 echo "========================================================="
@@ -80,7 +89,7 @@ S2D_N=4096 srun -N1 -n8 ./select_gpu $BIN --mpi 1.1.2.4 --grid 16.16.16.16 $OPTS
      2>&1 | tee s1.out
 S1RC=${PIPESTATUS[0]}
 S1CERT=$(grep -c "certificate" s1.out)
-echo "S1 exit code $S1RC, certificates printed $S1CERT (expect 2)"
+echo "S1 exit code $S1RC, certificates printed $S1CERT (expect 2, or 3 with getri)"
 
 ##############################################################################
 echo "========================================================="
@@ -94,7 +103,7 @@ echo "========================================================="
 echo "S3: THE comparison, 36 nodes, 288 ranks, N=138240 (nb=480, grid 16x18)"
 echo "========================================================="
 ##############################################################################
-if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -eq 2 ]
+if [ "$S1RC" -eq 0 ] && [ "$S1CERT" -ge 2 ]
 then
   S2D_N=138240 srun -N36 -n288 ./select_gpu $BIN --mpi 3.6.4.4 --grid 48.48.48.96 $OPTS1
   echo "S3 exit code $?"
@@ -105,6 +114,6 @@ fi
 echo "========================================================="
 echo "summary (both legs, all stages)"
 echo "========================================================="
-grep -h -E "Grid-vs-SLATE|GRID  :|SLATE :" slurm-$SLURM_JOB_ID.out 2>/dev/null
+grep -h -E "Grid-vs-SLATE|GRID  :|SLATE :|SLATE-getrs" slurm-$SLURM_JOB_ID.out 2>/dev/null
 grep -h "git commit hash" slurm-$SLURM_JOB_ID.out 2>/dev/null | sort -u
 echo "========================================================="

tests/debug/Test_schur2d_vs_slate.cc

@@ -55,7 +55,19 @@ Author: Peter Boyle <pboyle@bnl.gov>
 // grids, nb dividing N (48|720) and ragged (N=730, nb=50); both legs certify
 // max|A.Ainv-I| ~ 1e-15.
 //
+// Frontier (HIP, Target::Devices, Cray MPICH, 2026-08-25): run with
+// OMP_NUM_THREADS=1.  With 7 threads per rank SLATE's getrf deadlocks in its
+// tile broadcast (concurrent device-buffer MPI from OpenMP tasks); with one
+// thread it completes (8 GCDs, N=4096: certificate 2.7e-15).  One host
+// thread per rank is the like-for-like anyway: Grid's GPU build uses none.
+// Also: the site `slate` module (cpu env) must NOT be loaded -- its host-only
+// libblaspp shadows the ROCm one via LD_LIBRARY_PATH and throws
+// "device BLAS not available" from host_malloc_pinned.
+//
 //   S2D_N, S2D_NB as in Test_schur2d_scale (default nb = N/P).
+//   S2D_SKIP_GETRI=1 skips the getri leg (host loop; ~4 min at N=138240).
+//   S2D_NOWARM=1 skips the warm-up.
+// A third leg, getrf+getrs(I), is SLATE's device-resident inverse route.
 //////////////////////////////////////////////////////////////////////////////
 
 #include <Grid/Grid.h>
@@ -153,13 +165,22 @@ int main(int argc, char **argv)
   // a hang shows WHERE even through block-buffered stdout.
   auto Stage = [&](const char *s){ std::cout << GridLogMessage << "stage: " << s << std::endl << std::flush; };
   if ( !getenv("S2D_NOWARM") ) {
-    int64_t Nw = 8*P; int64_t nbw = 8;
+    // Fixed tiny size independent of P: the purpose is handle creation and
-    std::vector<int64_t> rs(P+1); for(int r=0;r<=P;r++) rs[r]=8*r;
+    // kernel loading, not work.  (8*P at P=288 was N=2304 -> a 122 s SLATE
-    std::vector<ComplexD> hw(8*Nw); for(int64_t jj=0;jj<Nw;jj++) for(int64_t i=0;i<8;i++) hw[i+jj*8]=Fill(8*me+i,jj);
+    // warm-up dominated by 288-way tile broadcasts.)  Ranks beyond the first
+    // Nw/nbw own nothing in the rows layout, which RowsToCyclic handles.
+    int64_t Nw = 64; int64_t nbw = 8;
+    // same partition rule as the main leg: Nw/P rows each, remainder to the
+    // first ranks; at P > Nw most ranks contribute zero rows.
+    std::vector<int64_t> rs(P+1); rs[0]=0;
+    for(int r=0;r<P;r++) rs[r+1] = rs[r] + Nw/P + ( r < (int)(Nw%P) ? 1 : 0 );
+    int64_t rw = rs[me+1]-rs[me];
+    std::vector<ComplexD> hw(std::max<int64_t>(rw,1)*Nw);
+    for(int64_t jj=0;jj<Nw;jj++) for(int64_t i=0;i<rw;i++) hw[i+jj*rw]=Fill(rs[me]+i,jj);
     deviceVector<ComplexD> dw(hw.size()); acceleratorCopyToDevice(&hw[0],&dw[0],hw.size()*sizeof(ComplexD));
     BlockCyclicMatrix W(grid,Nw,nbw,Pr,Pc);
     Stage("warm-up Grid RowsToCyclic");
-    BlockCyclicRedistribute::RowsToCyclic(grid,rs,&dw[0],8,W);
+    BlockCyclicRedistribute::RowsToCyclic(grid,rs,&dw[0],rw,W);
     Stage("warm-up Grid Invert");
     BlockCyclicSchurInverse RSIw; RSIw.Invert(W);
 #ifdef HAVE_SLATE
@@ -203,7 +224,7 @@ int main(int argc, char **argv)
       acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],A.data.size()*sizeof(ComplexD));
     double t2=usecond();
     Stage("Grid Invert");
-    RSI2.Invert(A);
+    { GRID_TRACE("GridInvert"); RSI2.Invert(A); }
     double t3=usecond();
     BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows);
     double t4=usecond();
@@ -218,7 +239,7 @@ int main(int argc, char **argv)
   // LEG 2: SLATE, every layout step timed and charged.
   ////////////////////////////////////////////////////////////////////////
 #ifdef HAVE_SLATE
-  {
+  if ( !getenv("S2D_SKIP_GETRI") ) {   // S2D_SKIP_GETRI=1: getri is a host loop, minutes at N=138240
     typedef std::complex<double> scalar_t;
     acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
     BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc);
@@ -262,10 +283,10 @@ int main(int argc, char **argv)
     slate::Pivots pivots;
     Stage("SLATE getrf");
     double t4=usecond();
-    slate::getrf(S, pivots, opts);          // LU, partial pivoting
+    { GRID_TRACE("SLATE_getrf"); slate::getrf(S, pivots, opts); }   // LU, partial pivoting
     double t5=usecond();
     Stage("SLATE getri");
-    slate::getri(S, pivots, opts);          // in-place inverse from the factor
+    { GRID_TRACE("SLATE_getri"); slate::getri(S, pivots, opts); }   // in-place inverse from the factor
     double t6=usecond();
 
     // back onto the device, in our layout (getri applies the pivots itself)
@@ -281,6 +302,72 @@ int main(int argc, char **argv)
               << "  TOTAL " << ((t1-t0)+(t8-t2))/1e6 << " s"
               << "   certificate " << cert << std::endl;
   }
+
+  ////////////////////////////////////////////////////////////////////////
+  // LEG 3: SLATE getrf + getrs(I) -- SLATE's device-resident route to an
+  // explicit inverse.  getri's L-side loop is hard-coded Target::HostTask
+  // (src/getri.cc: copy/gemmA/trsm/permuteRows all <HostTask>), so under
+  // Target::Devices it runs as a 288-step host loop on one thread per rank.
+  // getrs is two trsm's that honour the target (plus host row permutes).
+  // The identity RHS is built on the host in the same ScaLAPACK layout and
+  // its construction is on SLATE's clock.
+  ////////////////////////////////////////////////////////////////////////
+  {
+    typedef std::complex<double> scalar_t;
+    acceleratorCopyToDevice(&h[0], &rows1d[0], h.size()*sizeof(ComplexD));
+    BlockCyclicMatrix A(grid,N,nb,Pr,Pc), A0(grid,N,nb,Pr,Pc);
+    double t0=usecond();
+    BlockCyclicRedistribute::RowsToCyclic(grid,rowStart,&rows1d[0],myrows,A);
+    double t1=usecond();
+    if ( A.data.size() )
+      acceleratorCopyDeviceToDevice((void *)&A.data[0],(void *)&A0.data[0],A.data.size()*sizeof(ComplexD));
+    BlockCyclicLayout &L = A.layout;
+    uint64_t nloc = (uint64_t)L.mloc*L.nloc;
+    std::vector<scalar_t> hA(nloc ? nloc : 1), hB(nloc ? nloc : 1);
+    double t2=usecond();
+    if ( nloc ) acceleratorCopyFromDevice(&A.data[0], (void *)&hA[0], nloc*sizeof(ComplexD));
+    for(int64_t lj=0;lj<L.nloc;lj++){
+      int64_t gj = BlockCyclicLayout::LocalToGlobal(lj, nb, L.pcol, Pc);
+      for(int64_t li=0;li<L.mloc;li++){
+        int64_t gi = BlockCyclicLayout::LocalToGlobal(li, nb, L.prow, Pr);
+        hB[li+lj*L.mloc] = (gi==gj) ? scalar_t(1.0,0.0) : scalar_t(0.0,0.0);
+      }
+    }
+    double t3=usecond();
+    auto S = slate::Matrix<scalar_t>::fromScaLAPACK(N, N, &hA[0], (int64_t)std::max<int64_t>(L.mloc,1),
+                                                    nb, nb, slate::GridOrder::Row, Pr, Pc, grid->communicator);
+    auto B = slate::Matrix<scalar_t>::fromScaLAPACK(N, N, &hB[0], (int64_t)std::max<int64_t>(L.mloc,1),
+                                                    nb, nb, slate::GridOrder::Row, Pr, Pc, grid->communicator);
+#if defined(GRID_HIP) || defined(GRID_CUDA) || defined(GRID_SYCL)
+    slate::Target target = slate::Target::Devices;
+#else
+    slate::Target target = slate::Target::HostTask;
+#endif
+    slate::Options opts = {
+      { slate::Option::Target,        target },
+      { slate::Option::Lookahead,     1 },
+      { slate::Option::InnerBlocking, 16 },
+    };
+    slate::Pivots pivots;
+    Stage("SLATE getrf (getrs leg)");
+    double t4=usecond();
+    { GRID_TRACE("SLATE_getrf"); slate::getrf(S, pivots, opts); }
+    double t5=usecond();
+    Stage("SLATE getrs(I)");
+    { GRID_TRACE("SLATE_getrs"); slate::getrs(S, pivots, B, opts); }   // B <- A^{-1} I
+    double t6=usecond();
+    if ( nloc ) acceleratorCopyToDevice((void *)&hB[0], &A.data[0], nloc*sizeof(ComplexD));
+    double t7=usecond();
+    BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows);
+    double t8=usecond();
+    double cert = Certify(grid,A0,A,nb,Pr,Pc);
+    std::cout << GridLogMessage << "SLATE-getrs : redist->2D " << (t1-t0)/1e6
+              << "  D2H+I " << (t3-t2)/1e6 << "  wrap " << (t4-t3)/1e6
+              << "  getrf " << (t5-t4)/1e6 << "  getrs " << (t6-t5)/1e6
+              << "  H2D " << (t7-t6)/1e6 << "  redist->rows " << (t8-t7)/1e6
+              << "  TOTAL " << ((t1-t0)+(t8-t2))/1e6 << " s"
+              << "   certificate " << cert << std::endl;
+  }
 #else
   std::cout << GridLogMessage << "SLATE : leg not built (compile with -DHAVE_SLATE and link -lslate -lblaspp -llapackpp)" << std::endl;
 #endif