diff --git a/driver/spral_ssids.F90 b/driver/spral_ssids.F90 index 2e82f927..a244d172 100644 --- a/driver/spral_ssids.F90 +++ b/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 diff --git a/meson.build b/meson.build index 86a024af..c8a150d9 100644 --- a/meson.build +++ b/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' diff --git a/meson/blas_block_probe.cxx b/meson/blas_block_probe.cxx new file mode 100644 index 00000000..1e040a39 --- /dev/null +++ b/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= to a positive number, which bypasses this probe. + */ +#include +#include +#include +#include +#include + +// 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 A(static_cast(b) * b, 1.0); + std::vector B(static_cast(b) * b, 1.0); + std::vector C(static_cast(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(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(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 sizes; + for (int b = 32; b <= 384; b += 32) sizes.push_back(b); + + double peak = 0.0; + std::vector 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; +} diff --git a/meson_options.txt b/meson_options.txt index 2edd0bb3..c29afc19 100644 --- a/meson_options.txt +++ b/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') diff --git a/src/ssids/cpu/NumericSubtree.hxx b/src/ssids/cpu/NumericSubtree.hxx index 980765f0..45505fbe 100644 --- a/src/ssids/cpu/NumericSubtree.hxx +++ b/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 (ni, symb_[ni], nodes_[ni], options, thread_stats[this_thread], work, - pool_alloc_); + pool_alloc_, num_threads); if(thread_stats[this_thread].flag& 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(node.nelim); @@ -141,7 +142,8 @@ void factor_node_posdef( SymbolicNode const& snode, NumericNode &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& 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 */ diff --git a/src/ssids/cpu/kernels/block_size.hxx b/src/ssids/cpu/kernels/block_size.hxx new file mode 100644 index 00000000..9691127d --- /dev/null +++ b/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 +#include + +/* 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(fixed_div) + : static_cast(SPRAL_SSIDS_BLOCK_TILES_PER_THREAD) * nt; + + int b = static_cast(static_cast(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 */ diff --git a/src/ssids/cpu/kernels/cholesky.cxx b/src/ssids/cpu/kernels/cholesky.cxx index 970a482d..befbbc36 100644 --- a/src/ssids/cpu/kernels/cholesky.cxx +++ b/src/ssids/cpu/kernels/cholesky.cxx @@ -10,6 +10,7 @@ #include // 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); diff --git a/src/ssids/cpu/kernels/cholesky.hxx b/src/ssids/cpu/kernels/cholesky.hxx index c1d71322..d911dd01 100644 --- a/src/ssids/cpu/kernels/cholesky.hxx +++ b/src/ssids/cpu/kernels/cholesky.hxx @@ -5,7 +5,10 @@ */ namespace spral { namespace ssids { namespace cpu { -void cholesky_factor(int m, int n, double* a, int lda, double beta, double* upd, int ldupd, int blksz, int *info); +// nthreads is the size of the OpenMP team assigned to this subtree; it is only +// consulted when blksz<=0 (front-size-adaptive block size). It defaults to 1 +// (serial) for callers that always pass an explicit positive blksz. +void cholesky_factor(int m, int n, double* a, int lda, double beta, double* upd, int ldupd, int blksz, int *info, int nthreads = 1); void cholesky_solve_fwd(int m, int n, double const* a, int lda, int nrhs, double* x, int ldx); void cholesky_solve_bwd(int m, int n, double const* a, int lda, int nrhs, double* x, int ldx); diff --git a/src/ssids/cpu/kernels/ldlt_app.cxx b/src/ssids/cpu/kernels/ldlt_app.cxx index aa419436..325fa0e8 100644 --- a/src/ssids/cpu/kernels/ldlt_app.cxx +++ b/src/ssids/cpu/kernels/ldlt_app.cxx @@ -29,6 +29,7 @@ #include "ssids/cpu/cpu_iface.hxx" #include "ssids/cpu/Workspace.hxx" #include "ssids/cpu/kernels/block_ldlt.hxx" +#include "ssids/cpu/kernels/block_size.hxx" #include "ssids/cpu/kernels/calc_ld.hxx" #include "ssids/cpu/kernels/ldlt_tpp.hxx" #include "ssids/cpu/kernels/common.hxx" @@ -38,7 +39,8 @@ namespace spral { namespace ssids { namespace cpu { namespace ldlt_app_internal { -static const int INNER_BLOCK_SIZE = 32; +// INNER_BLOCK_SIZE is defined once in block_size.hxx (included above) and shared +// with the Cholesky kernel; it resolves here via the enclosing namespace. /** \return number of blocks for given n */ inline int calc_nblk(int n, int block_size) { @@ -2503,15 +2505,15 @@ size_t ldlt_app_factor_mem_required(int m, int n, int block_size) { } template -int ldlt_app_factor(int m, int n, int* perm, T* a, int lda, T* d, T beta, T* upd, int ldupd, struct cpu_factor_options const& options, std::vector& work, Allocator const& alloc) { - // If we've got a tall and narrow node, adjust block size so each block - // has roughly blksz**2 entries - // FIXME: Decide if this reshape is actually useful, given it will generate - // a lot more update tasks instead? - int outer_block_size = options.cpu_block_size; - /*if(n < outer_block_size) { - outer_block_size = int((int64_t(outer_block_size)*outer_block_size) / n); - }*/ +int ldlt_app_factor(int m, int n, int* perm, T* a, int lda, T* d, T beta, T* upd, int ldupd, struct cpu_factor_options const& options, std::vector& work, Allocator const& alloc, int nthreads) { + // Block size selection. A non-positive cpu_block_size opts into the + // front-size-adaptive rule; a positive value is used verbatim (the historic + // fixed-block behaviour). nthreads is the size of the OpenMP team assigned to + // this subtree (see NumericSubtree), so the ramp adapts to the parallelism + // actually available to this front. + int outer_block_size = (options.cpu_block_size > 0) + ? options.cpu_block_size + : adaptive_block_size(m, INNER_BLOCK_SIZE, nthreads); #ifdef PROFILE Profile::setState("TA_MISC1"); @@ -2532,7 +2534,7 @@ int ldlt_app_factor(int m, int n, int* perm, T* a, int lda, T* d, T beta, T* upd outer_block_size, beta, upd, ldupd, work, alloc ); } -template int ldlt_app_factor>>(int, int, int*, double*, int, double*, double, double*, int, struct cpu_factor_options const&, std::vector&, BuddyAllocator> const& alloc); +template int ldlt_app_factor>>(int, int, int*, double*, int, double*, double, double*, int, struct cpu_factor_options const&, std::vector&, BuddyAllocator> const& alloc, int nthreads); template void ldlt_app_solve_fwd(int m, int n, T const* l, int ldl, int nrhs, T* x, int ldx) { diff --git a/src/ssids/cpu/kernels/ldlt_app.hxx b/src/ssids/cpu/kernels/ldlt_app.hxx index 3e0ba66b..95bb4940 100644 --- a/src/ssids/cpu/kernels/ldlt_app.hxx +++ b/src/ssids/cpu/kernels/ldlt_app.hxx @@ -12,7 +12,7 @@ namespace spral { namespace ssids { namespace cpu { template -int ldlt_app_factor(int m, int n, int *perm, T *a, int lda, T *d, T beta, T* upd, int ldupd, struct cpu_factor_options const& options, std::vector& work, Allocator const& alloc); +int ldlt_app_factor(int m, int n, int *perm, T *a, int lda, T *d, T beta, T* upd, int ldupd, struct cpu_factor_options const& options, std::vector& work, Allocator const& alloc, int nthreads); template void ldlt_app_solve_fwd(int m, int n, T const* l, int ldl, int nrhs, T* x, int ldx);