SSIDS: front-size- and thread-adaptive CPU block size

driver/spral_ssids.F90

@@ -36,9 +36,9 @@ program run_prob
 
   real :: smanal, smfact, smaflop, smafact
 
-  integer, parameter :: nfact = 1
-  ! integer, parameter :: nfact = 50
-  ! integer, parameter :: nfact = 100
+  integer :: nfact = 1     ! override with env var SPRAL_NFACT for timing loops
+  character(len=32) :: nfact_env
+  integer :: nfact_stat
 
   integer, parameter :: nslv = 1
   ! integer, parameter :: nslv = 10
@@ -59,6 +59,9 @@ program run_prob
        flat_topology, ngpus)
   if (nrhs .lt. 1) stop
 
+  call get_environment_variable("SPRAL_NFACT", nfact_env, status=nfact_stat)
+  if (nfact_stat .eq. 0) read(nfact_env, *) nfact
+
   ! Read in a matrix
   write (*, "(3a)") "Reading '", filename, "'..."
   if (force_psdef) rb_options%values = -3 ! Force diagonal dominance

meson.build

@@ -111,6 +111,43 @@ if build_gpu and nvcc_available and libcuda.found()
      add_global_arguments('-DHAVE_NVCC', language : 'cpp')
 endif
 
+# SSIDS front-size-adaptive block size tuning constants (see
+# src/ssids/cpu/kernels/block_size.hxx). Set at configure time, e.g.
+#   meson configure -Dssids_block_div=0 -Dssids_block_min=128 -Dssids_block_max=512
+ssids_block_max = get_option('ssids_block_max')
+
+# The block-size floor (ssids_block_min): a positive option value is used as-is;
+# 0 (the default) auto-selects it with a configure-time BLAS efficiency probe
+# (meson/blas_block_probe.cxx), which reports the smallest square tile at which
+# this machine's BLAS runs near peak. The probe is skipped when cross-compiling
+# and any failure falls back to a safe default, so it never blocks configure.
+ssids_block_min_opt = get_option('ssids_block_min')
+ssids_block_min_fallback = 128
+if ssids_block_min_opt > 0
+  ssids_block_min = ssids_block_min_opt
+elif meson.is_cross_build()
+  ssids_block_min = ssids_block_min_fallback
+  message('SSIDS block-size floor: cross build, using fallback @0@'.format(ssids_block_min))
+else
+  block_probe = cxx.run(files('meson/blas_block_probe.cxx'),
+                        dependencies : [libblas, lm],
+                        name : 'ssids blas block-size floor probe')
+  if block_probe.compiled() and block_probe.returncode() == 0
+    ssids_block_min = block_probe.stdout().strip().to_int()
+    message('SSIDS block-size floor: BLAS probe selected @0@'.format(ssids_block_min))
+  else
+    ssids_block_min = ssids_block_min_fallback
+    message('SSIDS block-size floor: probe did not run, using fallback @0@'.format(ssids_block_min))
+  endif
+endif
+if ssids_block_max < ssids_block_min
+  error('ssids_block_max (@0@) must be >= ssids_block_min (@1@)'.format(ssids_block_max, ssids_block_min))
+endif
+add_global_arguments('-DSPRAL_SSIDS_BLOCK_DIV=@0@'.format(get_option('ssids_block_div')), language : 'cpp')
+add_global_arguments('-DSPRAL_SSIDS_BLOCK_TILES_PER_THREAD=@0@'.format(get_option('ssids_block_tiles_per_thread')), language : 'cpp')
+add_global_arguments('-DSPRAL_SSIDS_BLOCK_MIN=@0@'.format(ssids_block_min), language : 'cpp')
+add_global_arguments('-DSPRAL_SSIDS_BLOCK_MAX=@0@'.format(ssids_block_max), language : 'cpp')
+
 # OpenMP
 if build_openmp
   if fc.get_id() == 'nvidia_hpc'

meson/blas_block_probe.cxx

