arXiv CSopen access

Reducing Data Movement in the Galerkin Product of Block Algebraic Multigrid on GPUs

Reducing Data Movement in the Galerkin Product of Block Algebraic Multigrid on GPUs Mark F. Adams Lawrence Berkeley National Laboratory Berkeley, CA 94720 [email protected]

Abstract

arXiv:2607.28891v1 [cs.SE] 30 Jul 2026

T

The Galerkin triple product Ac = P AP dominates the recurring per-solve setup cost of algebraic multigrid (AMG). For AMG on systems of PDEs, where the null space is larger than the number of unknowns per mesh point, the product is a rectangular-block sparse matrix triple product. For 3D elasticity, the fine operator has 3 × 3 blocks, the prolongator 3 × 6 blocks, and the coarse operator 6 × 6 blocks, a shape for which no vendor sparse library provides a matrix–matrix product. We map the algorithm space of this operator, from the classical two-pass product through fused-recompute, schedule-reordered, and shared-memory-tiled variants to an inspector–executor form, under an explicit DRAM/L2 traffic model, and implement the leading variants in both portable Kokkos (CUDA) and native CUDA backends using new PETSc blocked matrix types. The model’s request-volume accounting, validated on an NVIDIA A100, predicts per level which variant moves the fewest bytes. Guided by it, a shared-memory-tiled kernel with a sorted, search-free schedule moves fewer bytes and less than half the time of the portable Kokkos team kernels on the fine-level product (10.5 versus 17.4 GB of DRAM, 45 versus 82 ms), within 2.5× of the model’s streaming floor for the full product and 1.9× on its A · P stage. We further present prolongator filtering, a new algorithm in PETSc’s GAMG that drops small blocks from the coarse grid space under a Frobenius criterion and applies a kernel-preserving projection. It reduces PtAP traffic, coarse-operator fill, and memory use and cuts the hot PtAP time 2.9× on the fine grid with iteration counts unchanged. The driving application is a fully GPU-resident blocked pipeline in PETSc: finite-element assembly writes directly into the blocked device matrix, and the AMG setup, Galerkin products, and solve all operate on primary blocked data with no scalar expansion and no operator-sized device–host transfers in the recurring phases.

Keywords: algebraic multigrid, matrix triple product, sparse matrix–matrix multiplication, block sparse matrices, smoothed aggregation, performance modeling, CUDA, Kokkos, PETSc

1

Introduction

Algebraic multigrid (AMG) solves the sparse linear systems of elliptic PDEs in near-linear work and is a mainstay of extreme-scale simulation. Its cost structure splits into a solve phase dominated by sparse matrix–vector products (SpMV), a matrix setup phase dominated by sparse matrix–matrix products, the Galerkin triple product Ac = P T AP (hereafter PtAP), and a mesh setup phase dominated by graph work and a matrix-matix product in smoothed aggregation [16]. For nonlinear and time-dependent applications the coarse grid spaces (columns of P ) are reused (re-smoothing is rare), but the fine and coarse operator values change, so the numeric PtAP recurs on every Newton or time step and becomes a recurring cost of the solve itself [2]. 1

For systems of PDEs with large null spaces AMG gives PtAP a structure that vendor sparse libraries do not serve: the null space of 3D elasticity (six rigid-body modes) makes the prolongator P rectangular-blocked: 3 × 6 blocks between 3 × 3-blocked fine operators and 6 × 6-blocked coarse operators. Vendor library support is at best partial: rocSPARSE [4] stores rectangular blocks (its GEBSR format) and provides sparse-times-vector and sparse-times-dense products on them, cuSPARSE [10] restricts block-sparse compute routines to square BSR, and Kokkos Kernels [12] is square-block only. No vendor library offers the sparse–sparse product over rectangular blocks that the Galerkin triple product requires. In previous work [1] we introduced a portable Kokkos-backed blocked matrix type in PETSc [5] (MATBAIJKOKKOS) and a conversion-free blocked GAMG setup. Profiling of that path shows the PtAP numeric is L2-bandwidth-bound: on an A100 its blocked A · P stage runs at 89 % of L2 bandwidth while DRAM sits at 20 % and the floating-point pipes below 20 % (Section 3.5), so its dense block contractions are cheap relative to the bytes they move. That observation frames this paper’s question: across the space of algorithms for the rectangularblock Galerkin product, which methods are most effective at reducing data movement cost? The space runs from the classical two-pass product (W = AP materialized, then P T W ) through fusedrecompute forms that never materialize W , traversal reorderings that improve the cache locality of the fused form’s redundant reads, and tiled kernels that stage coarse-row accumulators and W -panels in shared memory, to an inspector–executor form that flattens the numeric phase into a branch-free stream of block contractions. A companion question, whether FP64 tensor cores can accelerate the dense block contractions at the core of every variant, we answer analytically, in the negative: one FP64 tensor-core instruction on the A100 computes an 8×8×4 product (256 multiply–adds), the largest contraction our block shapes admit per fragment is 6×6×3 (a 42 % utilization, obtained by packing two A-blocks that multiply the same P -block), and 0.42 of the 19.5 TF/s tensor-core peak falls below the 9.7 TF/s FMA peak [9], before counting the cost of staging operands into the fragment’s fixed thread-to-element layout. Repacking these small dense blocks for tensor cores therefore cannot pay at full precision; their reduced-precision paths retain a large advantage even at poor fragment utilization, an accuracy-for-cost trade we leave as future work (Section 7). We predict per-level performance with an explicit traffic model, implement the promising variants in two equally supported back-ends, portable Kokkos and native CUDA (explicit shared memory, atomic-free owner-computes contracts, sorted schedules), and report measured-vs-model results per hierarchy level. The kernels stay general in that there is no symmetry exploitation, and the scalar case bs=1 is served by the same type: the Galerkin product runs the same native kernels, and because at bs=1 the block-CSR arrays are a scalar CSR matrix, the matrix–vector product dispatches to the vendor CSR kernel over the resident device arrays (Section 6.5), so the blocked type is safe as a general-purpose default. This paper makes two further contributions. First, prolongator filtering, a new algorithm in PETSc’s GAMG (the -pc gamg prolongator filter option): small blocks of the coarse grid space are dropped under a Frobenius criterion with a kernel-preserving projection that restores the near-null space. Because filtering perturbs the smoothed prolongator, it must be accounted as a trade: reduced work and memory against possible convergence degradation. In practice significant cost reductions are realized with negligible degradation of convergence rates (Section 6): nnzb(P ), and with it every term of the PtAP traffic model shrinks. Second, the driving application is a fully GPU-resident blocked pipeline in PETSc: finite-element assembly emits block COO (Coordinate format)directly into the blocked device matrix, and aggregation, prolongator construction, smoothing, filtering, PtAP, and the solve all operate on primary blocked data, with no scalar AIJ expansion and no device–host round trips in the recurring phases. Contributions.

2

1. A taxonomy and DRAM/L2 traffic model of the rectangular-block PtAP algorithm space (twopass, fused-recompute, scheduled, SMEM-tiled, inspector–executor), with per-level predictions and measured validation on A100 GPUs. 2. Low-level CUDA and Kokkos implementations of the leading variants using PETSc’s new blocked matrix types, using an atomic-free owner-computes contract and tile-sorted schedules, with a Kokkos-vs-CUDA comparison. 3. Prolongator filtering in PETSc: a block-Frobenius filter with a kernel-preserving projection that is evaluated as a trade of work and memory against convergence. 4. A fully GPU-resident blocked FEM-to-AMG pipeline in PETSc, from finite-element block-COO assembly through the linear solve, nonlinear solve and time integration [3]. The remainder is organized as follows. Section 2 reviews smoothed aggregation, block structure, and the driving application. Section 3 develops the algorithm space and traffic model. Section 4 presents prolongator filtering. Section 5 describes the CUDA and Kokkos implementations. Section 6 reports results on NVIDIA A100 GPUs. Section 7 concludes.

2

Background

2.1

Algebraic multigrid by smoothed aggregation

Multigrid methods damp high-frequency error with a cheap smoother and damp the smooth remainder on a coarser level, recursively [14, 2]. Algebraic multigrid builds the coarse levels from the matrix itself. The smoothed-aggregation (SA) variant [16] aggregates strongly-coupled nodes, builds a tentative prolongator P̃ that spans the near-null space exactly on the coarse level, and improves it with one damped-Jacobi step, P = (I − ωD−1 A)P̃ . The coarse operator is the Galerkin product Ac = P T AP , and the construction recurses. The near-null space is the null space of the operator without its essential (Dirichlet) boundary conditions and is provided by the user; SA is thus not a pure algebraic solver. For 3D elasticity these are the six rigid-body modes, and the theory requires the coarse space to span this space.

2.2

Block structure and the rectangular-block Galerkin product

Vector PDEs associate several unknowns with each mesh point, so the fine matrix carries dense 3 × 3 blocks (3D elasticity) sharing one column index per block. The SA coarse space carries one degree of freedom per null-space mode per aggregate: six in 3D. Consequently the block structure is not uniform across levels: P has rectangular 3 × 6 blocks on the finest level and the coarse operators carry 6 × 6 blocks. The Galerkin numeric is thus a chain of small dense contractions, (6×3) · (3×3) · (3×6) on the finest level and (6×6) · (6×6) · (6×6) below it, threaded through sparse block-graph traversal. Storing blocks natively amortizes one index over bsr bsc values and raises arithmetic intensity. The companion paper [1] develops the blocked storage argument and the conversion-free blocked GAMG setup that this paper builds on.

2.3

Driving application and test problem

The driving application is 3D linear elasticity on a cantilever beam, discretized with Q2 hexahedral finite elements (PETSc snes/tutorials/ex56k, ex56cu), clamped at one end with combined torsional and transverse end load, solved with conjugate gradients preconditioned by AMG. The 3

