GPU Backends and When They Help

Sentinel can offload parts of the zone decomposition inverse solve to a GPU: CUDA on NVIDIA hardware, Metal on Apple Silicon. This page discusses what the GPU is actually used for, why the two backends accelerate different things, and how to decide whether a GPU will help your problem. For the concrete steps, see Set Up GPU Acceleration.

Where the time goes

A zone decomposition inverse solve spends its time in a few distinct places: assembling and factorizing each zone's stiffness matrix for the forward and adjoint solves, computing the material gradient, and re-solving zones during line search. Zone decomposition produces hundreds of independent subproblems per global iteration — each far too small to saturate a GPU on its own, but together exposing substantial data parallelism.

GPU acceleration in Sentinel is therefore built around batching: rather than making any single zone faster, the driver restructures the iteration so that one operation can be launched across all zones at once.

The three execution modes

Sentinel offers three execution modes for zone decomposition. All three produce equivalent numerical results; they differ only in how zones are scheduled and where compute-intensive kernels run.

Mode	Entry point	Parallelism	Hardware
Serial	zone_decomposition_solve!	None — zones solved sequentially	Any CPU
Distributed	zone_decomposition_solve_parallel!	pmap_func (e.g., Distributed.pmap)	Multi-core CPU
Batched GPU	zone_decomposition_solve_batched_gpu!	Phased BSP: threaded CPU + GPU kernels	CPU + Metal or CUDA

In per-zone mode (serial and distributed), each zone runs a complete CG optimization loop independently — forward and adjoint solves, gradient, line search — before its result is consolidated. The zones are embarrassingly parallel, which is exactly what distributed mode exploits, but nothing inside a single zone is big enough for a GPU.

The batched driver restructures this loop into synchronized phases using a Bulk Synchronous Parallel (BSP) model, so that all zones advance through each phase together:

The synchronization is what buys GPU efficiency: Phase 2 becomes a single kernel launch covering every Gauss point of every zone (see GPU Gradient Kernels for the kernel and its data layout), and on CUDA, Phase 1 becomes one batched factorization of all zone matrices. The price is lockstep execution — but since zones within a global iteration do comparable amounts of work, little is lost to waiting.

Why Metal and CUDA accelerate different things

The two backends look similar from the API but occupy different points in a precision trade-off.

CUDA offloads everything. NVIDIA GPUs execute Float64 arithmetic in hardware, and Sentinel's stiffness matrices are ComplexF64. The cuDSS library factors and solves those systems entirely in GPU memory, so the dominant cost of the pipeline — the linear solves — moves to the GPU along with the gradient kernel. This is why CUDA is the production backend for large inversions.

Metal offloads only the gradient. Apple GPUs support Float32 only, and the two workloads tolerate that very differently:

• The gradient kernel accumulates products of shape function derivatives and displacement fields. Float32 rounding errors in each Gauss point contribution average out over many Gauss points and CG iterations, so Float32 gradients are accurate enough to drive the optimizer to the same answer.

• The linear solves are not so forgiving. The stiffness matrices are complex-valued and ill-conditioned; LU factorization in Float32 produces roughly 18% solution error, far beyond what a forward solver can accept. Apple's Metal Performance Shaders enforce Float32 ("Only MPSDataTypeFloat32 is supported"), which rules out GPU linear solves on Metal entirely.

So on Apple Silicon the solves stay on the CPU in Float64, and the win comes from two other places: the Metal gradient kernel, and the batched CPU line search (CPUBatchedSolver), which drives the per-zone factorizations as a lockstep batch with symbolic reuse. The latter is the larger effect — roughly 2.5× on the 2-iteration MGH benchmark versus the per-zone fallback (see Tune Performance).

When acceleration helps

• NVIDIA GPU available — use the batched driver with CUDSSBatchedSolver. Every phase is accelerated; the full 100-iteration MGH solve runs 2.1× faster than the 10-thread CPU baseline on the same machine (CUDA Backend has the full benchmark table).

• Apple Silicon — use the batched driver with KAGradientBackend(MetalBackend()) and a CPUBatchedSolver. Expect the batched line search, not the GPU, to provide most of the speedup.

• Multi-core CPU server, no GPU — use distributed mode with BLAS.set_num_threads(1) per worker. The zones are embarrassingly parallel; plain process parallelism captures that without any GPU.

• Small problems (roughly under 20 zones) — stay with the serial driver. Batching and GPU launch overhead dominate at small zone counts, and the gradient kernel only saturates a GPU once it has thousands of work items.

Design principles

The GPU support is deliberately additive:

1. Additive, not replacement. GPU code paths sit alongside CPU paths; no CPU functionality is modified or removed, and the CPU paths remain the reference implementation for correctness.

2. Opt-in via keyword arguments. Acceleration is activated by passing gradient_backend or cudss_batched_solver; the default is always CPU-only.

3. GPU code lives in package extensions. SentinelMetalExt and SentinelCUDAExt load only when Metal or CUDA/CUDSS are imported. The core module has no GPU dependencies, and the test suite runs entirely on CPU.

4. Dispatch-based backend selection. Backends extend a small set of functions (_gpu_float_type, _to_gpu_array) rather than branching on backend identity, so adding a backend does not touch driver code.

Direct or iterative solvers on the GPU?

Sparse direct factorization is inherently sequential, which sits awkwardly on GPU hardware. An investigation into iterative alternatives found that ILU(0) preconditioning reduces GMRES to a single iteration on every zone of the MGH dataset — effectively a parallel direct solve. In practice cuDSS batched factorization still wins at Sentinel's zone sizes, but the investigation explains the trade-off and when that could change: see Iterative Solvers on the GPU.