@@ -0,0 +1,111 @@
+/* blas_block_probe.cxx
+ *
+ * Configure-time probe that chooses the SSIDS front-size-adaptive block-size
+ * floor (SPRAL_SSIDS_BLOCK_MIN) from this machine's BLAS. It measures the
+ * single-threaded double-precision GEMM rate for a range of small square tile
+ * sizes and reports the smallest tile at which the BLAS runs near peak -- i.e.
+ * the smallest tile for which the tiled dense factorization is BLAS-efficient.
+ * Below that size the kernels run below peak and per-task/per-call overhead
+ * dominates, so it is a sensible lower clamp for the adaptive ramp.
+ *
+ * The program prints a single integer (a multiple of 32, clamped to a sane
+ * range) on stdout; the Meson build captures it via compiler.run(). It is
+ * deliberately self-contained and quick (well under a second). If anything is
+ * unexpected it prints a conservative fallback and still exits 0, so the build
+ * is never blocked by the probe.
+ *
+ * Note: because it times the machine at configure time, the chosen value can
+ * vary slightly between configurations run under different load. Builds that
+ * need a fixed, reproducible floor should set the Meson option
+ * -Dssids_block_min=<value> to a positive number, which bypasses this probe.
+ */
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include <chrono>
+#include <algorithm>
+
+// Reference-BLAS Fortran symbol (OpenBLAS, Netlib, Accelerate, MKL all export
+// this). Column-major, no transpose.
+extern "C" void dgemm_(const char* transa, const char* transb,
+                       const int* m, const int* n, const int* k,
+                       const double* alpha, const double* a, const int* lda,
+                       const double* b, const int* ldb,
+                       const double* beta, double* c, const int* ldc);
+
+namespace {
+
+// GEMM rate (flop/s) for a square b-by-b-by-b product, best of a few trials.
+double gemm_rate(int b) {
+   std::vector<double> A(static_cast<size_t>(b) * b, 1.0);
+   std::vector<double> B(static_cast<size_t>(b) * b, 1.0);
+   std::vector<double> C(static_cast<size_t>(b) * b, 0.0);
+   const char N = 'N';
+   const double one = 1.0, zero = 0.0;
+   const double work = 2.0 * b * b * b; // flop per gemm
+
+   // Enough repetitions that each timed burst lasts a few milliseconds.
+   int reps = static_cast<int>(3.0e8 / work);
+   if (reps < 5) reps = 5;
+
+   // Warm up (first call may allocate BLAS buffers / spin up threads).
+   dgemm_(&N, &N, &b, &b, &b, &one, A.data(), &b, B.data(), &b, &zero, C.data(), &b);
+
+   double best = 0.0;
+   for (int trial = 0; trial < 3; ++trial) {
+      auto t0 = std::chrono::steady_clock::now();
+      for (int r = 0; r < reps; ++r)
+         dgemm_(&N, &N, &b, &b, &b, &one, A.data(), &b, B.data(), &b,
+                &zero, C.data(), &b);
+      auto t1 = std::chrono::steady_clock::now();
+      double secs = std::chrono::duration<double>(t1 - t0).count();
+      if (secs > 0.0) best = std::max(best, work * reps / secs);
+   }
+   return best;
+}
+
+} // anonymous namespace
+
+int main() {
+   // Measure single-threaded BLAS: the tiled kernels each run on one thread
+   // inside the task scheduler, so single-thread efficiency is what sets the
+   // useful floor. Set before the first BLAS call so the runtime picks it up.
+   setenv("OPENBLAS_NUM_THREADS", "1", 1);
+   setenv("OMP_NUM_THREADS", "1", 1);
+   setenv("MKL_NUM_THREADS", "1", 1);
+   setenv("VECLIB_MAXIMUM_THREADS", "1", 1); // Apple Accelerate
+
+   const int gran = 32;   // report a multiple of the LDLT inner block size
+   const int lo   = 64;   // never floor below this
+   const int hi   = 256;  // never floor above this
+   const double frac = 0.90; // "near peak" threshold
+
+   // Sample square tile sizes from small to a size comfortably past the knee.
+   std::vector<int> sizes;
+   for (int b = 32; b <= 384; b += 32) sizes.push_back(b);
+
+   double peak = 0.0;
+   std::vector<double> rate(sizes.size(), 0.0);
+   for (size_t i = 0; i < sizes.size(); ++i) {
+      rate[i] = gemm_rate(sizes[i]);
+      peak = std::max(peak, rate[i]);
+   }
+
+   // Smallest tile reaching frac of peak.
+   int chosen = hi;
+   if (peak > 0.0) {
+      for (size_t i = 0; i < sizes.size(); ++i) {
+         if (rate[i] >= frac * peak) { chosen = sizes[i]; break; }
+      }
+   } else {
+      chosen = 128; // BLAS produced no timing; fall back conservatively
+   }
+
+   // Round to a multiple of the granularity and clamp to [lo, hi].
+   chosen = ((chosen + gran / 2) / gran) * gran;
+   if (chosen < lo) chosen = lo;
+   if (chosen > hi) chosen = hi;
+
+   std::printf("%d\n", chosen);
+   return 0;
+}