Figure 1: Deformed cantilever beam under combined torsional and transverse end load, colored by displacement magnitude, on an isotropic 80 × 8 × 8 Q2 mesh (139K degrees of freedom, unit-cube elements). The beam extends along x, and the clamped face is at x = 0, with the y and z axes drawn through the origin. The torsional component of the load is increased here to make the twist visible. beam is 10 × 1 × 1 with 10k × k × k elements, so the elements are isotropic (unit cubes), and the mesh is sized large per GPU so that the PtAP working sets exceed the L2 cache, the regime the traffic model addresses. The single-device base case fixes k = 16: 160 × 16 × 16 elements and 1,045,440 displacement degrees of freedom, with the convergence gate of Section 6.1 (19 CG iterations, tip displacement 1.23166). A dedicated study of element anisotropy and aggregation quality is future work (Section 7).

2.4

The full-GPU blocked pipeline

All measurements run inside a fully GPU-resident blocked pipeline in PETSc (Figure 2) with two equally supported device back-ends: one built on Kokkos (MATBAIJKOKKOS) and one directly on CUDA (MATBAIJCUDA). The back-end is selected at run time by a command-line option (-dm mat type baijkokkos or -dm mat type baijcuda, with the matching -dm vec type), with no change to application code. In the pipeline, a device backend for PetscFE assembles the Jacobian and emits block COO directly into the blocked matrix, using the blocked coordinate-assembly primitive of the companion paper [1]. GAMG’s aggregation (a device maximal-independent-set coarsener), tentative prolongator, prolongator smoothing, filtering (Section 4), and Galerkin products all operate on the block format, and the preconditioned solve (block-Jacobi smoothing on bs × bs blocks, Chebyshev/CG) never leaves the device. The Kokkos back-end was developed first [1], which allowed development on a CPU and facilitated a fairly mechanical translation to correct native-CUDA matrix, finite-element, and GAMG back-ends. Primary blocked data flows from element integration to the coarsest AMG level without scalar expansion, and Section 6.6 verifies with the device–host transfer counters that the recurring phases run without transfers.

3

The Rectangular-Block PtAP Algorithm Space

3.1

Notation and granularity

Let mF /mC be fine/coarse block rows; ZA = nnzb(A), ZP = nnzb(P ), ZW = nnzb(W ), ZC = nnzb(Ac ) the block nonzero counts, with per-row averages nA = ZA /mF , nP = ZP /mF , nC = 4

residual graph work

host

mesh (DMPlex); FE tabulation B, D maps, B, D (one-time)

libraries

device

FEM assembly Ae integration, block-COO scatter

A0

one-time AMG mesh setup graph, MIS aggregation, P̃ , smooth, filter: emits Pi per level

native blocked device kernels (MATBAIJCUDA / MATBAIJKOKKOS)

Pi A0

re-assemble A (FEM stage) A0 matrix setup Ac = P TAP per level (variants V0–V4)

KSP/SNES/TS control (iteration logic) scalars

coarse LU (option) one-time factor

x0

blocked SpGEMM (native); bs=1: cuSPARSE, KK

KSP solve CG on blocked vectors

AMG PC apply Chebyshev + pb-Jacobi smoother, P , R, coarse solve

SpMV: native blocked, bs=1 → vendor CSR; coarse solve: cuDSS or host LU

recurring Newton / time-step solve

Figure 2: The blocked FEM–AMG workflow. Stages run left to right; the middle row is the deviceresident pipeline, the top row the host work, and the bottom row the sparse linear-algebra realization of each stage, with the vendor alternatives where they exist. The host row is nearly empty by design: the Krylov iteration control runs on the host on scalar reduction results, the dashed edges are the only host–device transfers in the recurring phases (Section 6.6), and no operator-sized transfer crosses the bus. DMPlex/PetscFE build the closure maps and the reference-element tabulation (basis B and derivative D arrays) on the host, uploaded once at setup. The mesh-setup stage abbreviates the strength graph, aggregation, tentative-prolongator, smoothing, and filtering steps and emits the prolongators Pi of every level, with a residual of one-time graph work on the host. The dashed enclosure marks the recurring Newton/time-step solve: each step re-assembles A (the shadow FEM stage, running the same device assembly as stage one). ZC /mC ; and BA , BP , BW , BC the bytes of one dense block of each operand (72, 144, 144, 288 bytes at fine-level block sizes; all 288 bytes on coarse levels; 8 bytes at bs=1). Two derived quantities recur. First, nP , the average number of coarse aggregates each fine block row couples to after prolongator smoothing and filtering, turns out to control the choice of algorithm (Section 3.4). Second, the fill growth of the intermediate product: w = ZW /ZP measures the growth in coarse coupling from applying A to P (a row of W = AP unions the P -patterns of an A-neighborhood, reaching one A-layer beyond P ’s support). We measure ZW directly, with a counter in the product symbolic phase: w = 4.4, 7.3, and 3.3 on the three products of the test hierarchy (Table 2). The model’s conclusions are insensitive to w (the fused-variant crossover of Section 3.4 moves only from n∗P = 2.4 to 2.5 as w ranges from 2 to the measured 4.4) because w raises both the two-pass variant’s avoidable traffic and the fused variant’s recompute. Three granularities structure every variant: • the b-block : the dense dof block, the micro-kernel operand held in registers; • the tile: a contiguous set of coarse block rows (ideally whole aggregates) plus their fine-row fan-out, sized to fit shared memory plus an L2 slice; • the fusion grain: whether W = AP is materialized globally (two-pass), consumed per-nonzero or per-row as produced (fused), or staged per-tile (tiled fusion).

3.2

Traffic accounting: floors and request volumes

Two counting rules generate every prediction in this section. The stream rule: a pass that touches every block of an operand exactly once costs its full byte size in DRAM traffic, and an intermediate 5

that is written and later read back costs double. Summing streams gives each variant a DRAM floor, the bytes it must move if every distinct byte is fetched once. The touch rule: in a sparse traversal every block touch is a memory request, served by whatever tier of the memory hierarchy holds the block, and only misses become DRAM traffic. Counting touches gives a request volume, and the tier that serves the touches determines which bandwidth limits the kernel. The distance between an algorithm’s floor and its request volume is redundancy, and the organizing question of the design space is where each algorithm puts its redundancy. We study five classes of algorithms, labeled V0–V4 and defined in Section 3.3: • V0, the two-pass baseline: the intermediate W = AP is materialized in DRAM; the redundancy is the write and read-back of W and the gathered requests for P rows. • V1, fused recompute: W is never formed. The redundancy becomes repeated reads of the A and P rows, served by cache. • V2, a tiled family built on V1’s fused traversal: V2a reorders the traversal, V2b holds output tiles in shared memory, and V2c additionally stages per-tile W -panels on chip (reintroducing V0’s intermediate, but tile-local), moving the redundancy progressively into on-chip staging. • V3, inspector–executor: a schedule precomputed at symbolic time replaces all traversal and search work. The redundancy is the streaming of the schedule itself. • V4, transpose-free outer-product scatter variant of V0: the redundancy moves to the write side, as atomic read-modify-writes into Ac . Table 1: A100 machine parameters, from the whitepaper [9] and microbenchmark literature (L2). Ratio checks that need timings free of GPU launch and occupancy effects use a host control, an Apple M3 (six P-core threads, 85 GB/s measured triad DRAM bandwidth, 16 MB shared-P-cluster L2), introduced in Section 3.5. NVIDIA A100-SXM4-40GB DRAM bandwidth L2 cache L2 bandwidth Shared memory FP64 peak

1555 GB/s 40 MB ∼5 TB/s 164 KB/SM × 108 SMs 9.7 TF/s FMA / 19.5 TF/s tensor core

The A100’s memory tiers (Table 1) set the terms of the analysis. The DRAM : L2 : SMEM bandwidth ratio is roughly 1 : 3.2 : 12.5: the binding resource is L2 bandwidth, a kernel whose request volume is large gets pinned at the L2 tier even when its DRAM traffic is modest, and moving redundant requests from L2 into shared memory or registers buys another factor of ∼4.

3.3

The variants

V0: two-pass baseline. Materialize W = AP in DRAM, then form Ac = R W with R = P T stored explicitly: the transpose structure is built once at symbolic time and its numeric refresh is a pure value gather, measured at ∼5% of the fine-level numeric. (We prefer the cached explicit transpose to an implicit one because coalesced row access beats atomic column-scatter on GPUs, and library SpGEMM kernels require CSR inputs.) The DRAM floor is TV 0 = BA ZA + 2BP ZP + 2BW ZW + BC ZC . 6

