GitHub

Fused pipeline

FusedSolver runs the fused factorization and partitioned-inverse solve of SparseDirectSolver.Fused (src/fused.jl) behind the verbs of DirectSolver: the same execute! phases, setparam!, getparam and update!, for one right-hand side, on a GPU backend, for "SPD" and no-pivot "S" real matrices. It is the inner solver of the MadNLP extension's SDSSolver (MadNLP); Performance says when to use it.

using SparseDirectSolver, CUDA, CUDA.CUSPARSE

Ad = CuSparseMatrixCSR(tril(A))       # A: SPD, on the host
F = FusedSolver(Ad, "SPD", 'L')
analyze!(F)                           # ordering, symbolic analysis, the fused plans
factorize!(F)                         # the fused factorization
bd = CuVector(b)
x = F \ bd                            # partitioned-inverse solve
getparam(F, "info")                   # 0 on success

Ad.nzVal .*= 2                        # new values in the same buffer (the IPM loop)
refactorize!(F)
solve!(F, x, bd)
SparseDirectSolver.Fused.FusedSolver — Type
FusedSolver(A, structure, view; tuning = nothing, index = 'O')
FusedSolver(rowptr, colval, nzval, structure, view; tuning = nothing, index = 'O')

The fused pipeline of SparseDirectSolver.Fused (dependency-counter factorization mega-kernel and partitioned-inverse solve, the SDSSolver path of the MadNLP extension) as an AbstractDirectSolver, driven by the same verbs as DirectSolver: execute!, analyze!, factorize!, refactorize!, solve!, setparam!, getparam, update!, ldiv! and \ for one right-hand side. A is anything DirectSolver accepts with these arguments (a CSR, CuSparseMatrixCSR/CuSparseMatrixCSC, ROCSparseMatrixCSR, …) on a GPU backend; the values alias the caller's nzval, so writing new values into it and calling refactorize! is the interior-point loop (update! with the same pattern is forwarded).

structure is "SPD" (fused Cholesky) or "S" (fused LDLᵀ without pivoting: pivot_type = 'N', pivot_pairs = "none"). tuning is an SDSTuning, nothing meaning default_tuning(backend), resolved by resolve_tuning; the constructor applies it to the inner DirectSolver (amalgamation, regime-C width and rows, subtree parallelism and budgets, reordering_alg = "algo4", nd_nseps, nd_seed). The inner solver is the field inner.

Phases (execute!(phase, F, X, B; asynchronous = true)):

  • "reordering", "symbolic_factorization", "analysis": the inner solver's phase; the symbolic part then factorizes and solves once with the inner solver (stock kernels, it sizes the solve workspace) and builds the fused factorization plan, the per-front inverses and the fused solve segments. X and B are not used.
  • "factorization", "refactorization" (needs a previous "factorization"): the fused factorization of the current values, then the inverses of the diagonal blocks; sets "info" (0, or the original column of the first failed pivot: a non-positive Cholesky pivot, or for "S" a zero or non-finite pivot). A failed factorization does not throw. "S" factors its wide fronts in the fused kernel under either wides strategy (a front's pivots may take both signs), so the in-kernel wide-front limit applies to it on every part: a front with more than 1023 · 32 rows below a panel raises NotSupportedError at analysis.
  • "solve": X = A⁻¹ B with the partitioned-inverse solve, for a vector (or an n × 1 matrix) on the solver's backend; X === B is allowed.

setparam! accepts "user_perm", "reordering_alg", "nd_nseps", "nd_seed" and "pivot_epsilon" (the Cholesky pivot epsilon; the fused LDLᵀ factors with epsilon 0), all read by the next analysis, and "deterministic_mode", read at every solve. The factorization is deterministic under either mode: the children of a parent add in child_list order and the panels of an in-kernel wide front apply in order, with plain adds, so factor and D are bitwise reproducible for the same values. deterministic_mode = 0 (the default) scatters the forward-solve updates with atomic adds, so the solution can differ in the last bits from run to run; deterministic_mode = 1 runs the forward sweep without atomics (every front writes its update to its own buffer, the parent pulls its children's buffers in order), and the solution is bitwise reproducible too. getparam returns "info", "inertia" ((npos, nneg); "SPD": the positive pivots and 0; "S": the counts of positive and negative entries of D), "npivots" (0 for "S": the fused LDLᵀ perturbs no pivot), the analysis data "perm_reorder_row", "perm_reorder_col", "perm_row", "perm_col", "lu_nnz", "flops", "nsuperpanels" of the inner solver, "memory_estimates" (the inner solver's, plus after the analysis the fused plans' device buffers: slot 13 the private contribution blocks of the fused factorization, slot 14 every other buffer of the plans, both added to slots 1, 2 and 5), and the configuration parameters.

Everything the pipeline does not do raises NotSupportedError, at construction or in the phase, never a silent switch to the stock path: the KernelAbstractions CPU() backend (the dependency-counter kernels spin on resident workgroups, which the CPU backend runs one after another), complex element types, the structures "HPD", "H" and "G", a uniform batch (nbatch > 1), nrhs > 1, a transposed solve (solve_mode ≠ 0), the solve sub-phases ("solve_fwd_perm" … "solve_refinement", the Schur phases), Schur mode, matching_alg ≠ "default", ir_n_steps > 0, a factor_precision other than the working precision, "S" with pivot_type ≠ 'N', a matrix without a front of width ≤ 256 for the fused solve, and every setparam! name not listed above. The options of the inner solver are checked again at the analysis. DirectSolver remains the general solver.

source

The performance knobs of the pipeline live in one bundle, a platform's tuning being a single value; how each field was measured and the per-architecture champions are in bench/ARCH-TUNING.md. default_tuning gives oneAPI, Metal and any other backend that reports no CUDA thread count per multiprocessor the AMD profile.

SparseDirectSolver.Fused.SDSTuning — Type
SDSTuning

Every performance knob of the fused pipeline in one immutable bundle, so a platform's tuning is a single value (see bench/ARCH-TUNING.md for how each field was measured and the per-architecture champions). All fields are runtime parameters — no environment variables, no recompilation.

The T34 fields: solve_workgroup is the workgroup size of the fused solve sweeps (a multiple of 32 in 32:256); wide_host_min puts the wide fronts wider than it on the host vendor path in a :merged Cholesky plan (the per-front wide strategy; typemax(Int), the default, keeps every wide front in-kernel; ignored for "S"); graph_max_n replays the refactorization and the solve sweeps from a CUDA graph when n <= graph_max_n and no vendor call falls in the captured range (0, the default, never captures).

source
SparseDirectSolver.Fused.default_tuning — Function
default_tuning(backend) -> SDSTuning

The measured per-architecture defaults — the single place they live. 2048-threads/SM NVIDIA parts (sm70/80/90): host wides, 16K budgets, fused-solve workgroup 128 (T34, measured on an H200; sm70/80 not re-measured). 1536-thread parts (sm86/89, Blackwell workstation): merged wides, 8K budgets, solve workgroup 64 (T34 did not measure them). AMD (gfx906-measured): merged, 8K, workgroup 256, spin cap 16. A backend without a CUDA resident-thread count (`maxthreadspermultiprocessor` returns 0: AMDGPU, oneAPI, Metal, the CPU backend, anything but CUDA) gets this AMD profile.

source