meson_options.txt

@@ -88,3 +88,27 @@ option('openmp',
        type : 'boolean',
        value : true,
        description : 'option to compile SPRAL with OpenMP')
+
+option('ssids_block_div',
+       type : 'integer',
+       min : 0,
+       value : 0,
+       description : 'SSIDS front-size-adaptive block size: fixed ramp divisor (block ~= front_rows / div); 0 selects the thread-scaled divisor ssids_block_tiles_per_thread * nthreads')
+
+option('ssids_block_tiles_per_thread',
+       type : 'integer',
+       min : 1,
+       value : 2,
+       description : 'SSIDS front-size-adaptive block size: target block-rows per thread when ssids_block_div is 0 (divisor = this * nthreads)')
+
+option('ssids_block_min',
+       type : 'integer',
+       min : 0,
+       value : 0,
+       description : 'SSIDS front-size-adaptive block size: minimum (floor) block size; 0 auto-selects it at configure time with a BLAS efficiency probe (meson/blas_block_probe.cxx)')
+
+option('ssids_block_max',
+       type : 'integer',
+       min : 32,
+       value : 512,
+       description : 'SSIDS front-size-adaptive block size: maximum (ceiling) block size')

src/ssids/cpu/NumericSubtree.hxx

@@ -160,7 +160,7 @@ public:
             auto* this_lcol = &nodes_[ni]; // for depend
             auto* parent_lcol = nodes_.data() + symb_[ni].parent; // for depend
             #pragma omp task default(none) \
-               firstprivate(ni) \
+               firstprivate(ni, num_threads) \
                shared(aval, abort, child_contrib, options, scaling, \
                       thread_stats, work) \
                depend(inout: this_lcol[0:1]) \