The A, P/R, and Ac streams are irreducible (every algorithm reads the operands and writes the result), so the write-then-read of W is the largest avoidable term (49 % of fine-level bytes for the hierarchy of Table 2 at the measured w), and it grows with w nP while the A stream is fixed. V0 is the correctness anchor and the model’s reference stream. Measured, however, V0 runs far from this floor: profiled on the A100, its blocked A · P stage saturates 89 % of L2 bandwidth while DRAM sits at 20 % and the FP pipes at 19 %. The reason is the touch rule: the gather of P rows generates ∼nP ZA block requests, an order of magnitude above the stage’s DRAM volume, plus hash work to place each contribution in W . That measurement motivates the rest of the space: the binding resource is L2 requests and per-contribution discovery work, not DRAM streams or flops. V1: fused recompute. Never materialize W : the owner of coarse row I walks R-row I → A-rows → P -rows and accumulates directly into Ac row I. At runtime -mat product algorithm names a variant explicitly (memory for V1, speed for V0). Left at its default, the implementation chooses (Section 3.5). Fusing the Galerkin product on GPUs traces to Bell, Dalton, and Olson [7]. The rectangular-block setting changes the costs because every touch moves a whole dense block. V1’s DRAM floor drops below V0’s (BA ZA + 2BP ZP + BC ZC , no W stream) and, equally important, the W allocation disappears: about a gigabyte per device at our test sizes (Table 3), growing linearly with problem size. The price is redundancy with a precise coefficient: fine row i is walked once per coarse row that contains it, i.e. nP (i) times, so A-block touches are ∼nP ZA , P -block touches ∼n2P ZA , and the flop count is exactly nP × V0’s. The size of the working set determines which tier serves those touches: one aggregate’s rows of A and P plus its Ac row total ∼0.3 MB at our fine-level constants (too large for shared memory, but small against a 40 MB L2, which holds roughly 130 such working sets), so the redundant re-reads are served by L2, and V1 is, at these working-set constants, limited by the L2 request rate. On CUDA an owner-computes contract makes V1 atomic-free and (for register-resident rows) search-free. V2a: scheduling only. Identical data structures and touch counts to V1, but aggregate-contiguous coarse numbering (optionally an RCM ordering of the coarse graph) makes consecutive owners share fine rows, so redundant touches find their blocks resident. By construction V2a can change which tier serves a touch, never the touch count: it is a zero-format-change experiment that cleanly separates scheduling effects from capacity effects, and its measured value below is exactly that separation. V2b: shared-memory Ac -tile accumulators. A thread block owns a tile of coarse rows whose Ac blocks live in shared memory. A symbolic pass emits a per-tile contribution list sorted by output slot, eliminating per-hit binary search and all global atomics, and the tile is written out once, coalesced. V2b removes the output-side discovery work but input touches still go to L2. V2c: W -panel staging (full tile fusion). V2c merges the central ideas of V0 and V1: it reintroduces V0’s materialized intermediate, but per tile and on chip. The tile’s W -panel (the rows of AP restricted to the tile’s fine fan-out) is computed once into shared memory or an L2 slice and consumed in place, so V1’s ×nP recompute disappears without W ever reaching DRAM. Traffic approaches TV 2c ≈ (1 + h) (BA ZA + BP ZP ) + BP ZP + BC ZC , where h is the halo fraction of fine rows shared by neighboring tiles (we carry h = 0.25): each operand is read approximately once per tile. The first term is the fine-side reads of A and P that build the tile’s W -panel. The second is the read of R = P T in the reduction Ac = R W (the same 7

bytes as P , read a second time in transposed layout, with no halo factor because each coarse row belongs to exactly one tile), and the last is the Ac write. Feasibility at our fine-level constants: a one-aggregate tile’s W -panel plus one Ac row occupy ∼110 KB, inside the A100’s 164 KB shared memory, at the cost of occupancy, with an L2-resident panel as the fallback tier. Tile size is a tunable per level. V2c holds the lowest floor of any variant (Table 3), yet we do not implement it: the measured tiled kernels that would host its staging turn out to run latency- and occupancy-bound, well below the DRAM roofline, so the floor advantage does not convert to time on this hardware. Section 6.3 closes that question by measurement. V3: inspector–executor. The symbolic phase flattens the entire numeric into one tile-sorted stream of (R-block, A-block, P -block, output-slot) quadruples. With tile-local 16-bit offsets each quadruple packs into 8 bytes, and there is one quadruple per block multiply, so the schedule has cnt ≈ nP ZA + nP ZW entries. The numeric is then a branch-free sweep of fixed-shape block contractions. The traffic accounting per numeric pass: the operands are read once each (V1’s floor, since the tile-sorted order keeps the redundant operand touches in cache) plus the sequential read of the schedule itself. At our fine-level constants the schedule stream (1.09 GB at 8 bytes per entry) is about half the write-and-read of W that it replaces (2.09 GB at the measured w), so V3’s total lands between V1’s floor and V0’s stream: a modest traffic gain. What it buys beyond that is the elimination of all discovery work (no hash, no search, no atomics, no graph traversal) with the schedule streaming sequentially at full DRAM rate. This matches the dominant production regime: in nonlinear and time-dependent solves the hierarchy structure and the prolongator values are typically reused while only A’s values change per step (re-aggregation is rare), so one symbolic inspection amortizes over many numeric refreshes. Our cold/hot cost accounting (Section 6) matches this regime split directly. The implementation improves on the 8-byte estimate. The schedule is sorted by output slot, so all contributions to one slot are consecutive; storing each slot once, with the length of its run, removes the slot from the per-entry record, which then holds only the two operand offsets, and because those offsets are relative to the tile’s operand window they fit in 16 bits, resulting in 4 bytes per entry. Section 6.3 reports the measured host result: the schedule occupies 0.73–0.88× the W value storage it complements, and its construction cost is repaid within the first numeric pass (the measured reuse crossover is N =1). V4: transpose-free outer-product scatter. A single matrix–transpose–matrix primitive MTM(A, B) = AT B (PETSc’s MatTransposeMatMult) applied twice gives W ′ = MTM(A, P ) = AT P and Ac = MTM(W ′ , P ) = P T AP , exact for general nonsymmetric A. R is never formed, saving its storage and its per-step value refresh. The natural device realization is an outer-product scatter: iterate over fine rows k, read W ′ (k, :) and P (k, :) exactly once each (two perfectly coalesced streams) and scatter the |P (k, :)| · |W ′ (k, :)| block contributions into Ac by atomic accumulation. This trades the gather variants’ read amplification for write amplification: ∼nP ZW atomic read– modify–write updates of BC -byte blocks, each of which reads and writes its block. In its data movement V4 is therefore a variant of V0: a two-pass algorithm with a materialized intermediate and the same stream floor (it reads P once per pass instead of P then R). Its advantages are structural. Parallelism ranges over fine rows in every pass, which is uniform across levels and eliminates the coarse-level occupancy pathology measured below (kernels that parallelize over coarse rows have too few work units to fill the GPU on the small coarse matrices, Section 3.5). Finally, V4 connects the two ends of the taxonomy: replacing its atomic scatter with a precomputed slot-sorted schedule removes the write amplification at the cost of the schedule’s stream traffic, and the result is exactly V3’s executor. A cost of the atomic form is a nondeterministic summation order. Section 6.3

8

measures V4 on both backends.

3.4

Model predictions

Table 3 evaluates the model on the measured hierarchy of Table 2. Three predictions follow. Table 2: Measured hierarchy of the 1.05M-DOF single-device test case (Section 3.5): three Galerkin products, block units. P is the smoothed, filtered (-pc gamg prolongator filter 0.03) prolongator. ZW is the measured block count of the intermediate W = AP , and w = ZW /ZP its fill growth. product A/P blocks l0 (fine) 3×3 / 3×6 l1 6×6 / 6×6 l2 6×6 / 6×6

mF

ZA

nA

ZP

nP

ZW

w

mC

ZC

nC

348 480 21.23M 60.9 1.66M 4.76 7.27M 4.38 7 243 627K 86.6 7 243 627K 86.6 51.6K 7.12 376K 7.29 1 125 127K 112.9 1 125 127K 112.9 8 970 7.97 29.4K 3.27 98 3 770 38.5

Table 3: Model-predicted traffic and footprints (GB) per Galerkin product for the hierarchy of Table 2, using the measured per-product ZW (w = 4.38, 7.29, 3.27; h = 0.25). “Requests” count block touches (read side for V1, atomic read-modify-write side for V4, counted at twice the block bytes for the read and the write). Floors count each distinct byte once. l0

l1

l2

V0 stream (floor) V1 floor (no W ) V2c floor V4 floor V1 read requests V4 write requests

4.28 2.19 2.63 4.28 76.6 19.9

0.46 0.25 0.30 0.46 10.4 1.54

0.060 0.043 0.053 0.060 2.6 0.14

W footprint V3 schedule footprint (8 B/entry)

1.05 1.09

0.108 0.057

0.008 0.010

18.4 / 87.5

3.1 / 22.0

0.5 / 4.3

flops V0 / V1 (GF)

First, the floors do not decide: the redundancy does. V1 has the lowest floor on every level, yet its read-request volume is 35× its floor on the fine level (76.6 vs. 2.19 GB) and its flop count is nP × V0’s. Whether V1 can come out ahead on time is therefore an nP question: equating V1’s compute floor with V0’s DRAM floor gives a crossover at n∗P ≈ 2.5 for the fine-level constants of Table 2.1 Across our measured hierarchies the band is n∗P ≈ 2.2–2.5. The measured nP = 4.76 sits above the whole band, so the model predicts V1 loses the fine level by a factor of order nP , and by more on the coarse levels, where nP = 7–8. V1’s value at these nP lies in the deleted W footprint and in low-nP hierarchies, and nP is an algorithmic choice, not a constant. Stronger prolongator filtering reduces it (Section 4), and unsmoothed aggregation, a valid point in the design space that uses the tentative prolongator (P = P̃ , hence nP = 1), sits far below the crossover, where the model predicts V1 is the faster variant. 1

Both sides of the equation depend on nP : ZP = nP mF and ZW = w nP mF , so V0’s flop count F0 (nP ) is quadratic in nP and V1’s compute floor nP F0 (nP ) cubic, while V0’s DRAM floor TV 0 (nP ) is affine in nP . Solving nP F0 (nP )/(9.7 TF/s) = TV 0 (nP )/(1555 GB/s) with the fine-level constants of Table 2 at the measured w = 4.38 gives n∗P = 2.5 (2.4 at w = 2: the crossover is nearly insensitive to w). Holding F0 at its measured value in Table 3 would instead give 1.5.

9

Second, scheduling alone has limited value. V2a cannot reduce touch counts, so its best case is converting DRAM refetches into cache hits. If the redundant touches already hit L2, as the working-set arithmetic above predicts, the model expects a real but marginal traffic delta and little wall-clock change. A small V2a effect is itself information: it means the redundancy must be staged away (V2b/c), not reordered. Third, the achievable floor belongs to tiled fusion. V2c holds the lowest reachable bytes on every level (unlike V1’s floor, which its own request volume prevents it from approaching) while retaining the no-global-W memory saving. V4 matches V0’s streams while moving its amplification to the write side (19.9 vs. 76.6 GB of fine-level requests) and parallelizing uniformly over fine rows. V3 converts the whole numeric into schedule-driven streaming for the reuse-dominated production regime.

3.5

Model vs. measured

We validate the model’s ordering predictions with the first two variants implemented in the Kokkos and CUDA blocked backends, V0 and V1: the isotropic elasticity base case at 160×16×16 Q2 elements (1,045,440 DOF) on one A100 (Table 1), with the hierarchy of Table 2; every configuration converges identically (19 CG iterations). We record GPU-fenced end-to-end MatPtAPNumeric times and per-kernel DRAM/L2 bytes with Nsight Compute. As a control, the same hierarchy also runs on a host machine (an Apple M3, six P-core threads, per-product stage timers), whose timings carry no GPU launch or occupancy effects. The host control’s absolute times sit more than two orders of magnitude above its DRAM floor (the host path serializes team loops and pays the SpGEMM hash), so it enters the analysis only through variant ratios, never absolute times. Table 4: Measured A100 wall time per Galerkin product, V0 (two-pass) vs. V1 (fused), with the per-product nP of Table 2. Per-product times are Nsight Compute kernel durations (profiler replay, and their sum differs from the end-to-end total by ∼15%). Totals are GPU-fenced end-to-end MatPtAPNumeric. product

nP

V0 (ms)

V1 (ms)

ratio

l0 l1 l2

4.76 7.12 7.97

82 12.3 4.1

570 106 139

7.0 8.6 33.9

87

684

7.9

total

Table 5: A100 measured bytes per Galerkin product (Nsight Compute, all kernels in the numeric product), against the model floors and request volumes of Table 3. V0 (GB) product l0 l1 l2

V1 (GB)

DRAM

L2

DRAM

L2

L2 hit

20.3 1.9 0.2

199 16.7 2.8

184.7 5.4 0.11

1246 78 16.6

89.1 % 95.1 % 99.4 %

The fused variant loses by a factor of order nP . Table 4 gives the central result: at fine-level nP = 4.76, above the model’s crossover band, V1 is 7.9× slower than V0 end-to-end (684 vs. 87 ms). 10

The host control reproduces the verdict with per-product V1/V0 ratios that track nP level by level (6.1×, 7.3×, 10.7× against nP = 4.76, 7.12, 7.97), and the A100 fine level matches at 7.0×: the same recompute cost, with the same coefficient, decides the outcome on both memory systems (the A100 and the M3 host control), so the prediction is a property of the algorithm, not of the GPU. The bytes agree (Table 5): V1’s fine-level DRAM traffic is 9× V0’s (184.7 vs. 20.3 GB): the fused variant, whose floor is the lowest in Table 3, moves the most DRAM bytes, because request volume, not floor, governs. Model validation and model limits. The model says V1’s redundant touches should land in L2 and its DRAM traffic should be the L2 miss stream. Measured, the L2 hit rate rises exactly as the per-level working set shrinks into the 40 MB L2 (89.1 %, 95.1 %, 99.4 %: hit rates, distinct from the V0 kernel’s 89 % L2 bandwidth utilization quoted in Section 1), and the implied miss stream (1 − hit) × L2 reproduces 70–90 % of the measured DRAM traffic on every level (the remainder is write-backs). The V1 fine kernel meanwhile runs at 324 GB/s of DRAM and ∼2.2 TB/s of L2 bandwidth (saturating neither), so with its requests largely absorbed by L2 it presents as latencyand occupancy-limited rather than bandwidth-bound. On the absolute scale the touch model is a lower bound: measured L2 bytes exceed the block-touch volume of Table 3 by roughly an order of magnitude (sector granularity, index and search traffic, and accumulation read-modify-writes that the operand-touch count excludes). The model is reliable for ratios, orderings, and crossovers, and those are what the measurements confirm. Scheduling alone is insufficient. The V2a experiment behaves as predicted. On the A100, RCM reordering of the fused kernel’s coarse rows produces a real but marginal traffic improvement on the fine level (DRAM 184.7 → 174.4 GB, −5.6 %, L2 hit 89.1 → 89.8 %) and is wall-neutral (569.9 vs. 574.7 ms). On the host control, all three orderings (natural, minimum-fine-index, RCM) agree within run-to-run variance. Part of the explanation is that the natural ordering is already good: aggregate roots are discovered in fine-row order, so consecutive coarse rows already sweep nearly contiguous fine neighborhoods on this regular mesh, and RCM’s level-set renumbering has little structure left to recover. Scheduling alone is therefore insufficient (not useless, and worth revisiting on unstructured meshes, since the good natural ordering is a property of this regular mesh), and by the model’s dichotomy the redundancy must be staged away (V2b/c), not reordered. Coarse levels fail differently on the GPU. The A100’s l2 product is pathological for V1: 139 ms, more than the much larger l1 product’s 106 ms, at a 99.4 % L2 hit rate on a matrix of only 1 125 block rows. This is not traffic: launching one team per coarse row leaves most of the 108 SMs idle, and long serial per-team loops do the rest. The host control shows normal scaling on the same products (l2 ≪ l1), so the inversion is a GPU occupancy artifact of the launch shape, not a property of the algorithm. Two design consequences follow: any fused or tiled variant needs a different launch geometry on coarse levels (multiple teams per coarse row, or L2-resident accumulation), and V4’s fine-row parallelism, uniform across levels, eliminates exactly this failure mode. Per-level selection and what follows. The measured verdict is per-level and nP -dependent, not global. Among the variants measured in this section (V0, V1, V2a), at production filtering levels (nP ≈ 5–8) V0 remains the reference on every level. V1’s case is memory capacity (no W ), low-nP hierarchies (strong filtering, tentative prolongators), and its role as the correctness contract for the fused family. The forward path the model and measurements select jointly: V2b/c tiling on the fine level, where staged fusion holds the lowest reachable floor; a different launch shape on coarse 11

levels, where the occupancy evidence, not the traffic model, governs; and the V3 schedule for the reuse-dominated production regime, with V4 addressing the coarse-level geometry by construction (V2b, V3, and V4 are measured in Section 6.3). The same coarse-level evidence bears on GAMG’s processor consolidation (-pc gamg process eq limit), which gathers coarse operators onto fewer processes as levels shrink. The driver of consolidation is MPI message traffic, which is outside our model: on coarse levels the per-process work is so small that idling processors loses little compute while reducing the number of messages and the latency-dominated collectives. The measurements above supply the work side of that trade (coarse products are far from any bandwidth limit, so concentration costs essentially nothing) and suggest a rule for selecting the limit: each level should operate at the knee of its strong-scaling curve, before the plateau where added processes contribute only messages. To first order the knee is a fixed amount of work per process, not a property of the level, which is why a single equations-per-process value serves the whole hierarchy. Making the selection quantitative requires a message-cost model [6] joined to the traffic model, which we leave as future work. The product-variant choice itself should not burden the user: each level’s Galerkin product is a separate operator product, so when -mat product algorithm is left at its default the implementation selects the variant per level from the level’s own local shape (block size, row count, and memory footprint) following the measurements above. Naming a variant explicitly forces it on every level, which is how the comparisons in this paper are produced. Asking users to choose per level would be a significant burden for a small gain, and the local operator each process owns on a coarse level is itself a product of the consolidation the previous paragraph describes, which the automatic per-product choice adapts to without user involvement.

4

Prolongator Filtering

Smoothing the tentative prolongator densifies it: P = (I − ωD−1 A)P̃ inherits the fill of AP̃ , and that fill is squared through the Galerkin product. We introduce prolongator filtering in PETSc’s GAMG: after smoothing, blocks of P whose Frobenius norm falls below a threshold (relative to their block row) are dropped, and a kernel-preserving projection restores the null space (the row-sum, or lumping, correction of scalar AMG generalized to blocks). Dropping small entries with a compensating lumping to control operator complexity is established practice in smoothed aggregation [16, 15]. Lumping, however, preserves only the row sum, the action on constant vectors, that is, the translational modes. The projection here preserves the entire null space, including the rotational modes. To our knowledge this local post-filtering projection is new. Energy-minimizing prolongators [11] also preserve the full null space, but by imposing it as a constraint in a global minimization over a prescribed sparsity pattern rather than by a local correction after dropping. Filtering approximates the mathematics: dropping blocks perturbs the smoothed prolongator, so the correct accounting is a trade of possible convergence degradation for reduced work and memory. Measured, the trade is strongly favorable: on the base case the hot PtAP time falls 2.9× with iteration counts unchanged (Table 9, including the unfiltered row). Meanwhile nnzb(P ) and nnzb(W ), which enter every term of Section 3.3’s traffic model, and the fill of the coarse operators all shrink; at large enough problem sizes (the eight-GPU weak-scaling case of Section 6.4, and the scalar case of Section 6.5) the unfiltered Galerkin product exhausts GPU memory, so there filtering is also a capacity requirement. Filtering is one instance of a principle that recurs in this paper (Sections 1 and 7): spending controlled accuracy in the setup operators to buy work and memory.

12

5

Implementation

5.1

CUDA kernels: owner-computes, sorted schedules, tiles