@@ -192,7 +192,7 @@ public:
                   factor_node<posdef>
                      (ni, symb_[ni], nodes_[ni], options,
                       thread_stats[this_thread], work,
-                      pool_alloc_);
+                      pool_alloc_, num_threads);
                   if(thread_stats[this_thread].flag<Flag::SUCCESS) {
 #ifdef _OPENMP
                      #pragma omp atomic write

src/ssids/cpu/factor.hxx

@@ -42,7 +42,8 @@ void factor_node_indef(
       struct cpu_factor_options const& options,
       ThreadStats& stats,
       std::vector<Workspace>& work,
-      PoolAlloc& pool_alloc
+      PoolAlloc& pool_alloc,
+      int nthreads
       ) {
    /* Extract useful information about node */
    int m = snode.nrow + node.ndelay_in;
@@ -59,7 +60,7 @@ void factor_node_indef(
       // Use an APP based pivot method
       node.nelim = ldlt_app_factor(
             m, n, perm, lcol, ldl, d, 0.0, contrib, m-n, options, work,
-            pool_alloc
+            pool_alloc, nthreads
             );
       if(node.nelim < 0) {
          stats.flag = static_cast<Flag>(node.nelim);
@@ -141,7 +142,8 @@ void factor_node_posdef(
       SymbolicNode const& snode,
       NumericNode<T, PoolAlloc> &node,
       struct cpu_factor_options const& options,
-      ThreadStats& stats
+      ThreadStats& stats,
+      int nthreads = 1 // serial by default (e.g. small-leaf subtree nodes)
       ) {
    /* Extract useful information about node */
    int m = snode.nrow;
@@ -153,7 +155,8 @@ void factor_node_posdef(
    /* Perform factorization */
    int flag;
    cholesky_factor(
-         m, n, lcol, ldl, beta, contrib, m-n, options.cpu_block_size, &flag
+         m, n, lcol, ldl, beta, contrib, m-n, options.cpu_block_size, &flag,
+         nthreads
          );
    if(flag!=-1) {
       node.nelim = flag+1;
@@ -178,10 +181,11 @@ void factor_node(
       struct cpu_factor_options const& options,
       ThreadStats& stats,
       std::vector<Workspace>& work,
-      PoolAlloc& pool_alloc
+      PoolAlloc& pool_alloc,
+      int nthreads
       ) {
-   if(posdef) factor_node_posdef(0.0, snode, node, options, stats);
-   else       factor_node_indef(ni, snode, node, options, stats, work, pool_alloc);
+   if(posdef) factor_node_posdef(0.0, snode, node, options, stats, nthreads);
+   else       factor_node_indef(ni, snode, node, options, stats, work, pool_alloc, nthreads);
 }
 
 }}} /* end of namespace spral::ssids::cpu */

src/ssids/cpu/kernels/block_size.hxx

@@ -0,0 +1,88 @@
+/** \file
+ *  \copyright 2026 The Science and Technology Facilities Council (STFC)
+ *  \licence   BSD licence, see LICENCE file for details
+ *
+ *  Front-size-adaptive outer block size for the dense node kernels.
+ */
+#pragma once
+
+#include <algorithm>
+#include <cstdint>
+
+/* Tuning constants for the front-size-adaptive block size. These are fixed at
+ * build time via the Meson options ssids_block_div / ssids_block_tiles_per_thread
+ * / ssids_block_min / ssids_block_max (which the build system turns into the
+ * macros below). SPRAL_SSIDS_BLOCK_MIN in particular is, by default, determined
+ * by a configure-time BLAS efficiency probe (see meson/blas_block_probe.cxx).
+ * The fallback defaults here keep the header usable in a standalone compile. */
+#ifndef SPRAL_SSIDS_BLOCK_DIV
+#define SPRAL_SSIDS_BLOCK_DIV 0            /* 0 => thread-scaled (see below) */
+#endif
+#ifndef SPRAL_SSIDS_BLOCK_TILES_PER_THREAD
+#define SPRAL_SSIDS_BLOCK_TILES_PER_THREAD 2
+#endif
+#ifndef SPRAL_SSIDS_BLOCK_MIN
+#define SPRAL_SSIDS_BLOCK_MIN 128
+#endif
+#ifndef SPRAL_SSIDS_BLOCK_MAX
+#define SPRAL_SSIDS_BLOCK_MAX 512
+#endif
+
+namespace spral { namespace ssids { namespace cpu {
+
+/** Inner block size of the dense LDLT kernel: the granularity of its recursive
+ *  diagonal-block factorization, and the unit to which adaptive_block_size()
+ *  rounds. Defined here so the LDLT (ldlt_app.cxx) and Cholesky (cholesky.cxx)
+ *  kernels share a single definition rather than each hard-coding 32. */
+static const int INNER_BLOCK_SIZE = 32;
+
+/** Front-size-adaptive outer block size.
+ *
+ *  The optimal block size for the dense factorization of a supernode is not a
+ *  global constant: it tracks both the front's size and how many threads are
+ *  working it.
+ *
+ *  - Absolute clamps [MIN,MAX] bound single-block efficiency (BLAS shape,
+ *    cache footprint, per-task overhead) and are independent of thread count.
+ *    MIN is, by default, chosen at build time by a BLAS efficiency probe: the
+ *    smallest square tile at which this machine's BLAS runs near peak.
+ *  - The ramp rate governs how finely a front is tiled, i.e. how much parallel
+ *    work it exposes -- so it scales with the number of threads. We aim for
+ *    roughly TILES_PER_THREAD block-rows per thread, giving a divisor of
+ *    TILES_PER_THREAD * nthreads. Small fronts and/or many threads therefore
+ *    pull the block size down (more tiles to keep everyone busy) until the MIN
+ *    floor; very large fronts and/or few threads push it up to the MAX ceiling.
+ *
+ *  \param m front row count
+ *  \param granularity block size is rounded to a multiple of this; both kernels
+ *         pass INNER_BLOCK_SIZE (the LDLT inner block size)
+ *  \param nthreads number of threads collaborating on this front, i.e. the size
+ *         of the OpenMP team assigned to this subtree by SSIDS' topology model
+ *         (captured once at subtree entry). Serial callers should pass 1, for
+ *         which the rule selects large tiles. Values \f$\le 0\f$ are treated as
+ *         1 defensively.
+ *
+ *  A build may instead pin a fixed divisor (thread-independent) by setting the
+ *  Meson option ssids_block_div to a positive value; 0 (the default) selects
+ *  the thread-scaled divisor above.
+ */
+inline int adaptive_block_size(int m, int granularity, int nthreads) {
+   int const min_blk = std::max(granularity, SPRAL_SSIDS_BLOCK_MIN);
+   int const max_blk = std::max(min_blk, SPRAL_SSIDS_BLOCK_MAX);
+
+   int const nt = std::max(1, nthreads);
+
+   // Ramp divisor: a positive build-time override, else thread-scaled.
+   int const fixed_div = SPRAL_SSIDS_BLOCK_DIV;
+   int64_t const div = (fixed_div > 0)
+      ? static_cast<int64_t>(fixed_div)
+      : static_cast<int64_t>(SPRAL_SSIDS_BLOCK_TILES_PER_THREAD) * nt;
+
+   int b = static_cast<int>(static_cast<int64_t>(m) / div);
+   b = ((b + granularity/2) / granularity) * granularity; // round to nearest
+   if(b < min_blk) b = min_blk;
+   if(b > max_blk) b = max_blk;
+   return b;
+}
+
+}}} /* namespaces spral::ssids::cpu */

src/ssids/cpu/kernels/cholesky.cxx

@@ -10,6 +10,7 @@
 #include <cstdio> // FIXME: remove as only used for debug
 
 #include "ssids/profile.hxx"
+#include "ssids/cpu/kernels/block_size.hxx"
 #include "ssids/cpu/kernels/wrappers.hxx"
 
 namespace spral { namespace ssids { namespace cpu {
@@ -30,7 +31,11 @@ namespace spral { namespace ssids { namespace cpu {
  * \param info is initialized to -1, and will be changed to the index of any
  *    column where a non-zero column is encountered.
  */
-void cholesky_factor(int m, int n, double* a, int lda, double beta, double* upd, int ldupd, int blksz, int *info) {
+void cholesky_factor(int m, int n, double* a, int lda, double beta, double* upd, int ldupd, int blksz, int *info, int nthreads) {
+   // A non-positive block size opts into the front-size-adaptive rule; a
+   // positive value keeps the historic fixed-block behaviour. nthreads is the
+   // size of the OpenMP team assigned to this subtree.
+   if(blksz <= 0) blksz = adaptive_block_size(m, INNER_BLOCK_SIZE, nthreads);
    if(n < blksz) {
       // Adjust so blocks have blksz**2 entries
       blksz = int((int64_t(blksz)*blksz) / n);