The CUDA blocked types assign each output block row to one owner. In the single-thread-per-row reference kernels that ownership makes the two-pass and fused variants atomic-free with plain register accumulation, and those kernels are kept as the byte-reproducible accumulation references; the default team forms split the owned row across a sub-warp, register-stage the operand blocks ahead of the accumulation (Section 6.3 attributes the fine-level traffic to this staging), and accumulate with atomic adds (V4’s outer-product scatter is atomic by construction, Section 3.3). The tiled variant V2b replaces row ownership with output-slot ownership. A symbolic pass walks the same nested (i, A-col, P -col) loops as the two-pass kernel and stably buckets each block contribution by its output block-slot in C, emitting a per-slot contribution list (the two lists of A- and B-block operand offsets, ordered by output slot then by the two-pass iteration order). The numeric kernel then owns a contiguous tile of cap output block-slots per thread block: each thread accumulates one slot’s contribution run in a register block, with no per-hit binary search (the slot is precomputed) and no global atomics (each slot has one owner), stages the finished block into shared memory, and the block writes the tile out once, coalesced. Because every output block sums its contributions in exactly the two-pass order, V2b reproduces V0’s per-block accumulation order. A host-side bitwise check confirms this (on device, nvcc’s FMA contraction differs between the two separately compiled kernels, a floating-point accuracy effect in the last bits, and the device setup itself is not bit-reproducible run-to-run at the 10−12 level, so end-to-end agreement is to round-off rather than bit-for-bit). Tiling by output slot rather than by row means a fat coarse row simply spans several tiles: the shared-memory accumulator never exceeds cap, so there is no fat-row fallback and nothing is truncated. The one fallback is to the two-pass V0 kernel when a single output block exceeds the register accumulator bound (36 scalars, i.e. blocks larger than 6 × 6). The whole PtAP runs as the two-pass product W = A P then Ac = R W with each A · B tiled. The two schedules are built once and cached, and only operator values refresh under hierarchy reuse. The tile width cap is set so a tile’s output blocks fit the static shared-memory budget: cap = ⌊S/(bsr bsc · 8)⌋ for an FP64 build (Table 6). At S = 48 KB (the per-block shared-memory limit CUDA imposes unless a kernel explicitly requests more) a 6 × 6 tile holds 166 output blocks and one thread block occupies a full SM’s shared memory, capping occupancy near 12 %. Requesting the larger 64 KB or 164 KB per-block limits (dynamic shared memory) raises cap and occupancy proportionally. As Section 6.3 shows, low occupancy is not the binding constraint here: the tile count exceeds the SM count at every level, so the device stays filled at one block per SM, and DRAM/L2 request volume is what the schedule collapses. For output blocks too large to stage (or when higher occupancy is wanted at fixed tile width), an L2-resident accumulator is the alternative tier: A100 L2 atomics are cheap, and the sorted schedule makes the atomic-free register path the default because it avoids the atomic traffic whenever the tile fits.

5.2

MPI orchestration

The parallel product runs the selected on-rank kernel on a merged local operator: off-process rows of P are gathered with a block-granular broadcast over a dedicated PetscSF (one MPI datatype per bsr × bsc block), the on-rank triple product runs in the global coarse block-column space, and the result assembles through the blocked COO path [1]. All gather maps, communication plans, and transpose permutations are cached for values-only reuse. When P is unchanged (hierarchy reuse) only operator values refresh. Both backends implement this native parallel product for the Galerkin triple product and for the prolongator-smoothing product A · P , and both supply the distributed 13

Table 6: Shared-memory tile budget (FP64). cap output blocks per tile at the default 48 KB static shared memory, and at the 64/164 KB dynamic opt-ins, for the block shapes in the elasticity hierarchy. The register accumulator bounds a single output block at 6 × 6. output block

B/block

cap @48 KB

cap @64 KB

cap @164 KB

1 × 1 (scalar) 3×3 3 × 6 (W ) 6 × 6 (Ac )

8 72 144 288

6000 666 333 166

8000 888 444 227

20500 2277 1138 569

subset MatAXPY the smoothing sweep requires (the off-diagonal blocks of the two operands carry different column compressions, which are reconciled through their sorted global column maps on device). No stage of the distributed setup converts the blocked operands to a scalar type, so the operators stay blocked and device-resident across the whole hierarchy, with the off-process exchange the only host-visible step (and only when GPU-aware MPI is unavailable). Communication-avoiding variants of the MPI product are analyzed in [6]. Our focus is the on-rank kernel, and we adopt the standard row-wise decomposition.

5.3

Kokkos backend and performance portability

The blocked variants also express in Kokkos [13]: the two-pass team kernel and the transpose-free product carry over directly, and the search-free single-owner schedule that eliminates discovery work on the host expresses as a Kokkos range-parallel executor. We use the two backends as a controlled portability comparison of the same algorithms: raw-CUDA explicit warp ownership against Kokkos team and range parallelism. The comparison shows that device performance is a property of the control the backend exposes, not of the algorithm alone. On the A100 base case, the search-free schedule that prevails on the host is also the lowest-time CUDA kernel once it carries explicit warp ownership (its tiled form, Section 6.3), but its direct Kokkos expression, the range-parallel executor that wins on the host, is 2.7× slower than the Kokkos team kernels (217 vs. 82 ms on the fine-level product, Table 7). The deficit is structural, not a tuning matter: a range policy gives one thread each output block, so there is no warp-level sharing of the reread operand blocks, and the executor moves 2.7× the device memory traffic of the team kernels (46.2 against 17.1 GB on the fine product), its dominant A · P kernel running at a third of peak DRAM bandwidth at 11 % occupancy. Expressing that sharing in Kokkos requires the team-and-scratch machinery, which is the tiled executor measured below. The team and transpose-free Kokkos kernels instead run the base-case Galerkin product in 82 ms, within a factor of 1.8 of the tiled CUDA kernel (45.5 ms of kernel time, Table 8). The search-free schedule therefore stays a host-only variant in Kokkos, reachable explicitly but never selected automatically on the device. The tiled executor itself makes the same point from the other side. Its Kokkos expression (a TeamPolicy kernel with the accumulator tile in team scratch, one thread per output slot, the same schedule arrays) reproduces the CUDA kernel’s traffic — 13.3 against 10.5 GB on the fine-level product, both below the team kernels’ 17.4 — but not its time: the 48 KB scratch tile admits one team per SM, achieved occupancy falls to 4.6 %, and the kernel runs 168.8 ms against CUDA’s 45.5 and the Kokkos team kernels’ 82.4. The same schedule, the same staging strategy, and byte-level output identity across backends (Section 6.3) still leave a 3.7× time gap, because at one resident team per SM the portable kernel cannot keep enough loads in flight where the raw-CUDA kernel can. The tiled variant is therefore reachable in the Kokkos backend (-mat product algorithm tile) but excluded from its automatic default, which keeps 14

the transpose-free and team kernels. Selection is per product. Any variant can be forced globally through -mat product algorithm, but the default (the option left unset) selects the variant automatically for each MatProduct, and therefore for each multigrid level, since each level’s Galerkin product is its own product built from that level’s operator. The CUDA default is the tiled kernel when the output block fits the per-slot register accumulator (Section 5.1), the fused single-pass path when materializing the intermediate would exhaust free device memory, and the two-pass path otherwise. The feasibility estimate accounts both the W values and the schedule build’s transient, which is linear in the contribution count of the product graph (four index arrays plus the sort’s temporary, about 24 bytes per contribution) and dominates the estimate at small block sizes; the count is one cheap device reduction per schedule. Section 6.5 shows this fallback extending the blocked type’s capacity past the scalar backends at bs=1. The Kokkos default is the transpose-free product when the operator is structurally symmetric (which the aggregation coarsening preserves from the fine operator down the hierarchy) and the two-pass team kernel otherwise, on both host and device. The range-parallel search-free schedule is excluded on the device for the traffic reason above. Each per-level choice is recorded through PetscInfo, and because the default selects the same variant the sweeps identify per level, the automatically selected run reproduces the forced-algorithm results of Section 6 with no per-level regression.

6

Performance

Section 3.5 validated the traffic model’s orderings against the two variants that exist in the blocked backend today. This section reports the surrounding systems results: the cost of building the hierarchy once the setup runs entirely on the device (Section 6.2), the quality of the CUDA numeric kernels relative to their Kokkos expression and the resulting motivation for tiling (Section 6.3), the prolongator-filtering trade (Section 6.4), the scalar bs=1 case (Section 6.5), and the full-GPU finite-element-to-solve pipeline (Section 6.6).

6.1

Experimental setup

The single-device base case is the isotropic elasticity problem of Section 2.3 at 160×16×16 Q2 hexahedral elements (1,045,440 displacement DOF: the 349,569-node Q2 mesh minus the 1,089-node clamped face, times three components), a 3 × 3 fine block structure carrying a six-dimensional rigid-body null space, on one A100 (Table 1). The coarsening builds the four-level hierarchy of Table 2 with 3 × 6 and 6 × 6 coarse blocks. The solver configuration is held fixed across every variant and backend (conjugate gradient with the unpreconditioned residual norm to a relative residual of 10−8 , preconditioned by smoothed aggregation with one Chebyshev point-block-Jacobi smoother per level and a direct coarse solve), and every configuration must reproduce the frozen convergence gate of 19 iterations to a beam-tip displacement of 1.23166 exactly, so that a variant can only change time and memory, never the mathematics. Absolute A100 kernel measurements are Nsight Compute counters (dram bytes, lts t bytes, achieved occupancy, kernel duration). End-to-end product times are GPU-fenced MatPtAPNumeric. The study is single-device by design. Multi-rank runs enter as convergence and parity gates for the distributed implementation, and the weak-scaling study of this problem family is deferred to follow-on work (Section 6.6). The canonical run command (PETSc ex56cu, native-CUDA backend; the Kokkos twin substitutes -dm mat type baijkokkos -dm vec type kokkos) is ./ex56cu -dm_plex_dim 3 -dm_plex_shape zbox -dm_plex_simplex 0 -dm_plex_box_faces 160,16,16 -dm_plex_box_lower 0,-0.5,-0.5 -dm_plex_box_upper 10,0.5,0.5 -lx 10 -petscspace_degree 2 -run_type 4

15

-dm_mat_type baijcuda -dm_vec_type cuda -snes_type ksponly -snes_max_it 1 -n_solves 2 -ksp_type cg -ksp_norm_type unpreconditioned -ksp_rtol 1.e-8 -pc_type gamg -pc_gamg_aggressive_coarsening 1 -pc_gamg_threshold 0.05 -pc_gamg_threshold_scale 0.5 -pc_gamg_coarse_eq_limit 2000 -pc_gamg_process_eq_limit 1000 -pc_gamg_prolongator_filter 0.03 -pc_gamg_prolongator_filter_scale 0.5 -mg_levels_pc_type pbjacobi -mg_coarse_ksp_type preonly -mg_coarse_pc_type bjacobi

with -mat product algorithm naming the product variant under study (Section 5.3 describes the default, per-level selection).

6.2

Baseline: building the hierarchy on the device

The blocked pipeline is designed so that no operator returns to the host between assembly and solve: aggregation, the tentative prolongator, prolongator smoothing, the prolongator filter, and the Galerkin products all run natively on the blocked device operator, on both backends. On the 1.05M-DOF base case of Section 6.1 the one-time hierarchy construction (PCSetUp: coarsening, prolongator construction, the three Galerkin products with their symbolic phases, and the coarse factorization) takes 2.16 s on the native-CUDA backend and 2.30 s on the Kokkos backend, alongside a 0.63 s finite-element Jacobian assembly and a 0.37 s solve (19 iterations). Under hierarchy reuse the recurring setup falls to 0.10–0.15 s, about half of it the numeric Galerkin products (Section 6.4). The MatConvert count is the host-detour indicator for the residency claim: no setup stage converts the blocked operator to a scalar host type. With a device (cuDSS) coarse factorization the setup runs with no host detour at all; a run that instead factors the coarsest level on the host shows seven MatConvert calls totaling 0.02 s, all belonging to that host factorization (4.3 MB of coarse-operator traffic), not blocked-operator detours. displacement 1.23166). Where the numeric time goes. Within the solve, the Galerkin triple product and the finiteelement Jacobian assembly are the two heavy operators. Section 3.5 placed the triple product against the traffic model, and Section 6.6 reports the assembly. The remaining performance question for the triple product is not which algorithm to run (the model and Section 3.5 answer that per level) but how well the chosen algorithm is realized as a GPU kernel, which the next section quantifies.

6.3

Kernel realization: accumulation strategy, team kernels, and V4

The Kokkos and native-CUDA backends express the same algorithms, so comparing their kernels on the fine-level product isolates kernel quality from algorithm choice. The comparison spans three accumulation strategies. The row-parallel kernels of both backends (V0, V1, V4) use the same decomposition: a team per output block row whose lanes split the row’s blocks, sum each block product in registers, and accumulate into the output in global memory with atomic adds. Both backends register-stage the operand blocks ahead of the atomic accumulation; by default they differ only in launch shape, the Kokkos kernels running one 16-lane team per thread block and the native-CUDA kernels packing sixteen such teams into a 256-thread block. Table 7 runs the CUDA kernels at the Kokkos shape so the comparison isolates the backend. The tiled kernel V2b (Section 5.1) stages the output tile itself in shared memory and accumulates each slot exactly once, in schedule order. Table 7 reports the fine-level numeric product for all variants on both backends, product kernels alone. The totals of Table 5 additionally include the transpose value gather and helper kernels, which accounts for the small differences between the two tables. Two facts stand out. First, the algorithm ordering is identical across backends (V4 has the lowest time of the row-parallel

16

variants, V1 the highest), as the model requires. Second, operand staging and launch geometry, not the backend, set the traffic of the row-parallel A · P kernel. An attribution sweep varying operand staging (on or off), sub-warp width (16 or 32 lanes), and teams per block (1 to 16) decomposes the gap. An unstaged launch at the CUDA default geometry moves 99.8 GB, 17× the Kokkos team kernel, at 91.5 % occupancy and 414 GB of L2 request volume. Register-staging the operand blocks ahead of the atomic accumulation — now the default, matching the staging the Kokkos team kernel already carried — cuts this to 11.7 GB (8.5×), and matching the Kokkos launch geometry (one 16-lane team per block row, the geometry of Table 7) removes the remaining factor, to 5.6 GB at 47 % occupancy and 119 GB of L2 — at parity with the Kokkos kernel’s 5.5 GB. The 17× thus factors into staging (8.5×) and launch geometry (2×); neither is a backend property, since Kokkos generates the CUDA machine code. Staging trades occupancy for locality: fewer block rows in flight keep the live output and operand working set L2-resident, so the atomic read–modify–writes and operand re-reads no longer spill to DRAM. At matched geometry the row-parallel backends are therefore at parity on the A · P kernel (5.6 vs 5.5 GB). They part on the full products: the two-pass V0 total and the fused V1 move 1.9–2.2× the Kokkos traffic on CUDA, and the excess sits in the R W reduction (21.8 against 11.6 GB) and V1’s A-block recompute rather than in launch geometry. We have not isolated its cause (the layout of the explicit transpose R and the W re-read pattern are the candidates) and record it as an open observation. Table 7: Fine-level numeric PtAP on both backends (A100, Nsight Compute, DRAM GB / kernel ms). Kokkos generates CUDA kernels on this platform, so to isolate the backend from the launch shape the warp-cooperative CUDA team kernels are run at the Kokkos launch geometry (one 16-lane team per block row, operand staging on both). On CUDA the scheduled kernels V2b (tiled) and V3 (flat) are the lowest-time products, within noise of each other; the two-pass V0 and the fused V1 move more in their R W reduction and recompute passes, and V4 is at parity. On Kokkos the tiled V2b reaches the same low traffic, but its flat form V3 is the host-oriented range executor, which over-fetches on the device (46.2 GB, Section 5.3). Kokkos (team) variant V0 (two-pass) V1 (fused) V4 (mtm) V3 (flat) V2b (tiled)

CUDA (warp-cooperative)

DRAM (GB)

ms

DRAM (GB)

ms

17.1 184.3 15.5 46.2 13.3

81.6 573.8 75.6 217.0 168.7

27.4 407.8 13.3 10.5 10.5

85.8 678.0 74.5 48.1 45.4

V4 is at parity, and the scheduled kernels (V2b tiled, V3 flat) are the lowest-time products, within noise of each other on CUDA. V1 is the negative result on both backends: its traffic is dominated by re-reading the A blocks across the prolongator fan-out, which staging cannot hoist. The scheduled kernels, where the CUDA backend’s effort went, remove the working-set exposure by construction (a sorted schedule and single-touch accumulation). On the A·P kernel itself, V2b moves 5.3 GB in 24.7 ms, below even the Kokkos team kernel (5.5 GB, 46.7 ms) and at half its time. Against the stage-consistent floor (the A · P pass alone must stream BA ZA + BP ZP + BW ZW = 2.81 GB at the measured ZW ) V2b sits at 1.9× and the team kernel at 2.0×; the matched-geometry warpcooperative CUDA kernel also sits at 2.0× (an unstaged default-geometry launch, at 99.8 GB, sits at 36×). On the full product, V2b’s 10.5 GB is 2.5× the 4.28 GB two-pass floor of Table 3. The Kokkos expression of V2b (same schedule, same staging) reaches the same DRAM class (13.3 GB total) but not the same time (168.8 ms, an occupancy limit); Section 5.3 reports that portability 17

result. Profiled, the tiled CUDA kernel is not limited by any bandwidth tier: its busiest unit is the shared-memory/L1 pipe at 60 % of peak, with DRAM at 14 % and L2 under 25 %, the intended effect of staging the accumulators in shared memory, which moves the pressure the team kernel places on L2 (89 % of L2 bandwidth) off the bandwidth tiers entirely. Its optimality evidence is therefore the byte count against the floor rather than a saturation figure, which redundant traffic could inflate. V4 realizes the coarse-level fix. The mtm variant (V4), which parallelizes over fine block rows uniformly across levels, matches V0’s fine-level wall at fewer bytes and runs 7.5× faster than V1 on the Kokkos backend, while eliminating the coarse-level occupancy inversion that Section 3.5 identified for V1: the level-2 product falls from V1’s 138.8 ms to 2.3 ms (60×), because fine-row parallelism keeps the device filled on the small coarse matrices where one-team-per-coarse-row starves it. The warp-cooperative CUDA V4 behaves the same way (level-2 at 1.6 ms, below level-1’s 6.6), so fine-row parallelism delivers the occupancy fix on both backends; what it cannot deliver is the low traffic of the tiled kernel, which is that kernel’s contribution. V2b per level. The tiled kernel’s launch geometry is one thread block per tile of output blockslots rather than per block row: the grid is the tile count ⌈nnzC /cap⌉ (thousands of tiles even on the coarsest operators), so the device stays filled independent of level, and each output block is accumulated once in shared memory from an L2-resident operand working set. Table 8 reports the per-level triple product, tiled CUDA V2b against the portable Kokkos team baseline of Table 7 (the in-session CUDA tile numbers reproduce the earlier committed profile within one percent of session drift). V2b improves on the team baseline at every level, by 1.8× at the fine level and 2.3–2.5× on the coarse levels, and both columns are monotone in problem size — neither kernel family exhibits a coarse-level inversion, the tile because of its slot-based grid, the team baseline because V4’s fine-row parallelism serves its coarse levels. The tile runs at low achieved occupancy (12–25 %: the 48 KB shared-memory tile admits one block per SM) because DRAM/L2 request volume, not occupancy, is the binding resource and the schedule collapses it. Table 8: Per-level fine-to-coarse Galerkin triple product (W +Ac ), portable Kokkos team kernels vs. the tiled CUDA V2b (A100, Nsight Compute, DRAM GB / kernel ms). The tile improves time at every level and both columns are monotone in problem size (no coarse-level occupancy inversion in either kernel family). Kokkos (team) level l0 (fine) l1 l2

V2b (tiled, CUDA)

DRAM (GB)

ms

DRAM (GB)

ms

speedup

17.4 1.41 0.062

82.4 12.3 4.0

10.5 0.81 0.064

45.5 5.4 1.6

1.8× 2.3× 2.5×

The inspector–executor (V3) pays for itself where discovery work is expensive. The inspector–executor schedule (Section 3.3) is implemented in the Kokkos backend as a scheduled form of the two-pass V0: the symbolic phase buckets every block multiply of both products by output slot, in V0’s iteration order, so the executor, a branch-free sweep with no hash, search, atomics, or graph traversal, produces results byte-identical to V0, which we verify bitwise. On the host control (the canonical case of Section 3.5) the schedule repays its construction within the

18

first numeric pass (the measured reuse crossover is N = 1), and in the hot A-only refresh regime the numeric falls by 1.4–1.6× across independent run pairs: the eliminated discovery work of the host path (per-contribution binary search and serialized team loops), with the measured schedule footprint 0.73–0.88× the W value storage at 4 bytes per entry and its stream smaller than the W write-and-read it replaces, consistent with the model’s placement of V3. On the A100 we implement the same schedule as a device executor, staging each tile’s accumulators in shared memory exactly as the SMEM-tiled V2b, so that the schedule encoding is the only difference between the two. It is byte-identical to V0, verified bitwise. Against the staged warp-cooperative two-pass kernels the scheduled forms cut the fine-level numeric 1.9× by the profiler’s kernel meter (85.8 to 45.4 ms, Table 7) and cut the scalar Q2 case of Section 6.5 1.7× by the GPU-fenced wall meter (25.2 to 15.0 ms per product), the discovery- and atomic-elimination gain of a sorted schedule. Against V2b, however, the two are within a few percent on the hot numeric on both cases, a wash. Two measurements explain the parity. First, the 16-bit packing rarely engages at production sizes: the elasticity fine level and the dominant levels of the scalar Q2 case exceed the 216 operand span and fall back to 32-bit (only some coarse-level W schedules pack), so V3’s encoding equals V2b’s eight bytes on the levels that dominate the time. Second, the executor kernels are latency- and occupancy-bound rather than bandwidth-bound (their measured DRAM throughput is 7–42 % of peak, typically 10–18 %), so the entry width is second order regardless. V3’s symbolic phase, meanwhile, costs 1.7× V2b’s on the scalar case (1.89 versus 1.14 s): its extra schedule build (a per-slot minimum/maximum reduced off the device and a host partition of the tiles) comes on top of the contribution sort the two share. On the device V2b is therefore the better choice of the two, the same hot cost at a cheaper build. The inspector–executor’s compression advantage is therefore a property of the host path, where the discovery work it removes is expensive. On the device, where that work is already cheap, the simpler tiled schedule prevails. The executor measurements close the V2c question. V2c (Section 3.3), which stages per-tile W -panels on chip and holds the lowest reachable floor in Table 3, needs no implementation to be judged on this hardware: the tile and executor kernels that would host its staging are latency bound at 7–42 % of DRAM peak, so the further traffic reduction the W -panels would buy is second order. We close the question by measurement and leave V2c as the analytic floor of the model rather than an implementation. The verdict is itself an instance of the division of labor this paper argues for: the traffic model ranks the floors correctly, but once a sorted schedule has collapsed the request volume, the binding resource moves off the bandwidth tiers altogether, and chasing a still-lower floor buys nothing a timer can see. V2c becomes worth implementing only on a machine whose scheduled kernels are bandwidth-bound again, a re-tuning question rather than an open algorithmic one.

6.4

Prolongator filtering study

Filtering (Section 4) trades a controlled perturbation of the smoothed prolongator for reduced work and memory. Table 9 sweeps the threshold on the frozen base case (single A100, blocked Kokkos backend): the iteration count is 19 at every threshold, including the unfiltered product, while the hot MatPtAPNumeric time falls 2.9× (0.230 to 0.079 s at pf = 0.03, 3.7× at 0.05) as the filter thins P and, squared through the Galerkin product, the coarse operators. We have observed that pf = 0.05 is too high for some problems. At pf = 0.03 the filter removes 37 % of the smoothed fine-level P blocks and 68 % on the next level (the threshold halves per level), and the sustained device footprint falls from 10.2 to 9.4 GB; the peak is a threshold-independent transient of the fine-level product’s buffers. The choice inside the band is a mild time optimum, and the hierarchy quality does not move: operator complexity stays low (grid 1.05, operator 1.25 at pf = 0.03). At this problem size the 19

unfiltered product fits the 40 GB device; filtering becomes load-bearing for capacity at larger sizes and at bs=1, where Section 6.5 shows even filtered vendor SpGEMMs exhausting the device. On the eight-GPU elasticity weak-scaling case the same flat band appears (114–115 iterations from 10−2 to 5×10−2 , hot MatPtAPNumeric 0.266 to 0.162 s) and there the unfiltered product does exhaust the device at the base-case size; above the band, aggressive thresholds eventually perturb the operator enough to lose the convergence gate. Table 9: Prolongator-filter threshold sweep on the frozen base case (single A100, blocked Kokkos backend; hot regime, three Galerkin products per solve). Iterations are 19 at every threshold, including unfiltered, while the hot triple-product time falls 2.9× at the paper’s pf = 0.03 as the filter thins the smoothed P . Newton = hot KSPSolve + hot MatPtAPNumeric; memory is the 95th-percentile sampled device footprint. threshold (pf) 0 (off) 0.01 0.02 0.03 0.05

6.5

iterations

MatPtAPNumeric hot (s)

Newton (s)

fine-P nnz

p95 mem (GB)

19 19 19 19 19

0.230 0.112 0.091 0.079 0.061

0.706 0.532 0.501 0.486 0.448

47.4M 39.5M 33.8M 29.9M 24.4M

10.2 9.7 9.5 9.4 9.1

Scalar (bs = 1) results

The scalar case exercises a division of labor. The Galerkin product is native: the same kernel family serves every block size, with bs=1 compile-time instantiations that strip the block loops, and the taxonomy of Section 3 matters more here, not less, because the block contraction degenerates to a single multiply–add and index handling dominates — the hash, search, or atomic cost per 8-byte payload is exactly what the sorted, search-free tile schedule removes. The matrix–vector product is instead dispatched to the vendor : at bs=1 the block-CSR arrays are a scalar CSR matrix, so the mult family routes to the vendor CSR SpMV over the resident device arrays (cuSPARSE directly in the CUDA backend, KokkosSparse::spmv in the Kokkos backend). The benchmark is a 3D Q2 finite-element Poisson problem (PETSc’s snes/tests/ex13) on a 443 -element hex mesh: 893 nodes = 704,969 unknowns at roughly 60 nonzeros per row, a real assembled FEM operator rather than a stencil. Each backend assembles the operator through the same discretization and command line, differing only in the matrix type (-dm mat type baijcuda/baijkokkos/aijcusparse/aijkokkos); a prolongator filter of 0.03 is applied identically, and all four backends converge in 10 iterations with final residuals agreeing to round-off. The hot regime reuses the hierarchy (-pc gamg reuse interpolation) and re-runs the numeric Galerkin products of all six levels each solve. The mesh size is itself a measurement. The intended case was 643 elements (1293 nodes = 2,146,689 unknowns), matching the row count of common scalar benchmarks, but the Q2 operator’s Galerkin fill exhausts the 40 GB device there for every scalar path even with the filter applied: cuSPARSE’s SpGEMM buffers fail from 483 elements up, and the Kokkos Kernels path (which also serves the blocked Kokkos backend’s graphs) fails at 643 . The blocked CUDA backend alone completes the 643 case, because the automatic per-level selection (Section 5) predicts that the fine level’s tile schedule plus W would not fit and falls back to the fused V1 there, running the tiled kernel on the coarser levels — the fused variant’s no-materialization property, analyzed as a traffic loser in Section 3.3, is exactly what makes it the capacity winner. The cross-backend comparison below therefore uses 443 , the largest probed size all four paths complete.

20

Table 10 reports the backend-agnostic measure, the GPU-timed MatPtAPNumeric on the reused hierarchy, the hot path that repeats each Newton step (a per-kernel byte comparison is not meaningful here: the vendor SpGEMMs decompose into hundreds of primitives that also serve the prolongatorsmoothing product, so no kernel-name filter isolates the triple product). The native tile kernel at bs=1 runs the numeric product in 15.0 ms per product, below cuSPARSE’s 17.7 ms and the native Kokkos Kernels SpGEMM’s 16.1 ms of device time, because the sorted tile schedule is searchand atomic-free and rediscovers no structure on reuse. (The benchmark assembles the Jacobian on the host each solve; the GPU-timer bracket excludes the resulting upload wait. The Kokkos Kernels reuse path additionally re-uploads the operator values into its SpGEMM handle each solve, a device–host round trip folded into its 16.1 ms.) Table 10: Scalar (bs=1) PtAP on the Q2 Poisson case (893 nodes, A100, GAMG/CG, prolongator filter 0.03, per-backend assembly, 10 iterations everywhere; Kokkos Kernels built without its cuSPARSE TPL, so its column is the native KK SpGEMM/SpMV rather than a vendor re-dispatch). On the hot, reused-hierarchy numeric Galerkin product the blocked tile kernel is the fastest, below both cuSPARSE and native Kokkos Kernels; the one-time symbolic, built on the device, is within a factor of about three of the scalar paths; and the CUDA backend’s bs=1 SpMV dispatches to the vendor CSR kernel (at cuSPARSE parity), the Kokkos backend to KokkosSparse::spmv (native here). blocked tile blocked Kokkos cuSPARSE Kokkos Kernels (CUDA, bs=1) (mtm, bs=1) (native SpGEMM) hot MatPtAPNumeric / product (ms)

15.0

23.1

17.7

16.1

cold MatPtAPSymbolic (s) hot KSPSolve (ms) SpMV (GFlop/s)

1.17 99 153

0.98 67 128

0.40 102 156

0.47 76 119

The one-time symbolic (the block-graph products and the tile-schedule build) also runs on the device, at 1.17 s against cuSPARSE’s 0.40 and native Kokkos Kernels’ 0.47, the same order as the mature scalar paths. On the dispatched side, the CUDA backend’s bs=1 SpMV dispatches to the vendor CSR kernel and reaches 153 GFlop/s, at parity with cuSPARSE’s 156; the Kokkos backend dispatches to KokkosSparse::spmv, native here (no TPL), at 128 against native Kokkos Kernels’ 119. Either dispatch far exceeds the native block SpMV, whose index handling dominates the degenerate 1 × 1 contraction (61 CUDA, 32 Kokkos). The blocked type is thus a no-cost default at bs=1 on the Galerkin product and the SpMV, and at the largest sizes it is the only path of the four that completes at all.

6.6

End-to-end: the full-GPU blocked pipeline

The blocked types complete the loop from finite-element assembly to solve without leaving the device: the element Jacobian and residual are integrated on the GPU and scattered directly into the blocked operator through the block-granular COO path [1], and the smoothed-aggregation setup and solve then run on that operator with the setup entirely on-device (Section 6.2). Table 11 reports the assembly kernels for both backends: they are near-identical twins at kernel granularity (the Jacobian integration differs by ∼7 % in bytes, a tabulation-staging layout difference, and the COO scatter and residual match within a few percent), which is the portability claim made concrete: the same single-source pointwise physics compiles to equivalent kernels through the two programming models. The -log view transfer counters quantify the recurring (per-Newton, per-solve) phases 21

on the base case: the Jacobian block-COO assembly, every Galerkin numeric, and, with a device (cuDSS) coarse factorization, the entire hot solve run with zero copies in either direction. The only transfer that remains appears when the coarsest level is factored on the host instead of the device: a one-time 1.1 MB factor pull plus a 4.7 KB coarse-level vector each way per iteration. No operator-sized transfer crosses the bus in the recurring phases. Table 11: Finite-element assembly and residual on the base case (A100, Nsight Compute, DRAM GB / kernel ms). The two backends are near-identical twins, evidence that the single-source pointwise physics is performance-portable across the programming models. CUDA kernel Jacobian integration (Ae ) blocked scatter → COO MatSetValuesCOO (fine) residual (volume)

Kokkos

GB

ms

GB

ms

221.5 12.7 4.08 4.09

370.9 24.1 4.3 –

238.7 12.9 4.08 4.14

387.6 24.3 4.2 –

Scope of the measurements. The performance study in this paper is single-device by design: the algorithm space, the traffic model, and the kernel realizations are per-rank questions, and the distributed blocked product (Section 5.2) is implemented, tested for parity against the scalar path at up to four ranks, and exercised by the multi-rank convergence gates. A multi-device weak-scaling study of the same problem family (the 10k × k × k beam at fixed work per device, k = 16, 32, 64 on 1, 8, 64 devices), together with a sweep of the process-reduction control at scale, is deferred to follow-on work.

7

Conclusion and Future Work

This paper treated the rectangular-block Galerkin product as an algorithm space rather than a single kernel. A traffic model that separates a mandatory data floor from the request volume a schedule actually issues predicts, per hierarchy level, which variant has the lower cost as a function of a single geometric ratio, the average number of prolongator block columns per fine block row nP , and the predictions hold on the A100, with a host control confirming that the fused variant’s loss coefficient tracks nP free of GPU launch and occupancy effects (Section 3.5): the fused variant loses by a factor of order nP because it moves the most DRAM traffic despite the lowest floor, and the coarse levels fail on the GPU for occupancy, not traffic, reasons. The model’s reliability is for ratios, orderings, and crossovers rather than absolute byte counts (Section 3.5). The atomic fine-row variant (V4) matches the two-pass reference at the fine level while removing the coarse-level occupancy inversion by construction, and prolongator filtering, with a block null-space-preserving correction that retains the rotational modes a scalar row-sum lumping would drop, both controls operator complexity and is a capacity requirement at scale. The blocked types run the finite-element-to-solve pipeline entirely on the device, with the smoothed-aggregation setup free of host round-trips and the assembly kernels performance-portable twins across CUDA and Kokkos (Section 6.6). The tiled kernel closed the sequential realization gap and, measured latency bound well below the DRAM peak, closed the V2c question: on-chip W -panel staging has no first-order traffic left to save (Section 6.3). The distributed product’s on-rank stages now run the same tiled kernel. At block size one the type matches the vendor scalar paths on the numeric Galerkin product and the matrix–vector product

22

(the latter by dispatching to the vendor CSR kernel over the same arrays), and its fused fallback completes scalar problems whose Galerkin fill exhausts the device for every vendor SpGEMM path (Section 6.5), so the blocked type stands as a general default. The remaining kernel work is narrower still: the block sparse matrix–vector smoother remains latency-limited after reordering, with its tuning space measured as exhausted. Future work includes extending the accuracy-for-cost principle that prolongator filtering exemplifies (Section 4) to reduced-precision arithmetic: reduced-precision Galerkin products, where the tensor cores’ reduced-precision paths (FP16 or TF32 inputs with FP32 accumulation, 156–312 TF/s on the A100) retain a large advantage even at the poor fragment utilization that rules them out at FP64 (Section 1), and, further, coarse levels stored and applied entirely in reduced precision, which the algebraic error structure of multigrid makes particularly tolerable [8]. Other directions are aggregate-aligned additive-Schwarz smoothers with blocks spanning process boundaries, a dedicated study of element anisotropy and aggregation quality, communication-avoiding MPI formulations of the blocked product [6], the deferred weak-scaling study of the 10k × k × k beam on 1, 8, and 64 devices with the process-reduction sweep at 64 devices, and re-tuning the tile models for H100-class parts (larger L2 and shared memory) [9].

Code and Data Availability The implementation is developed in a fork of PETSc. A curated public release accompanies this paper at https://gitlab.com/markadams4/petsc-device-fem-amg: a single commit on an upstream PETSc base containing the complete source, the test problems and run configurations, the measurement data behind the figures and tables, and documentation for building and reproducing the experiments, so that the contribution is viewable as one diff against upstream PETSc. A snapshot of the release will accompany the final version of this paper with an archival DOI.

Use of AI Portions of this work, including code development and debugging, performance-data reduction and analysis, and preparation of this manuscript, were carried out with the assistance of Anthropic’s Claude (Claude Code). All design decisions, experiments, and results were directed, generated, and verified by the author, who takes full responsibility for the content.

Acknowledgments Thanks to the PETSc team for developing a well-engineered numerical library that provided an ideal basis for the extensions developed in this project. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) Program through the FASTMath Institute, under contract number DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory.

References [1] Mark F. Adams. A natively blocked, device-resident algebraic multigrid GPU path in PETSc, 2026.

23

[2] Mark F. Adams, Ravi Samtaney, and Achi Brandt. Toward textbook multigrid efficiency for fully implicit resistive magnetohydrodynamics. J. Comput. Phys., 229(16):6208–6219, 2010. [3] Mark F. Adams, Peng Wang, Jacob Merson, Kevin Huck, and Matthew G. Knepley. A performance portable, fully implicit landau collision operator with batched linear solvers. SIAM Journal on Scientific Computing, 47(2):B360–B381, 2025. [4] Advanced Micro Devices, Inc. rocSPARSE library, 2024. GEBSR (general block sparse row) format supports independent row/column block sizes. [5] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Vaclav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.22, Argonne National Laboratory, 2024. [6] Grey Ballard, Christopher Siefert, and Jonathan Hu. Reducing communication costs for sparse matrix multiplication within algebraic multigrid. SIAM J. Sci. Comput., 38(3):C203–C231, 2016. [7] Nathan Bell, Steven Dalton, and Luke N. Olson. Exposing fine-grained parallelism in algebraic multigrid methods. SIAM J. Sci. Comput., 34(4):C123–C152, 2012. [8] Stephen F. McCormick, Joseph Benzaken, and Rasmus Tamstorf. Algebraic error analysis for mixed-precision multigrid solvers. SIAM J. Sci. Comput., 43(5):S392–S419, 2021. [9] NVIDIA Corporation. NVIDIA A100 tensor core GPU architecture, 2020. Whitepaper; FP64 tensor-core (DMMA) peak 19.5 TF/s vs. 9.7 TF/s FMA, 40 MB L2. [10] NVIDIA Corporation. cuSPARSE library, 2024. Block compressed sparse row (BSR) format supports equal-size blocks only. [11] Luke N. Olson, Jacob B. Schroder, and Raymond S. Tuminaro. A general interpolation strategy for algebraic multigrid using energy minimization. SIAM J. Sci. Comput., 33(2):966–991, 2011. [12] Sivasankaran Rajamanickam, Mehmet Deveci, Christian Trott, Seher Kim, Nathan Ellingwood, Seyong Deveci, Mauro Perego, and Dan Sunderland. KokkosKernels: Performance portable sparse/dense linear algebra and graph kernels. In Proc. IEEE High Performance Extreme Computing Conf. (HPEC), 2021. [13] Christian R. Trott, Damien Lebrun-Grandié, Daniel Arndt, Jan Ciesko, Vinh Dang, Nathan Ellingwood, Rahulkumar Gayatri, Evan Harvey, Daisy S. Hollman, Daniel Ibanez, Nevin Liber, Jonathan Madsen, Jeff Miles, David Poliakoff, Amy Powell, Sivasankaran Rajamanickam, Mikael Simberg, Dan Sunderland, Bruno Turcksin, and Jeremiah Wilke. Kokkos 3: Programming model extensions for the exascale era. IEEE Trans. Parallel Distrib. Syst., 33(4):805–817, 2022. [14] U. Trottenberg, C. W. Oosterlee, and A. Schüller. Multigrid. Academic Press, London, 2001.

24

[15] Ray S. Tuminaro and Charles Tong. Parallel smoothed aggregation multigrid: Aggregation strategies on massively parallel machines. In Proceedings of the 2000 ACM/IEEE Conference on Supercomputing (SC ’00), Dallas, TX, 2000. IEEE Computer Society. [16] Petr Vaněk, Jan Mandel, and Marian Brezina. Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems. Computing, 56(3):179–196, 1996.

25

Related documents

• Shortcutting the Fix: Identifying and Categorizing Agentic Exploits in Software Engineering Benchmarks#668134
• CAROL: Context-Aware Online Learning for Fuzzer Scheduling#668135
• An Autonomy Aware Metamodel for Human AI Collaboration in Software Engineering#668136
• A Unified Policy Architecture (UPA): The Governance Kernel for Enterprise AI Operating Systems#668137
• From Reading Code to Reading Spec: A Verified Layer for LLM-Driven Codebase Maintenance#668138