Matrix-Free 3D SIMP Topology Optimization with Fused Gather–GEMM–Scatter Kernels Shaoliang Yang1 , Jun Wang∗1 , and Yunsheng Wang1
arXiv:2604.18020v1 [cs.CE] 20 Apr 2026
1
Department of Mechanical Engineering, Santa Clara University, Santa Clara, CA 95053, USA
Abstract The matrix-free gather–batched-GEMM–scatter pattern eliminates global stiffness assembly for three-dimensional SIMP topology optimization, but the conventional three-stage implementation forces avoidable DRAM traffic between stages. We present a single fused CUDA kernel, implemented through CuPy’s runtime compilation interface, that performs gather, per-element stiffness multiplication, and scatter accumulation in one pass. On a single RTX 4090 (24 GB), the fused path reaches a problem-size-dependent 4.6–7.3x end-toend SIMP wall-time speedup across 216k–4.9M cantilever elements and 4.4x on the 499,125element torsion benchmark. Against the same-precision FP32 three-stage baseline, the fused path still yields 2.3–4.6x on cantilever and 2.8x on torsion. Isolated CUDA-event cantileveroperator measurements reach 8.9–13.8x per matvec call, while separate instrumented boardpower traces at 216k and 1M show 3.2–4.9x lower energy than matched FP64 runs. A separate bridge stress test shows the same FP32-versus-FP64 three-stage trend under one distributed-load case; direct fused-kernel bridge benchmarks are not reported. We also evaluate a BF16 WMMA variant: a separate PyTorch BF16 GEMM proxy on matching tensor shapes yields 14.3x, but direct condition-number estimates of 6.1e5–2.3e6 across 64k–512k uniform-density test states imply BF16 conditioning products of 2.4e3–9.1e3, far above the 256 threshold, observed alongside BF16 iterative-refinement stagnation at the two tested inner tolerances.
Keywords: topology optimization; SIMP; matrix-free FEM; GPU kernel fusion; BF16 tensor cores; mixed precision; iterative refinement; conjugate gradient; CuPy
1
Introduction
Topology optimization (TO) is the canonical computational design method for lightweight structural engineering: given a design domain, boundary conditions, and a target volume fraction, it distributes material to minimize compliance (or other objective functions) subject to equilibrium and volume constraints. The Solid Isotropic Material with Penalization (SIMP) method [1–3] has become the standard formulation for three-dimensional structural TO because it admits a simple alternating update loop—iterate between a finite-element solve, a sensitivity computation, and an Optimality Criteria (OC) density update—that scales to the million-element range when the linear system is solved efficiently [4, 5]. The GPU is the natural accelerator for this loop: the finite-element assembly and matvec operations are embarrassingly parallel over elements, the sensitivity computation and OC update are elementwise, and the iterative Conjugate Gradient (CG) solver maps to the massively-parallel vector-operation model native to CUDA. ∗
Corresponding author. E-mail: [email protected]
1
Prior work on GPU-accelerated TO has demonstrated the viability of this approach at increasing scales [6–12], but has been constrained by the memory cost of explicit sparse stiffness matrix assembly. For trilinear hexahedral (Q1) elements, the global stiffness matrix K has 576 nelem element-level contributions, so explicit sparse assembly enters a low-million-element VRAM wall once values, indices, and temporary assembly buffers are included on a standard RTX 4090 GPU. The matrix-free baseline studied here addresses this VRAM wall by replacing explicit assembly with a matrix-free Kv operator that evaluates the product element-by-element without forming K globally, enabling the single-GPU Python/CuPy workflow used throughout this paper and avoiding the assembled-CSR memory footprint that would otherwise dominate at the million-element scale [13]. This direction is consistent with recent matrix-free TO implementations on both GPU and CPU platforms [12, 14]. That implementation decomposes the matvec into three sequential kernel launches: a gather that reads element DOF values from the global vector, a batched GEMM that applies the scaled element stiffness matrix, and a scatter that atomically accumulates element contributions back to the global output vector. While this three-stage design is straightforward to implement, it wastes memory bandwidth: the intermediate per-element arrays produced by gather and consumed by scatter must transit DRAM between kernel launches, adding two full DRAM round-trips per matvec call. The kernel-fusion opportunity. The roofline model [15] identifies memory bandwidth as the binding constraint for the gather–GEMM–scatter pattern: the arithmetic intensity of the trilinear hexahedral element matvec is approximately 3–6 FLOP/B, well below the FP32 ridge point of 82 FLOP/B on the RTX 4090. In this regime, the dominant cost is not floating-point operations but DRAM round-trips, and kernel fusion—combining the three stages into a single kernel launch that keeps intermediate data in registers and shared memory—directly reduces matrix (24 × 24 = 576 floating-point values, occupying 2.25 kB these round-trips. The Kunit e in shared memory) is identical for all elements of a uniform mesh and can be broadcast to shared memory once per thread block, so all 24 elements of the per-element matvec read Kunit e from shared memory rather than global DRAM. The expected throughput gain from fusion is approximately 2× in the element-data DRAM path (about 1.8× once index and density reads are included), and the realized end-to-end speedup can exceed this factor, consistent with the additional savings from collapsing three launches into one. Tensor cores as the next frontier. NVIDIA’s Ada Lovelace architecture (RTX 4090) exposes BF16 tensor cores via the WMMA (Warp-Level Matrix Multiply–Accumulate) API, offering 165.2 TFLOP/s in BF16 versus 82.6 TFLOP/s in FP32—a factor-of-2 hardware throughput advantage [16]. Tensor cores have also been examined in prior scientific- computing and mixed-precision solver studies, clarifying both their performance potential and their nontrivial numerical behavior outside deep learning [17–19]. For the 24 × 24 per-element dense GEMM at the heart of the fused kernel, BF16 WMMA offers the prospect of large acceleration at the GEMM stage relative to FP64 batched DGEMM. However, BF16’s 7-bit mantissa (unit roundoff εBF16 = 2−8 ≈ 3.9 × 10−3 under the Carson–Higham convention) raises the fundamental question of whether this precision is sufficient for the Conjugate Gradient solver applied to SIMP stiffness matrices, whose condition numbers can readily exceed the BF16 stability threshold once SIMP penalization and density contrast become large. The classical iterative refinement (IR) framework and later mixed-precision solver literature [20–22] give the standard sufficient condition εBF16 · κ(K) < 1, i.e., κ < 256 for BF16. We use this as a conservative reference threshold for the experiments below. To our knowledge, in the literature reviewed for this paper we did not identify a prior study that directly tests this theoretical prediction in the density-based TO setting considered here. Likewise, in that reviewed literature we are not aware of a prior study that combines BF16 tensor-core kernels with density-based 3D SIMP TO in this single-GPU setting. Contributions. This paper makes four contributions: 1. Fused CUDA kernel. We design and implement a single fused gather–GEMM–scatter 2
CUDA kernel (Algorithm 1) that reduces effective off-chip traffic and achieves 6.0–6.6× per-matvec speedup over the FP64 three-stage baseline in the synthetic hot-path microbenchmark; isolated CUDA-event measurements on the actual cantilever operator (Table 3) yield 8.9–13.8×, consistent with the additional Python/CuPy dispatch overhead avoided by replacing three separate API calls with one launch. The kernel is deployed through CuPy’s runtime kernel-compilation interface with no standalone CUDA build. 2. End-to-end SIMP speedup across load cases. The fused kernel achieves a problemsize-dependent 4.6–7.3× end-to-end SIMP-120 speedup on the cantilever benchmark from 216 k to 4.9 M elements, together with 4.4× on the 499,125-element torsion benchmark. Against the same-precision FP32 three-stage baseline, the cantilever scaling rows yield 2.3–4.6× and the torsion row 2.8×. The bridge hard-problem stress test shows the same FP32 three-stage gain under distributed loading, while the accompanying MBB hardproblem rows remain cap-limited under the current Jacobi-PCG configuration and are not promoted as a main quantitative claim. 3. BF16 WMMA throughput probe. We implement a CuPy-accessible BF16 WMMA matvec kernel for density-based topology optimization and evaluate it through full-kernel per-matvec timings plus a separate 14.3× BF16-versus-FP64 PyTorch GEMM proxy timing on matching tensor shapes and Ada tensor cores; we do not claim an end-to-end BF16 SIMP speedup. 4. Mixed-precision convergence analysis. We analyze the standard precision barrier εBF16 ·κ(K) ≳ 1 for BF16-in-CG on SIMP stiffness matrices, and the reported experiments provide an empirical test of this BF16 barrier in the present density-based 3D SIMP setting. The BF16 path stagnates under the tested Jacobi-preconditioned CG solve, suggesting multigrid preconditioning as the practical path to exploiting tensor-core throughput. Paper organization. Section 2 reviews related work on GPU TO, matrix-free FEM, kernel fusion, and mixed-precision iterative solvers. Section 3 presents the 3D SIMP formulation, the fused kernel design, the BF16 WMMA variant, a roofline characterization of the reported kernels, and the theoretical convergence analysis. Section 4 reports per-matvec profiling, end-to-end SIMP scaling, bridge and torsion stress tests, BF16 convergence experiments, and qualitative topology validation. Section 5 interprets the speedup trends, provides contextual comparison to prior implementations, discusses limitations, and outlines the solver requirements that motivate a separate geometric-multigrid follow-on. Section 6 concludes.
2
Related Work
2.1
Large-Scale GPU Topology Optimization
The memory and computational cost of three-dimensional SIMP topology optimization are dominated by the linear solve at each design iteration, motivating the extensive subsequent literature on solver architectures for this application [4, 23, 24]. The compact 99-line and 88-line MATLAB implementations of Sigmund and Andreassen et al. [25, 26] codified the pedagogical density-filtered SIMP workflow in 2D, while Liu and Tovar [24] extended the sparse-assembly approach to 3D and established a widely used 3D MATLAB baseline, while Aage, Andreassen, and Lazarov [4] scaled 3D SIMP to hundreds of millions of elements by coupling PETSc sparse solvers to FGMRES with algebraic multigrid (AMG)—a purely MPI-parallel approach requiring distributed-memory HPC infrastructure. The extreme limit of this paradigm was demonstrated by Aage et al. [5], who optimized a 1.1-billion finite-element aircraft-wing structure on a PRACE supercomputer. GPU acceleration of topology optimization emerged from 2D implementations (Wadbro and Berggren [6], Challis et al. [7]) and grew to multi-GPU distributed systems (MartínezFrutos and Herrero-Pérez [8]; Herrero-Pérez and Martínez Castejón [9]). At the single-GPU 3
scale most directly relevant to the present work, Hou et al. [10] employed CuPy vectorized sparse matrix–vector products (SpMV) to accelerate 2D and 3D TO problems, reporting sizes up to 63 million elements and demonstrating that high-level CuPy implementations can remain effective at very large scales. That scale comparison should be read cautiously against the present paper: our reported ceiling is set by end-to-end single-GPU SIMP runs with the current Jacobi-preconditioned matrix-free pipeline on a 24 GB RTX 4090, rather than by an isolated sparse-operator throughput study. The present paper differs in focusing on a matrixfree gather–GEMM–scatter operator for trilinear hexahedral 3D SIMP systems, with a fused single-kernel implementation and a mixed-precision conditioning analysis. Qi et al. [11] reported another recent single-GPU 3D solver for topology optimization of continuous fiber-reinforced composites, reaching 67.1 million Wilson’s incompatible elements and 201.3 million design variables on a Tesla V100 with a multigrid-preconditioned conjugate-gradient solver and Taylorapproximation- based element stiffness updates. This result is not directly comparable to ours because it targets a different material class and solver architecture. The matrix-free CuPy baseline used in the present implementation replaces explicit CSR assembly with a matrix-free Kv operator, eliminating the assembled-storage memory wall that otherwise dominates low-millionelement 3D problems on a single RTX 4090. The present paper advances that baseline by fusing the multi-kernel Python path into a single CUDA kernel and extending it with BF16 tensor-core arithmetic. The closest related work is that of Träff et al. [12], who present two matrix-free 3D SIMP implementations—one in Futhark and one in OpenMP-C—together with an asymptotically exact complexity analysis of the gather–contract–scatter operator. Their solver achieves 65.5 M elements on an NVIDIA A100 (80 GB HBM2e) with SSOR V-cycle multigrid preconditioning. The present paper shares the same mathematical operator but contributes a CuPy runtimecompiled realization of the fused CUDA path, a BF16 tensor-core variant, and a focused mixedprecision failure analysis in the TO linear-solve context. Wang et al. [14] recently demonstrated a matrix-free MATLAB implementation reaching 128 M elements on a CPU workstation (64 GB RAM) using geometric multigrid with non-dyadic Galerkin coarsening, indicating that matrixfree methods are an increasingly important approach for large-scale TO across CPU and GPU platforms.
2.2
Matrix-Free FEM Operators on GPUs
The principle of evaluating Kv products element-by-element without assembling K globally was introduced for structural mechanics by Hughes, Levit, and Winget [27] and extended to nonlinear schemes by Carey and Jiang [28]. Schmidt and Schulz [29] translated this elementby-element CG to GPU in 2,589 lines of CUDA C++, establishing viability for 3D TO at the hundred-thousand element scale. Wu, Dick, and Westermann [30] subsequently reported a highresolution GPU topology-optimization system that demonstrated the practical importance of GPU-resident solver pipelines for large structural problems. In the high-order finite-element community, Kronbichler and Kormann [31, 32] established the matrix-free tensor-product evaluation paradigm in deal.II, showing that sum factorization makes high-order operator application efficient on hexahedral elements. The CEED/MFEM project [33, 34] has since ported these matrix-free kernels to GPU. Cao et al. [35] separately reported 85–100% of the A100 roofline for matrix-free high-order stencil FEM via on-the-fly geometric-factor recomputation. These high-order results do not transfer directly to the trilinear hexahedral (k = 1) elements used in SIMP TO, where sum factorization offers no advantage and where the dominant cost is the 24-DOF per-element dense multiply—the regime the present paper targets.
4
2.3
Kernel Fusion for Memory-Bound GPU Kernels
GPU kernel fusion—combining multiple sequential kernel launches into a single launch that keeps intermediate data in registers or shared memory—is a well-studied technique for memorybound workloads. Filipovič et al. [36] demonstrated up to 2.6× speedup over cuBLAS by fusing BLAS-1 and BLAS-2 chains using a source-to-source compiler that analyzes kernel dependency graphs. Wahib and Maruyama [37, 38] built an end-to-end framework for fusing and transforming stencil kernels in production CFD and weather codes, reporting 1.2–1.75× speedup by eliminating inter-kernel DRAM traffic. The Halide language [39] decouples algorithm from schedule and provides a formal framework for the fusion–recomputation trade-off that underlies these approaches. The present fused kernel differs from the reviewed fusion literature in an important respect: the intermediate data that is fused is neither a BLAS-level chain nor a structured Cartesian stencil, but a three-stage pipeline in which (i) the gather stage accesses a globally irregular element-to-DOF table, (ii) the GEMM stage applies a dense 24×24 symmetric matrix scaled by a per-element scalar, and (iii) the scatter stage performs atomic-add reduction back to a globally irregular DOF array. In the fusion systems literature reviewed for this paper, we are not aware of an automated framework that handles this combination; the fused kernel is therefore hand-written in CUDA C and inlined in Python through CuPy’s runtime kernelcompilation interface. The present fused implementation builds directly on the unfused threestage matrix-free operator used as the FP64/FP32 baseline in this paper. Träff et al.’s Futhark implementation implicitly fuses a closely related gather–contract–scatter pattern at compile time; the present work makes that fusion explicit in CUDA C via CuPy and extends it with the BF16 WMMA path and the mixed-precision convergence analysis.
2.4
Mixed Precision and Tensor Cores in Iterative Solvers
The theoretical foundation for mixed-precision iterative refinement was established by Carson and Higham [20], who derived sufficient conditions for convergence and error bounds for threeprecision IR in terms of the working-precision unit roundoff and a problem-dependent constant. Higham and Mary [21] provide a comprehensive survey of mixed-precision algorithms in numerical linear algebra, cataloguing the regimes in which FP16 and BF16 arithmetic can safely accelerate factorization and iterative refinement. The five-precision GMRES-IR framework of Amestoy et al. [40] is the most general recent extension, achieving high-precision solutions via structured cascade of low-precision inner solves. On the hardware side, Markidis et al. [17] characterized the programmability, performance, and precision of NVIDIA tensor cores for non-ML workloads immediately after their introduction; Fasi et al. [18] subsequently provided a systematic empirical study of rounding modes, subnormal handling, and accumulation order for V100, T4, and A100 tensor cores, showing that hardware behavior deviates from IEEE 754 in subtle ways relevant to scientific computing. Henry, Tang, and Heinecke [22] specifically advocated BF16 for HPC iterative solvers, projecting convergence over “a large range of condition numbers” when an FP32 outer refinement loop is used; Kalamkar et al. [41] clarified BF16 semantics as a 16-bit format with a 7-bit mantissa and FP32 dynamic range. The most directly comparable mixed-precision solver application is the work of Haidar, Tomov, Dongarra, and Higham [19], who demonstrated up to 4× speedup over FP64 on dense linear systems using V100 FP16 tensor cores with FP32 accumulation in a MAGMA-based IR scheme. Clark et al. [42] showed that FP16/FP32 mixed-precision CG achieves production accuracy in lattice QCD using mixed-precision Krylov solvers with reliable updates, illustrating that low-precision bulk arithmetic can succeed when the solver strategy and operator class are favorable. McCormick, Benzaken, and Tamstorf [43] extended this to a rigorous algebraic framework for mixed-precision multigrid solvers on symmetric positive-definite elliptic systems. More re5
cently, Bai, Ootomo, and Yokota [44] demonstrated a tile-grained mixed-precision single-kernel CG solver on GPUs, further underscoring the relevance of kernel-integrated mixed-precision Krylov designs. The present paper contributes a focused BF16 mixed-precision analysis in the topology-optimization context. The reported study includes a separate 14.3× BF16-versusFP64 GEMM proxy timing at 512 k elements, but the classical iterative refinement scheme fails to converge to an acceptable solution state—stagnating at large compliance error—because direct power-iteration estimates (Section 4.8) indicate κ(K) ≈ 6.1 × 105 –2.3 × 106 across the tested mesh sizes (64 k–512 k elements), placing εBF16 · κ ≈ 2.4 × 103 –9.1 × 103 —well above the 1/εBF16 = 256 sufficient threshold implied by the standard IR bound εBF16 κ(K) < 1 from Carson and Higham [20]. Henry et al.’s optimistic BF16-IR projection [22] is thereby qualified for the specific tested engineering PDE class, and a path forward—using BF16 as a multigrid smoother where the coarse-grid spectrum is bounded—is identified.
2.5
Tensor Cores for PDE and FEM Computations
Beyond dense linear algebra, tensor cores have been applied to scientific computations through a series of operator-mapping innovations. Dakkak et al. [45] showed that arbitrary reductions and prefix scans can be reformulated as WMMA operations, achieving up to 100× speedup over stateof-the-art reduction kernels. Ootomo and Yokota [46] developed a split-FP16 compensation scheme that recovers single-precision accuracy from FP32-accumulate FP16 WMMA on V100 tensor cores, while Ootomo, Ozaki, and Yokota [47] extended this to FP64 emulation via INT8 tensor cores. Stencil computations on tensor cores—motivated by iterative PDE solvers—have been pursued via im2col reshaping in ConvStencil [48] and low-rank matrix approximation in LoRAStencil [49]. Recent tensor-core work on finite-element-style operators has also considered tensor-product kernels in scientific computing. Cui [50] studies acceleration of tensor-product operations with tensor cores, illustrating that scientific-operator mappings can benefit from specialized tensorcore layouts even outside machine learning. The present paper differs in both operator shape and problem setting: it targets the 24 × 24 trilinear-hexahedral SIMP matvec in density-based topology optimization, uses BF16 WMMA fragments, and analyzes the resulting BF16-in-CG convergence barrier. In the literature and public implementations and supplements reviewed for this paper, we are not aware of a published report matching this exact tensor-core operator mapping for density-based topology optimization on a single GPU.
2.6
Preconditioner Gap and the Path to Future Work
The CG iteration count per SIMP step is determined by the preconditioner and directly controls end-to-end wall time. The Jacobi (diagonal) preconditioner is the simplest choice—trivially p parallelizable and free to compute—but requires O( κ(K)) CG iterations for convergence. For the cantilever problem, Jacobi yields a few hundred iterations per SIMP step in the present work; for more challenging BVPs (torsion, column) it can approach the 1,000-iteration cap; for problems with near-rigid-body modes (3D MBB with the standard pin boundary condition) it fails entirely. Amir, Aage, and Lazarov [51] demonstrated that multigrid-preconditioned CG can substantially reduce per-step counts for 3D SIMP problems, with the hierarchy reused across SIMP iterations for further efficiency. Träff et al. [12] use SSOR smoothing within a Vcycle preconditioner and report large-scale single-GPU runs up to 65.5 million elements. Peetz and Elbanna [52] systematically compared AMG and geometric multigrid (GMG) for 3D TO, highlighting the trade-off between per-iteration efficiency (GMG advantage) and robustness to topology evolution (AMG advantage). McCormick et al. [43] provide a rigorous algebraic framework for mixed-precision multigrid on elliptic problems, which is directly applicable to the TO stiffness operator when the multigrid hierarchy is fixed. These references collectively define
6
the open problem that most directly limits the present solver’s applicability and motivate the geometric multigrid extension proposed as future work in Section 5.
3
Methodology
3.1
Three-Dimensional SIMP Topology Optimization
We consider the classical minimum-compliance topology optimization problem over a fixed hexahedral mesh of nelem trilinear eight-node (Q1) elements: min f ⊤ u(ρ) ρ
s.t.
K(ρ) u = f ,
1
nX elem
nelem e=1
ρe = Vf ,
0 < ρmin ≤ ρe ≤ 1,
(1)
where f ∈ Rndof is the external load vector, u ∈ Rndof the displacement vector, Vf the prescribed volume fraction, and ρe ∈ [ρmin , 1] the design density of element e. The global stiffness matrix is assembled from element contributions via the SIMP (Solid Isotropic Material with Penalization) material interpolation [1]: K(ρ) =
nX elem
p unit B⊤ e ρmin + (1 − ρmin ) ρe Ke Be ,
(2)
e=1
where Be ∈ {0, 1}24×ndof is the Boolean gather/scatter matrix that selects the 24 local DOFs of element e from the global vector, Kunit ∈ R24×24 is the element stiffness matrix for a unite modulus isotropic material, p is the penalty exponent (typically p = 3), and ρmin = 10−9 prevents singularity in void regions. For the three-field regularization, raw densities ρe are smoothed through a cone filter of radius rmin to obtain ρ̄e , following the standard density-filter/projection pipeline of topology optimization studies [26, 53, 54], and then projected through a smoothed Heaviside: tanh(β η) + tanh(β(ρ̄e − η)) ρ̃e = , (3) tanh(β η) + tanh(β(1 − η)) with projection threshold η = 0.5 and β controlled by a deterministic continuation schedule in the quantitative benchmark runs. Design sensitivities ∂c/∂ρe are computed analytically and used by the Optimality Criteria (OC) update rule [3] with bisection to enforce the volume constraint at each SIMP iteration. In the quantitative benchmark runs reported later, these nominal parameters are embedded in a deterministic continuation schedule that ramps p, β, rmin , and the OC move limit; the exact schedule is stated in Section 4.1. For reporting and validity checks, the grayness metric is defined as g=
4
nX elem
nelem e=1
ρe (1 − ρe ),
so g = 0 denotes a fully binary design and larger values indicate more intermediate-density material.
3.2
Matrix-Free Kv Operator: Three-Stage Baseline
Explicit CSR assembly of K is infeasible at the million-element scale on a single GPU: for a trilinear hexahedral mesh each element contributes 24 × 24 entries indexed by 24 degrees of freedom (DOFs), so the CSR index arrays alone occupy O(576 nelem ) integers—exceeding 24 GB vram once explicit sparse assembly reaches the low-million-element regime on an RTX 4090. Instead, the matrix–vector product w ← Kv is evaluated element-by-element without ever forming K globally: w=
nX elem
p unit B⊤ e ρmin + (1 − ρmin ) ρe Ke Be v,
e=1
7
(4)
Algorithm 1 Fused Gather–GEMM–Scatter kernel (per thread) Require: edof[nelem , 24], v[ndof ], Kunit e [24, 24], ρ[nelem ], p, ρmin Ensure: w[ndof ] (atomic-accumulated) 1: e ← (block index) · (block size) + local thread index 2: if e ≥ nelem then 3: return 4: end if 5: Shared memory: cooperatively load Kunit e [24, 24] once per block p 6: ke ← ρmin + (1 − ρmin ) ρe 7: Gather local DOFs: uj ← v[edof[e, j]] for j = 0, . . . , 23 8: for i = 0, . . . , 23 do 9: fi ← 0 10: for j = 0, . . . , 23 do 11: fi += ke · Kunit e [i, j] · uj 12: end for 13: Atomically accumulate fi into w[edof[e, i]] 14: end for
▷ SIMP scaling
where Be is the Boolean gather/scatter matrix from Equation (2). The baseline Python/CuPy implementation decomposes this into three sequential kernel launches: Stage 1. Gather: For each element e, copy ue = u[edof e ] ∈ R24 from the global displacement ×24 . vector using the precomputed DOF index table edof ∈ Znelem unit p Stage 2. Batched GEMM: Compute fe = ρmin + (1 − ρmin ) ρe Ke ue using CuPy’s batched matrix-multiplication path, which dispatches to cuBLAS over all nelem elements simultaneously. Stage 3. Reduction scatter: Accumulate w[edof e ] += fe with CuPy’s histogram-style reduction back to the global output vector. This three-stage decomposition requires three full DRAM round-trips per matrix–vector product: Gather reads 24 nelem doubles, batched GEMM reads the element vectors again and writes results, and Scatter reads those results and writes back to the global vector. The intermediate per-element arrays—uelem and felem —each occupy 24 nelem × 8 bytes; at nelem = 2 × 106 these two buffers total 768 MB, a non-trivial fraction of the 24 GB vram budget. More critically, these buffers are written by one kernel and read by the next, forcing them to transit off-chip DRAM. For the tested RTX 4090, GDDR6X bandwidth peaks at 1.008 TB/s; for reference, the 80 GB A100’s HBM2e bandwidth is up to 1.94 TB/s. This off-chip traffic is what makes the baseline memory-bound at all sizes tested.
3.3
Fused Gather–GEMM–Scatter Kernel
The fused kernel eliminates the intermediate arrays by performing all three stages within a single CUDA kernel launch, keeping ue and fe in thread registers and shared memory throughout. The implemented FP32 kernel uses one thread per finite element and 128 threads per block. The launch configuration is therefore (⌈nelem /128⌉, 1, 1) blocks × (128, 1, 1) threads, with each thread gathering one element, applying the 24 × 24 dense multiply, and atomically scattering its 24 contributions. Shared memory layout. The unit element stiffness Kunit ∈ R24×24 is the same for all e elements (it depends only on the element geometry and material constants, not on ρe ); it occupies 242 × 4 bytes = 2.25 kB in shared memory and is loaded once per thread block via a cooperative broadcast before the main loop. Threads then read rows of Kunit from shared e memory during the inner GEMM loop, achieving a shared-memory bandwidth amplification of 24× versus a naive global-memory approach.
8
Figure 1: Comparison of the three-stage baseline pipeline (left) and the fused single-kernel implementation (right). The baseline requires three separate CUDA kernel launches per matrix– vector product, each forcing the intermediate per-element arrays uelem and felem to transit DRAM. The fused kernel keeps these in thread registers and shared memory, substantially reducing effective off-chip traffic (theoretical element-data reduction ≈ 2×). Memory traffic analysis. Let Belem denote bytes per element per Kv call. The baseline three-stage path transfers: Bbaseline =
|24·8 {z }
gather out
+
|24·8 {z }
GEMM in (elem)
+
|24·8 {z }
GEMM out (elem)
+ 24·8 | {z } = 768 bytes/element. (5) scatter in
The fused kernel eliminates the two intermediate passes; the minimum traffic is: Bfused =
24·8 | {z }
gather global read
+
|24·8 {z }
scatter global write
9
= 384 bytes/element,
(6)
a factor of 2× reduction in element-data DRAM traffic. Including one read of the edof table (96 bytes/element) and the ρ array (4 bytes/element), the theoretical minimum DRAM reads rise to ≈ 480 bytes/element for the fused kernel versus ≈ 868 bytes/element for the baseline, still a factor of 1.8×. In practice, L2 cache captures repeated edof accesses and the runtime speedups observed in Section 4.2 are consistent with a substantially larger reduction in effective off-chip traffic. CuPy runtime compilation. The fused kernel is written as a CUDA C string and compiled at first use through CuPy’s runtime kernel-compilation interface, which invokes NVRTC (NVIDIA Runtime Compilation) without requiring a standalone CUDA toolchain. The kernel string is approximately 80 lines of CUDA C, embedded directly in the Python source file; no separate CUDA source file or offline compiler invocation is required. This design preserves the fully Python-native workflow of the three-stage baseline used throughout the paper while adding near-native CUDA performance.
3.4
BF16 WMMA Tensor-Core Variant
Modern NVIDIA GPUs expose tensor cores through the WMMA (Warp-Level Matrix Multiply– Accumulate) API, which computes C ← AB + C for small fixed-size tile shapes in hardwareaccelerated mixed-precision arithmetic. On the RTX 4090 (Ada Lovelace architecture) the supported BF16 WMMA shape is 16 × 16 × 16 (m,n,k), accumulating into FP32. BF16 represents values with a 7-bit mantissa and 8-bit exponent (same dynamic range as FP32) and achieves 2× the peak throughput of FP32 CUDA cores on the RTX 4090: 165.2 TFLOP/s (BF16 tensor cores) versus 82.6 TFLOP/s (FP32).
Figure 2: BF16 WMMA tensor-core tiling for the per-element stiffness multiply. Left: the 24 × 24 element stiffness matrix Kunit is zero-padded to 32 × 32 and covered by four 16 × 16 e output tiles, each accumulated from two 16 × 16 × 16 WMMA multiplies along the padded k dimension. Right: 128 elements are batched per thread block (8 warps), with 16 elements processed per warp and 4 tensor-core matrix-multiply calls per warp. The WMMA multiply uses BF16 fragments; FP32 accumulation and the final atomic scatter maintain output fidelity. Kernel design. The per-element GEMM fe = ke Kunit ue involves a 24 × 24 matrix times a e 24 × 1 vector. Directly mapping this to WMMA requires padding to multiples of 16: Kunit is e zero-padded to 32 × 32, and ue is padded to 32 × 16 (one column per element with 15 zero columns). Each CUDA thread block processes 128 elements, with 8 warps per block and 16 10
elements assigned to each warp. The launch configuration is therefore (⌈nelem /128⌉, 1, 1) blocks × (256, 1, 1) threads. Within each warp, 4 tensor-core matrix-multiply calls tile the 32 × 32 output with 16 × 16 fragments, yielding the product Kunit [ue0 | . . . | ue15 ] simultaneously for all e 16 elements in the block. The SIMP scaling ke is applied post-WMMA in FP32, and the scatter uses atomic accumulation with FP32 precision. Memory architecture implications. The BF16 WMMA path stores the gathered displacement values and the padded element matrix in FP32 and casts them to BF16 immediately before loading WMMA fragments. The SIMP penalty vector is maintained in FP32 throughout, and the scatter output is written in FP32. When this kernel is embedded in the mixedprecision solver, the resulting matvec is used inside an FP32 CG loop; in the BF16 experiments reported later, an FP32 outer residual-correction loop wraps that inner BF16 solve. This FP32→BF16→FP32 mixed-precision cascade is the defining characteristic of the BF16 integration path and the source of the convergence difficulties analyzed in the next section.
3.5
BF16 Arithmetic in the Conjugate Gradient Solver: Convergence Analysis
Residual-correction formulation used in the reported experiments. The reported BF16 experiments use a practical two-precision residual-correction loop inspired by iterative refinement [20]: r(k) = f − Ku(k) e(k) ≈ K−1 r(k) u
(k+1)
=u
(k)
+e
(FP32), (BF16 inner CG),
(k)
.
(7) (8) (9)
The classical convergence lens is still the standard εℓ · κ(K) ≤ c for a modest constant c < 1 and working-precision unit roundoff εℓ [20]. For BF16, we use the Carson–Higham unit-roundoff convention εBF16 = 2−8 ≈ 3.9 × 10−3 throughout. Condition number of the SIMP stiffness matrix. The condition number κ(K) grows with both the SIMP penalty exponent and the density contrast between solid and void elements. For a well-converged SIMP design with large penalization and near-void regions, the modulus contrast between solid (ρe = 1) and near-void (ρe ≈ ρmin ) elements grows as 1/ρpmin , driving κ(K) well above the BF16 threshold at moderate penalization. A companion condition-numberestimation study estimates κ(K) via matrix-free power iteration for the largest eigenvalue and inverse iteration for the smallest, evaluating the uniform-density initialization states (ρ = 0.5) across three mesh sizes (64 k, 216 k, and 512 k) and two penalization levels (p ∈ {3.0, 5.0}). No explicit K is assembled; both extremal eigenvalues are estimated entirely through matrixfree Kv applications (Algorithm 1), making the procedure available at all tested sizes. Power iteration targets a relative change tolerance of 10−6 on successive eigenvalue estimates, and inverse iteration uses a 10−4 relative change tolerance on the reciprocal Rayleigh estimate; the supplementary tabulation records the realized outer-iteration counts for each row. The current workflow uses fixed random start vectors (seed 42 for power iteration and seed 123 for inverse iteration); each inverse-iteration inner solve uses SciPy CG with relative and absolute tolerance 10−8 and a 2,000-iteration cap. In the reported rows, the power iteration reaches the 50-step cap before meeting the 10−6 target, so the reported λmax and κ values should be read as conservative estimates. The reported estimates indicate that κ(K) exceeds 1/εBF16 ≈ 256 for the benchmark systems studied here (see Section 4.8). Convergence barrier. The BF16 convergence condition requires: 1 εBF16 · κ(K) < 1 ⇐⇒ κ(K) < ≈ 256. (10) εBF16 Direct power-iteration estimates (Section 4.8) indicate that the tested systems lie firmly and deeply in the non-convergence zone: κ(K) ≈ 6.1×105 at 64 k, ≈ 1.3×106 at 216 k, and ≈ 2.3×106 11
at 512 k elements, giving εBF16 ·κ ≈ 2.4×103 –9.1×103 —more than an order of magnitude above the threshold across all tested sizes. The BF16 residual-correction experiments in Section 4.7 stagnate at large compliance error precisely because the inner BF16 solve is unable to reduce the outer FP32 residual below the BF16 noise floor in this regime. Path forward. The BF16 GEMM-stage throughput advantage (the separate 14.3× PyTorch BF16 GEMM proxy timing at 512 k elements discussed in Section 4.2) motivates using the BF16 kernel as a multigrid smoother rather than as a preconditioner for the full-system CG. In a V-cycle smoother, the spectrum of the residual presented to each level is bounded by the coarse-grid correction, so the effective κ seen by the BF16 smoother can be kept below the 1/εBF16 ≈ 256 threshold. This direction is identified as the primary future extension in Section 5.
3.6
Preconditioned Conjugate Gradient and Warm-Start Strategy
The global linear system Ku = f is solved at each SIMP iteration using Jacobi-preconditioned conjugate gradient (PCG) [55, 56]. The Jacobi preconditioner M = diag(K) is formed exactly in O(nelem ) via a secondary matrix-free pass that accumulates diagonal contributions element-byelement, without any explicit assembly. Despite its simplicity, Jacobi preconditioning requires p O( κ(K)) CG iterations for convergence; in practice this yields a few hundred iterations per SIMP step for the cantilever problem, while the torsion benchmark often approaches the 1,000iteration cap. The Jacobi PCG used in the core FP64/FP32/fused solver path admits a pure Python/CuPy implementation with no standalone CUDA build step. The separate tensor-core GEMM proxy benchmark reported later uses PyTorch only for timing those BF16/FP16 GEMM proxy calls; it is not part of the core solver path. Warm-start. Successive SIMP iterations produce design fields that evolve slowly once past the first 10–20 iterations; the displacement solution u(k) from iteration k is therefore used as the initial guess for iteration k + 1. Warm-start is enabled in all quantitative SIMP benchmarks and materially reduces the later-iteration CG workload relative to repeated cold starts. Convergence criterion. In the reported implementation, CG is run until the unpreconditioned relative residual satisfies ∥rk ∥/∥f ∥ ≤ 10−5 , or a maximum of 1,000 iterations is reached. This is the stopping rule used throughout the reported experiments and across the FP32/BF16 matrix-free solver variants.
3.7
Roofline Analysis
The roofline model [15, 57] bounds kernel throughput by the minimum of compute peak (Π, in FLOP/s) and memory bandwidth times arithmetic intensity (I, in FLOP/B): P ≤ min(Π, I ·b), where b is the memory bandwidth ceiling. Table 1 summarizes the RTX 4090 performance ceilings and the arithmetic intensity of the fused and BF16 kernels. The idealized arithmetic intensity of the fused FP32 kernel is estimated as Iideal,fused =
2 × 242 FLOP ≈ 6 FLOP/B, 2 × 24 × 4 bytes
and a coarse implementation-level accounting gives an effective Iideal ≈ 5.8 FLOP/B—well below the FP32 ridge point of 81.9 FLOP/B. The fused kernel is therefore DRAM-bandwidth-bound across all tested sizes, and the roofline-level speedup over the three-stage baseline is explained by lower DRAM traffic rather than by higher arithmetic intensity. Under the profiling traffic model that also counts edof index reads and per-element density reads, the reported values are Iprofile = 1.33 (FP64 three-stage) and 3.95 (fused FP32) in Table 3. The reported BF16versus-FP64 14.3× GEMM proxy timing at 512 k elements shows that the isolated GEMM sub-operation can be accelerated dramatically by BF16 tensor-core arithmetic in a separate PyTorch proxy benchmark. The paper does not report a separate stage-level roofline derivation 12
Table 1: RTX 4090 roofline ceilings and idealized arithmetic intensity (Iideal ) for the full-kernel variants reported in the paper [16]. Dashes indicate hardware-ceiling rows for which arithmetic intensity is not applicable. All three full-kernel variants remain DRAM-bandwidth-bound on the tested RTX 4090. Section 4.2 separately reports profiling-traffic arithmetic intensity Iprofile , which includes index and density reads.
Variant
Π (TFLOP/s)
b (TB/s)
Ridge (FLOP/B)
Iideal (FLOP/B)
1.29 82.6 165.2
1.008 1.008 1.008
1.29 81.9 163.9
— — —
— — —
1.008 1.008 1.008
— — —
≈ 3.2 ≈ 5.8 ≈ 5.8
FP64 CUDA cores FP32 CUDA cores BF16 tensor cores Baseline 3-stage (FP64) Fused FP32 kernel Fused BF16 full kernel
for that GEMM proxy. The only full-kernel roofline points used in the paper are those listed in Table 1, and under that accounting the full fused-BF16 kernel remains bandwidth-bound because gather and scatter preserve the same irregular global-memory traffic as the fused FP32 path.
3.8
Software Architecture
The complete solver is implemented in Python 3.11 and CuPy 13 and runs on a single consumer GPU (NVIDIA RTX 4090, 24 GB GDDR6X) without any standalone CUDA build step. The software layers are: 1. Problem setup: mesh generation, DOF table (edof) construction, boundary condition assembly, and filter kernel construction—all in NumPy/CuPy. 2. Fused CUDA kernel: the CUDA C string is compiled once at first use through CuPy’s runtime kernel-compilation interface and cached; subsequent calls reuse the compiled binary. 3. PCG solver loop: pure Python with CuPy array operations; the fused kernel is called as a function object for the Kv product. 4. SIMP outer loop: OC update, filter application, Heaviside projection, and convergence check in NumPy/CuPy. 5. Results I/O: density fields and benchmark metrics are serialized to NumPy array and CSV formats for post-processing. All data remain on the GPU throughout the inner PCG loop; only scalar convergence monitors and per-SIMP-iteration compliance values are transferred to the CPU through ordinary scalar conversion. In the reported SIMP scaling runs, the fused path uses about 5.0 GB at nelem ≈ 2 × 106 and 10.56 GB for an 8 M-element FEA-only solve, leaving substantial headroom within the 24 GB VRAM limit.
4
Experimental Results
4.1
Experimental Setup
Hardware. All GPU experiments are conducted on a single NVIDIA RTX 4090 [16] (Ada Lovelace; 24 GB GDDR6X; 1.008 TB/s bandwidth; 165.2 / 82.6 / 1.29 TFLOP/s for BF16/FP32/FP64 respectively). The reported results were generated on Microsoft Windows 10.0.26200.8037 with NVIDIA driver 595.71, in Python 3.11 under CuPy 13.6.0, PyTorch 2.5.1 built against CUDA 12.1, NumPy 2.2.6, SciPy 1.15.3, Matplotlib 3.10.7, scikit-image 0.25.2, Pandas 2.3.3, and PyVista 0.46.3; the CuPy runtime reports CUDA runtime version 12090, and the 13
board uses the default 450 W power limit. This CUDA 12.1/12.9 split reflects PyTorch’s bundled CUDA 12.1 user-space runtime alongside the installed CUDA 12.9 runtime used by CuPy; that mixed configuration is supported by the installed NVIDIA driver on this host. The core solver path uses CuPy; PyTorch is used only in the hot-path microbenchmark script for the separate BF16/FP16 GEMM proxy timings, while the render scripts additionally use PyVista and scikit-image. Benchmark problems. Two main quantitative benchmark families are used throughout (cantilever and torsion), together with two 187,500-element hard-problem stress tests (MBB and bridge): • Cantilever beam: domain 2×1×0.5 with fixed left face and a unit downward point load at the right-face midpoint. All cantilever runs use Vf = 0.30 and an initial filter radius rmin = 1.5. The hot-path microbenchmark covers 64 k, 216 k, 512 k, and 1 M elements. The canonical cold-start FEA scaling figure covers 216 k, 512 k, 1 M, 2.0 M, 4.9 M, and 8 M elements, while the end-to-end SIMP scaling study reports 216 k, 512 k, 1 M, 2.0 M, and 4.9 M elements. • MBB hard-problem stress test (150 × 50 × 25 = 187,500 elements; domain 3 × 1 × 0.5; ux constrained along the left edge, uy constrained at (3, 0, 0.25), and a unit downward point load applied at (0, 1, 0.25); Vf = 0.50; initial rmin = 1.5). • Bridge hard-problem stress test (150×50×25 = 187,500 elements; domain 3×1×0.5; fixed support points on the left lower edge, roller-x support points on the right lower edge, and a distributed downward load on the top edge; Vf = 0.30; initial rmin = 1.5). • Torsion shaft (165 × 55 × 55 = 499,125 elements; fixed left face, equal-and-opposite torsional loads on the top and bottom edges of the right face; Vf = 0.25; initial rmin = 1.5). An additional left-clamped, right-roller beam example is shown later as a qualitative render companion to the bridge load family; it is not the source of the bridge timing table. The main cantilever and torsion optimization studies use 120 outer iterations and the same fixed four-phase continuation schedule: iterations 1–15 use (p, β, m) = (1.5, 1.0, 0.20); iterations 16–40 use (3.5, 4.0, 0.15) while reducing rmin toward 1.35; iterations 41–65 use (4.5, 16.0, 0.08) while reducing rmin toward 1.25; and iterations 66–120 use (4.5, 32.0, 0.05) while reducing rmin toward 1.20, where m is the OC move limit. If compliance rises above 1.12× the current selected compliance, the run restarts from the currently selected design field. Table 6 is intentionally separate from this main protocol: its truncated hard-problem stress test uses 60 outer iterations, cold-start initialization, and no warm-start. For the main SIMP-120 studies, reported compliance values and rendered topologies therefore correspond to the selected design iterate from each run: the lowest-compliance iterate that passes the run’s validity checks, not necessarily the final iterate. Appendix B defines this selected-versus-final convention formally for all reader-facing tables and figures; the hard-problem table uses the same best-valid selection rule, but on the shorter 60-step cold-start schedule stated in its caption. In the present implementation, an iterate becomes eligible for this selection only once p ≥ 3.0 and the grayness metric satisfies g < 0.25. Here the grayness metric is g=
4
nX elem
nelem e=1
ρe (1 − ρe ),
so g = 0 denotes a fully binary design and larger values indicate more intermediate-density material. All SIMP runs use deterministic uniform-density initialization, so no random seed is required for the optimization path itself. Solver variants. Three implementations are benchmarked: fp64 (three-stage baseline: gather + CuPy-dispatched batched DGEMM + histogram-style scatter reduction, FP64 throughout), fp32 (same three-stage pipeline in FP32 state, with the histogram-style scatter reduction accumulated in float64 before cast-back), and fused (the single fused gather–GEMM–scatter kernel in FP32). A fourth variant bf16-wmma is benchmarked for full fused-kernel matvec through14
put, for the separate BF16 GEMM proxy timing reported in Table 2, and for representative mixed-precision linear-solve experiments. Warm-start is enabled for the main cantilever and torsion SIMP-120 benchmarks (Section 3.6); the hard-problem SIMP-60 rows in Table 6 are explicit cold-start runs without warm-start. Timing is measured with CUDA events (deviceside) for kernel profiling and wall-clock time for end-to-end SIMP benchmarks. The end-to-end SIMP wall times in Tables 4 and 5 are single representative runs rather than sample means; Tables 8–10 report the separate repeat, determinism, and high-cap validation studies, and the bridge hard-problem table states its mixed mean/single-run reporting explicitly in the caption.
4.2
Hot-Path Microbenchmark and Single-Solve Scaling
Table 2 reports a synthetic hot-path microbenchmark for the gather/GEMM/scatter operator shape at four cantilever mesh sizes. The source script seeds random synthetic data with the same tensor shapes as the real operator, including random DOF-index patterns rather than the structured cantilever connectivity, using a fixed seed of 42 and an adaptive repeat count that targets about 50 ms of total measurement per size. Figure 4 then reports cold-start wall-clock times for a full FEA solve at uniform density ρ = 0.5, providing a size-scaling check of the cold-start linear-solve cost. Table 2: Synthetic hot-path timing breakdown (µs) on the RTX 4090. tfull64 : three-stage FP64 pipeline (gather + batched DGEMM + scatter). tfused : single fused FP32 kernel. tbf16 : single fused BF16 WMMA kernel. Speedup is relative to tfull64 .
Size (nelem )
tgather
tDGEMM64
tscatter
tfull64
tfused
tbf16
Fused speedup
BF16 speedup
64 k 216 k 512 k 1M
28.6 70.8 183.4 351.3
140.8 452.3 1154.3 2280.7
167.1 315.6 731.2 1423.2
342.6 1052.2 2409.0 4615.8
56.9 168.3 364.1 765.3
58.3 166.3 361.9 739.0
6.0× 6.3× 6.6× 6.0×
5.9× 6.3× 6.7× 6.2×
In the synthetic operator microbenchmark, the fused kernel achieves a 6.0–6.6× speedup over the FP64 three-stage baseline across the full size range. Profiling-guided traffic accounting indicates that the fused kernel reduces effective off-chip traffic substantially, consistent with the theoretical analysis in Section 3.3. The remaining gap between the theoretical 2× elementdata reduction (or roughly 1.8× once index and density reads are included) and the observed 6× speedup is consistent with two additional effects: (i) elimination of kernel-launch overhead across the three GPU-resident stages and (ii) improved warp utilization from the regular 128thread block configuration of the fused FP32 kernel. Because the stage timings in Table 2 are measured independently, the stacked component sum in Figure 3 is a lower bound on the full FP64 pipeline time rather than a complete decomposition of every pipeline overhead. GEMM-stage BF16 proxy timing. While the total per-matvec speedup of the BF16 WMMA kernel is comparable to the fused FP32 kernel (both around 6× vs. the FP64 pipeline), the GEMM proxy benchmark reported by the separate profiling study shows dramatic improvement: at 512 k elements, the BF16 GEMM proxy takes 80.6 µs versus 1154.3 µs for the FP64 batched DGEMM—a 14.3× BF16-versus-FP64 GEMM proxy speedup. This proxy is measured separately with PyTorch BF16 GEMM on the same tensor shapes; it is not an instrumented stage timing extracted from the custom CuPy WMMA kernel itself. The reason the full-pipeline speedup saturates near 6× is that the gather and scatter stages (which still stage through the global-vector and reduction paths) become the dominant cost once the GEMM is accelerated by tensor cores: at 512 k elements, gather+scatter consumes 914.6 µs of the FP64 pipeline but only about 281.3 µs of non-GEMM work in the fused BF16 path
15
Synthetic hot-path breakdown
Time per matvec (µs)
4000 3000
FP64 full pipeline Gather GEMM (FP64) Scatter-add Fused FP32 (single kernel) Fused BF16 (WMMA)
2000 1000 0
6.0× 6.2× 6.0× 5.9×
6.3× 6.3×
64k
216k
6.6× 6.7×
512k
1M
Figure 3: Synthetic hot-path timing breakdown at four element counts. Bars show the FP64 gather, batched DGEMM, and scatter components from the three-stage microbenchmark, together with the corresponding fused FP32 and fused BF16 kernels measured on the same synthetic input shapes. This figure illustrates operator-level cost composition; it is not a full SIMP or full FEA timing plot. The stacked gather/DGEMM/scatter bars represent the sum of separately measured stage durations, while the black markers indicate the full FP64 pipeline time from Table 2; the stage sum is therefore a lower bound on the complete FP64 path. Its gather/scatter locality should be read as a synthetic stress test rather than as a literal replay of the structured cantilever connectivity pattern. (361.9 − 80.6 ≈ 281.3)1 —so even with zero-cost GEMM, the achievable speedup would be capped at 2409/281.3 ≈ 8.6×. This analysis motivates the future direction of a fully fused BF16 kernel that also eliminates the FP64 gather–scatter round-trip. Bandwidth utilization on the actual operator. To complement the synthetic microbenchmark, a separate CUDA-event measurement on the actual cantilever operator (50 timed iterations with 10 warm-up iterations) yields effective memory bandwidths of 61–179 GB/s (6–18% of the RTX 4090’s 1,008 GB/s peak DRAM bandwidth) for the fused FP32 kernel, versus 20– 43 GB/s (2–4% of peak) for the FP64 three-stage baseline (Table 3). Both paths sit far below the roofline ridge point of ≈ 82 FLOP/B (82.6 TFLOP/s ÷ 1,008 GB/s), indicating that the gather–GEMM–scatter operator is memory-bandwidth limited regardless of precision. The isolated per-call speedup on the actual cantilever operator is 8.9–13.8×—higher than the 6.0–6.6× from the synthetic benchmark—because the actual Python/CuPy dispatch overhead is also eliminated: the fused kernel replaces three separate CuPy API calls (one per stage) with a single runtime-compiled kernel invocation. 1
Computed from the canonical 512 k profiler row: full fused BF16 kernel minus the separate BF16 GEMM proxy timing. This subtraction is a useful upper-bound accounting argument, not an instrumented stage timing inside the CuPy WMMA kernel.
16
Mean FEA solve wall time (s)
Cold-start FEA solve scaling FP64 three-stage baseline FP32 three-stage baseline Fused FP32 kernel 101
100
106
Number of elements Figure 4: Cold-start FEA solve wall time versus element count for the three cantilever solver variants (RTX 4090, log–log scale) at uniform density ρ = 0.5 for 216,000 to 8,000,000 elements. Each point is a full Jacobi-PCG solve, not a single matvec timing. The 2 M, 4.9 M, and 8 M points reach the current 1,000-iteration cap and should be read as capped stress-test timings rather than as fully converged solves. The 8 M point is available only for the FP32 and fused paths in the current study. The fused path preserves the same linear-solve algorithm while reducing the cost of each operator application. Table 3: Memory bandwidth utilization of the actual cantilever operator, measured with CUDA events on the RTX 4090 (50 timed iterations, 10 warm-up). Iprofile = profiling-traffic arithmetic intensity (FLOP/byte), including index and density reads. Effective bandwidth is theoretical bytes moved divided by measured wall time under the profiling traffic model.
Size 64 k 216 k 512 k 1M
4.3
FP64 three-stage
Fused FP32
Speedup
t (µs)
BW (GB/s)
Iprofile
t (µs)
BW (GB/s)
Iprofile
2736 4871 11291 20285
20.3 38.5 39.4 42.8
1.33 1.33 1.33 1.33
307 353 833 1700
61.0 178.8 179.4 171.8
3.95 3.95 3.95 3.95
8.9× 13.8× 13.6× 11.9×
End-to-End SIMP Scaling: Cantilever Benchmark
Table 4 reports full SIMP-120 wall times and compliance values for the cantilever benchmark across five mesh sizes. The fused FP32 kernel achieves consistent speedups of 4.6–7.3× over the FP64 three-stage baseline end-to-end. Against the same-precision FP32 three-stage baseline, the same rows correspond to 2.3–4.6× speedups. These rows are representative single-run measurements from the dedicated scaling workflow; Section 4.9 reports a separate five-repeat study at 216 k and 512 k to quantify run-to-run variability.
17
Table 4: End-to-end SIMP-120 results on the cantilever benchmark (RTX 4090). FP64 and FP32 denote the three-stage baselines. Reported compliance values are the selected compliances from the 120-step continuation schedule. Speedup is FP64/Fused. These rows are representative single-run timings from the scaling workflow: each deposited (size, path) pair appears once in that workflow, while repeat variability at 216 k and 512 k is reported separately in Table 8. The 2 M fused compliance deviates by 2.46% from FP64; see Section 5. Size 216 k 512 k 1M 2M 4.9 M
nelem 216,000 512,000 1,000,000 2,000,376 4,913,000
FP64
FP32
Fused FP32
t (s)
c
t (s)
c
t (s)
c
Speedup
80.4 137.9 226.2 2610.2 6140.9
2.172 1.825 1.606 1.421 1.046
80.5 57.9 120.5 1190.4 3000.6
2.174 1.828 1.608 1.407 1.040
17.5 24.8 45.4 358.0 996.9
2.174 1.828 1.608 1.386 1.045
4.6× 5.6× 5.0× 7.3× 6.2×
Several observations merit discussion. First, the FP32 path shows negligible difference relative to FP64 at 216 k elements (80.5 s vs. 80.4 s) but achieves 2.4× at 512 k and 1.9× at 1 M. This is consistent with a bandwidth-bound workload in which FP32 halves the per-entry data footprint relative to FP64. At 216 k, the non-FEA SIMP overhead is still large enough to dilute that matvec-level advantage in the full end-to-end wall time. Second, the fused kernel shows a higher relative speedup at 2 M (7.3×) than at 1 M (5.0×), for which one plausible explanation is cache behavior: the Kunit matrix (2.25 kB shared memory) is reused across all 2 M thread blocks, e but at smaller sizes the per-block overhead of loading the kernel launch context dominates. We do not have cache-counter measurements for this interpretation, so this should be read as a hypothesis rather than as a verified mechanism. Third, the selected compliance values remain within about 0.2% of the FP64 reference through 1 M elements. At 2 M, the fused path deviates by 2.46% and the FP32 path by 0.97%, so the exact values are reported explicitly rather than summarized as decimal-place parity. External CPU baseline comparison. To provide limited timing context against a publicly available Python-based 3D SIMP code, Figure 7 compares the fused RTX 4090 path with local PyTopo3D reruns under its CPU/PyPardiso configuration [58] on the same cantilever geometry and mesh sizes. This is a timing-only comparison: the two implementations use different SIMP formulations and mesh-scaling conventions, so compliance values are not directly comparable. The deposited PyTopo3D rerun log records 20 CPU cores and 69.5 GB available RAM for that host, and the accompanying environment manifest records PyTopo3D 0.1.0, PyPardiso 0.4.7, and default host scheduling with no explicit thread pinning; the figure should therefore be read as contextual timing evidence rather than as a version-controlled head-to-head benchmark. At 64 k and 216 k elements, the reported PyTopo3D rows use completed 120-iteration runs; the fused GPU path is faster by approximately 58× and 531×, respectively, but these ratios should be read only as local timing context rather than as a formulation-matched benchmark.
18
SIMP-120 end-to-end wall time (s)
SIMP-120 end-to-end scaling FP64 three-stage baseline FP32 three-stage baseline Fused FP32 kernel 103
102
106
Number of elements Figure 5: End-to-end SIMP-120 wall time versus number of elements for the cantilever benchmark (log–log). The plotted points are the representative single-run timings from the dedicated scaling workflow used for Table 4; Section 4.9 reports the separate five-repeat variability study at 216 k and 512 k. All three paths show near-linear scaling (O(nelem )). The fused FP32 path achieves 4.6–7.3× speedup over FP64; the larger gap at 2 M and 4.9 M is an observed trend for which the paper does not report cache-counter evidence.
Baseline-to-fused end-to-end speedup 7
SIMP-120 speedup
6 5
7.3×
FP64 / Fused FP32 / Fused
6.2×
5.6× 5.0×
4.6× 4.6×
4
3.3×
3
2.7×
2.3×
2
3.0×
1 0
216k
512k
1.0M
2.0M
4.9M
Figure 6: Baseline-to-fused SIMP-120 speedup at five element counts. The blue bars report FP64/Fused and the orange bars report FP32/Fused. The bars are computed from the representative single-run scaling rows in Table 4; Section 4.9 gives the separate repeat-study variability at 216 k and 512 k. FP64/Fused grows from 4.6× at 216 k to 7.3× at 2 M and 6.2× at 4.9 M. 19
Contextual wall-time comparison on the cantilever benchmark 531× PyTopo3D rerun (CPU/PyPardiso, 20-core) Present work, fused FP32 (RTX 4090)
SIMP-120 wall time (s)
104
58×
103 102 101 100
64k
216k
Figure 7: Wall-time comparison against local CPU/PyPardiso PyTopo3D reruns [58] on the cantilever geometry and mesh sizes. This is a timing-only comparison: the two implementations use different SIMP formulations and mesh-scaling conventions, so compliance values are not directly comparable. The deposited PyTopo3D rerun log records 20 CPU cores and 69.5 GB available RAM, and the accompanying environment manifest records PyTopo3D 0.1.0, PyPardiso 0.4.7, and default host scheduling with no explicit thread pinning, so the figure is contextual only. The 64 k GPU bar comes from a separate deposited 64 k fused rerun; the 216 k GPU bar comes from the main cantilever scaling ladder. The 64 k and 216 k bars use completed 120-iteration runs; the 512 k PyTopo3D run was terminated before completion and is excluded from this bar chart.
(a) Side profile (XY plane)
(b) Isometric view
Figure 8: Selected cantilever topology, 120 × 60 × 30 = 216,000 elements (selected design iterate from a separate same-configuration fused FP32 SIMP-120 rerun, Vf = 0.30; smoothed marchingcubes isosurface at ρ = 0.5 with 80-step Laplacian smoothing, relaxation factor 0.08). Left: side-profile view (XY plane) showing the classic arrowhead-arch truss with two oval cutouts. Right: isometric view revealing the three-dimensional ribbed shell structure. Grayness = 0.000; compliance = 2.174.
20
(a) Side profile (XY plane)
(b) Isometric view
Figure 9: Selected cantilever topology at 200 × 100 × 50 = 1,000,000 elements (selected design iterate from a separate same-configuration fused FP32 SIMP-120 rerun, Vf = 0.30; smoothed marching-cubes isosurface at ρ = 0.5 with 80-step Laplacian smoothing, relaxation factor 0.08). Left: side-profile view. Right: isometric view. The higher resolution resolves a cleaner arrowhead arch at the same qualitative design family; the corresponding table-supporting fused run in Table 4 takes 45.4 s. Grayness = 0.000; compliance = 1.608.
4.4
Torsion Benchmark
The torsion benchmark serves as a difficult stress test with much higher CG iteration counts than the cantilever. Table 5 summarizes results for the 165 × 55 × 55 = 499,125 element mesh. Table 5: SIMP-120 torsion benchmark (nelem = 499,125, RTX 4090). FP64 and FP32 denote the three-stage baselines. Reported compliance is the selected compliance from the 120-step continuation schedule. In the reported histories, 108–110 of the 120 iterations hit the current 1,000-iteration CG cap. Path
Wall time (s)
VRAM (GB)
Compliance
Speedup vs. FP64
FP64 FP32 Fused
643.5 415.0 146.5
2.46 2.32 2.38
2.0357 2.0360 2.0360
1.0× 1.6× 4.4×
The torsion problem requires substantially more CG iterations per SIMP step than the cantilever at a comparable size, which is why the absolute wall times are larger than for the cantilever at a comparable 512 k scale. In fact, 108 of 120 FP64 iterations and 110 of 120 FP32/fused iterations reach the current 1,000-iteration cap, with the first cap appearing at iteration 2 in all three histories. The torsion history should therefore be read as a difficult capped stress test rather than as a clean low-iteration regime. Even with that harder linearsolve profile, the fused kernel still delivers a 4.4× end-to-end speedup over the FP64 torsion baseline.
21
Torsion SIMP-120 convergence (499,125 elements) FP32 three-stage baseline FP64 three-stage baseline Fused FP32 kernel
Compliance
102
101
1000
CG iters
900 800 FP32 three-stage baseline FP64 three-stage baseline Fused FP32 kernel
700 600 0
20
40
60
SIMP iteration
80
100
120
Figure 10: SIMP-120 convergence history for the 499,125-element torsion benchmark. Top: compliance versus SIMP iteration (log scale) — all three paths follow nearly identical trajectories throughout the continuation. Bottom: CG iteration count per SIMP step — the high CG counts reflect the elevated difficulty of the torsion linear systems, and 108–110 of the reported iterations hit the current 1,000-iteration cap. The compliance axis uses a log scale to show the full trajectory over the large early-to-late drop.
(a) End view (YZ projection)
(b) Isometric view
Figure 11: Selected topology for the 499,125-element torsion shaft benchmark (selected design from a separate same-configuration fused FP32 SIMP-120 rerun, Vf = 0.25; smoothed marchingcubes isosurface at ρ = 0.5 with 80-step Laplacian smoothing, relaxation factor 0.08). Left: end-view projection (YZ view) of the selected density field. Right: isometric surface rendering of the same selected field. This smoothed isosurface is a qualitative shape illustration rather than direct evidence of the internal cavity geometry. Grayness = 4.1×10−5 ; selected compliance = 2.0360.
22
4.5
MBB and Bridge Hard-Problem Stress Test
Table 6: Hard-problem stress test on the RTX 4090 (nelem = 187,500). FP64 and FP32 denote the three-stage baselines. The FEA-only rows report mean wall time over five timed cold-start solves after warmup. The SIMP-60 cold-start rows report single representative 60-iteration runs without warm-start. For the FEA-only rows, compliance parity is computed from the final uniform-density cold-start solve compliance relative to FP64. For the SIMP-60 rows, it is computed from the selected best-valid compliance relative to FP64. Problem
Mode
MBB Bridge MBB Bridge
FEA-only FEA-only SIMP-60 cold-start SIMP-60 cold-start
Unit
FP64
FP32
Speedup
Compliance parity
ms ms s s
1138.0 1130.8 74.13 75.24
754.7 890.0 49.75 49.78
1.51× 1.27× 1.49× 1.51×
2.61% 0.49% 1.74% 0.16%
The hard-problem study includes both the bending-dominated MBB case and the distributedload bridge case at 150 × 50 × 25 = 187,500 elements. These two presets share the same mesh size and Jacobi-PCG implementation, but they are not like-for-like boundary-condition matches: MBB uses Vf = 0.50 with a concentrated point load and mixed pin constraints, whereas the bridge uses Vf = 0.30 with fixed/roller supports and a distributed top load. The rows should therefore be read as two separate hard-case probes rather than as a controlled MBB-versus-bridge comparison. In the FEA-only rows, all four cold-start solves hit the current 1,000-iteration cap; the SIMP-60 rows capture the short continuation runs. The bridge-side rows indicate that the FP32-versus-FP64 three-stage trend is not specific to the cantilever boundary condition in the single distributed-load case tested here: the FP32 three-stage path remains faster than FP64 on both the cold-start FEA solve and the truncated SIMP continuation, while the compliance drift remains small. This evidence is limited to the three-stage paths; the paper does not report a direct fused-kernel bridge benchmark.
4.6
Additional Qualitative Structures
Figure 12 shows a left-clamped, right-roller beam with central load as a qualitative companion to the bridge stress test above. This run is included to illustrate topology diversity; it is not the source of Table 6 and should not be read as a duplicate quantitative measurement. An additional centered-patch cantilever topology at Vf = 0.10 from a separate fused-FP32 SIMP rerun is included in Appendix A as a supplemental qualitative example outside the main timing tables. Grayness metric. For the qualitative figures, grayness is reported using the definition given in Section 4.1; g = 0 denotes a fully binary design and larger values indicate more intermediatedensity material.
23
(a) Front elevation (XY plane)
(b) Isometric view
Figure 12: Separate qualitative exemplar from the bridge-family load class: selected topology for a left-clamped, right-roller beam with central load, 120 × 60 × 30 = 216,000 elements (selected iterate from a separate fused FP32 SIMP-120 run with Vf = 0.30 and initial filter radius rmin = 3.0 for this qualitative exemplar only; smoothed marching-cubes isosurface at ρ = 0.5 with 80-step Laplacian smoothing, relaxation factor 0.08). The solver recovers an arch-andstrut system transferring the center-top load to the fixed left wall and the roller right support — a structurally distinct topology from the quantitative cantilever, torsion, and 187,500-element bridge timing benchmarks, and a qualitative companion to the bridge hard-problem stress test reported above. This panel is qualitative only and is not the quantitative source for Table 6. Grayness = 0.000; compliance = 0.903.
4.7
BF16 Convergence Study
We investigate the convergence behavior of the BF16 WMMA kernel when integrated into the CG solver via the reported FP32/BF16 residual-correction experiments. Table 7 reports compliance and relative error for four solver configurations for a single cold-start linear solve on the cantilever benchmark at 64 k and 216 k elements, with uniform density ρ = 0.5 and p = 3. These are not full SIMP runs; they are representative linear-solve experiments designed to expose solver convergence behavior. Table 7: BF16 linear-solve convergence study for a single cold-start linear solve on the cantilever benchmark (ρ = 0.5, p = 3). FP32 is the reference. Plain BF16 CG uses no refinement. BF16-IR (inner tol 10−3 ) and BF16-IR (inner tol 10−5 ) denote BF16 inner solves with iterative-refinement outer loops (maximum 8 outer corrections) at two inner tolerances. Compliance relative error δc = |cref − c|/cref . Size
Solver
64 k
FP32 reference Plain BF16 CG BF16-IR (inner tol 10−3 ) BF16-IR (inner tol 10−5 )
Wall (ms) CG iters Compliance
FP32 reference Plain BF16 CG 216 k BF16-IR (10−3 ) BF16-IR (10−5 )
δc
456 186 1201 2144
345 230 1085 1857
13.431 — 6.178 0.540 7.014 0.478 7.014 0.478
591 281 1581 2646
512 234 1257 2044
11.022 — 6.102 0.446 6.424 0.417 6.424 0.417
The results are consistent with the convergence analysis of Section 3.5 in three ways. First, plain BF16 CG produces compliance values with 44–54% error, indicating that the CG solver does not converge to a valid FE solution. The compliance field is therefore numerically invalid, and BF16 CG is unusable as a direct solver for the TO linear systems encountered here. Second, BF16 iterative refinement partially mitigates but does not resolve the problem: compliance error drops from 54% to 48% at 64 k and from 45% to 42% at 216 k, but stagnates 24
there regardless of how tight the inner tolerance is set. Tightening the inner tolerance from 10−3 to 10−5 merely increases CG iteration count (1,085 → 1,857 at 64 k) and wall time (about 1.7–1.8×) without improving the outer residual—a clear signature of the stagnation predicted by the εBF16 · κ(K) ≥ 1 barrier. Third, the stagnation level is problem-size-invariant: at both 64 k and 216 k, BF16-IR converges to the same compliance plateau (δc ≈ 0.42–0.48), indicating that the dominant conditioning effect is more geometry- and penalization-driven than mesh-size-driven. Because the test state is the uniform-density configuration (ρ = 0.5, p = 3), it represents a lower-contrast probe rather than the hardest late-continuation SIMP system; the companion kappa estimation study in Section 4.8 shows the threshold crossing already at all tested uniform states and both reported penalization levels, directly supporting this interpretation. In the single-solve cold-start regime tested here, these findings indicate that the BF16 WMMA kernel is not a viable drop-in replacement for FP32/FP64 in the CG solve for the tested 3D SIMP systems, and provide a focused empirical characterization of this failure mode.
4.8
Condition Number Estimation
To ground the BF16 convergence analysis in a direct measurement rather than an inference from stagnation behavior, the paper includes a dedicated condition-number estimation workflow, which estimates κ(K) via matrix-free power iteration [59] (50 fixed steps for the largest eigenvalue in the current study) and inverse iteration (5–6 steps for the smallest eigenvalue). Both extremal eigenvalues are estimated entirely through matrix-free Kv applications (Algorithm 1), so no explicit K is assembled and the procedure is feasible at all tested mesh sizes. The current workflow uses fixed start vectors (seed 42 for power iteration and seed 123 for inverse iteration), and each inverse-iteration inner solve uses SciPy CG with relative and absolute tolerance 10−8 and a 2,000-iteration cap. The reported tabulation covers the uniform-density initialization (ρ = 0.5) at two penalization levels (p ∈ {3.0, 5.0}) for 64 k, 216 k, and 512 k elements. The reported estimates already indicate that κ(K) far exceeds the BF16 stability threshold at the lowest-contrast test state (uniform density, ρ = 0.5). At 64 k elements, power iteration yields κ ≈ 6.1×105 , giving εBF16 ·κ ≈ 2.4×103 —nearly 9× above the 1/εBF16 = 256 convergence threshold. At 216 k elements the conditioning is larger: κ ≈ 1.3 × 106 and εBF16 · κ ≈ 5.2 × 103 — more than 20× above the threshold. At 512 k elements the trend continues: κ ≈ 2.3 × 106 and εBF16 · κ ≈ 9.1 × 103 —more than 35× above the threshold. Because all six reported rows hit the 50-step power-iteration cap, these κ values should be read as conservative estimates 2/3 for the corresponding uniform states. The three estimates scale approximately as O(nelem ) (κ ratios of 2.18× for the 216 k/64 k pair and 1.75× for the 512 k/216 k pair, consistent with the theoretical h−2 growth for elliptic PDEs on uniform meshes), indicating that the conditioning barrier intensifies monotonically with mesh refinement. All three estimates are insensitive to the penalization exponent at uniform density (identical κ for p = 3 and p = 5): at ρe = 0.5 all element moduli are equal regardless of p, so the condition number is driven solely by the mesh geometry. At late-SIMP density contrasts, near-void elements (ρe ≈ ρmin = 10−9 ) introduce columns with modulus ratio ρpmin ≈ 10−27 , driving κ orders of magnitude higher still. These uniform-density measurements already exceed the BF16 stability threshold by wide margins. Late-SIMP states are expected to be more ill-conditioned still, but that stronger statement is an inference rather than a direct measurement here. A separate uniform-density BF16 extension at p = 4.5 shows the same qualitative stall as at p = 3.0 (δc ≈ 0.54, BF16; δc ≈ 0.48, BF16-IR), with no evidence of improvement at the higher penalization.
25
4.9
Statistical Reproducibility and Run-to-Run Determinism
The end-to-end SIMP tables in the main text report single representative runs. To quantify run-to-run variability, a dedicated repeat-run workflow runs SIMP-120 N = 5 times at each (size, path) pair for the cantilever benchmark and reports mean wall time, standard deviation, and coefficient of variation (CV%). Table 8 reports those deposited repeat-run summaries. Table 8: Five-repeat SIMP-120 variability study for the cantilever benchmark. Each row summarizes the deposited N = 5 reruns at a fixed (size, path) pair. Size
Path
N
Wall mean ± std (s)
CV (%)
c̄
cstd
216 k 216 k 512 k 512 k
FP64 Fused FP64 Fused
5 5 5 5
66.39 ± 6.30 39.58 ± 11.82 111.45 ± 9.34 26.48 ± 0.48
9.49 29.85 8.38 1.82
2.172194 2.173862 1.825360 1.827701
3.42 × 10−6 1.26 × 10−5 9.10 × 10−6 2.76 × 10−5
Table 8 shows near-zero compliance variation (cstd ≤ 3 × 10−5 , less than 0.002% of the mean compliance) across all four tested (size, path) combinations, indicating negligible run-to-run compliance variation in the tested cases. Wall-time CV is larger: the FP64 baseline shows 9.5% at 216 k (66.4 ± 6.3 s) and 8.4% at 512 k (111.5 ± 9.3 s), reflecting OS-level scheduling and thermal jitter between cold-start reinitializations. The fused path shows only 1.8% CV at 512 k (26.5 ± 0.5 s); the fused kernel was already compiled by the preceding 216 k run within the same process, so all five 512 k reps represent pure post-JIT execution. At 216 k the fused CV rises to 29.9% (39.6±11.8 s): the first cold-start rep includes CuPy NVRTC kernel-compilation overhead, but the observed five-run spread also contains later slow reps, so this row should be read as a mixed cold-start/thermal-variability measurement rather than a pure post-JIT timing study. The representative single-run values in Table 4 come from a separate scaling workflow rather than from this repeat study. In that workflow, each size is executed sequentially as FP64, FP32, then fused within one process, so Table 4 should be read as a representative scaling trajectory, whereas the present table quantifies run-to-run variability. A companion determinism workflow runs N = 10 identical cold-start fused FP32 solves on the same problem and measures the relative compliance spread ∆c/cref = (cmax − cmin )/cFP64 . This study quantifies non-determinism from the FP32 atomic scatter in Algorithm 1, whose floating-point reordering across GPU threads is theoretically non-associative [60]. Table 9 reports the deposited determinism summaries. Table 9: Determinism study for repeated cold-start fused FP32 solves on the cantilever benchmark from the dedicated determinism workflow (separate from the representative residualhistory solve shown in Figure 13). Here ∆c/cFP64 = (cmax − cmin )/cFP64 and δc = |cfused − cFP64 |/cFP64 . Size 64 k 216 k
N 10 10
cmin 13.431176 11.022279
cmax
∆c/cFP64
max δc
CG iters
13.431200 11.022296
1.78 × 10−6
5.94 × 10−4
345 512–513
1.56 × 10−6
1.09 × 10−3
Table 9 shows run-to-run compliance spread of 1.78 × 10−6 relative at 64 k elements and 1.56 × 10−6 at 216 k elements—both well below 10−4 —indicating that the fused FP32 path is deterministic to well within engineering accuracy in practice despite the theoretically nonassociative atomics. The FP32-versus-FP64 absolute compliance offset (δc ≈ 6 × 10−4 at 64 k and 1.1 × 10−3 at 216 k) arises from the precision difference itself and is consistent with the values reported in Table 4. 26
Table 10: Selected-iterate high-cap validation for the cantilever benchmark. Canonical rows use the current 1,000-iteration CG cap; the validation reruns raise that cap to 5,000. Relative difference is |ccanonical − chigh |/ccanonical . Size
Path
Canonical cap
High cap
ccanonical
chigh
Relative difference
216 k 216 k 512 k 512 k
FP64 Fused FP64 Fused
1000 1000 1000 1000
5000 5000 5000 5000
2.172196 2.173877 1.825346 1.827643
2.172196 2.173838 1.824751 1.827066
3.61 × 10−9 1.78 × 10−5 3.26 × 10−4 3.16 × 10−4
Selected-iterate validity. A further validation study compares each path’s selected compliance from the standard 1,000-iteration CG cap against a high-cap 5,000-iteration rerun. This validation currently covers only the 216 k and 512 k FP64 and fused cantilever rows. Table 10 reports the deposited high-cap reruns. Reported results show a maximum deviation of |∆c|/c ≤ 3.26 × 10−4 across all four (size, path) combinations tested (216 k and 512 k; FP64 and Fused): the 216 k FP64 canonical and high-cap compliances are identical to 3.6 × 10−9 relative error; the worst case (512 k FP64) is 3.26 × 10−4 , while the 512 k fused row is 3.16 × 10−4 . All four cases fall below the current 5 × 10−4 acceptance threshold, indicating that the selected compliances reported in the tested 216 k and 512 k FP64/fused rows of Table 4 are not artifacts of the 1,000-iteration cap but represent well-converged optima for the tested SIMP continuation schedule. The larger 2 M and 4.9 M scaling rows are not part of this high-cap validation set.
4.10
VRAM Scaling and Memory Efficiency
Table 11 summarizes end-of-run VRAM allocation snapshots at select sizes. Table 11: End-of-run VRAM allocator snapshots on the RTX 4090 (24 GB total). Values are the deposited used-memory snapshots recorded immediately after the canonical SIMP scaling runs; they are not instrumented peak-memory measurements. The fused path is modestly more memory-efficient than the three-stage FP64 baseline, but the observed savings depend on allocator behavior and on other live buffers in the SIMP loop. nelem
FP64 (GB)
FP32 (GB)
Fused (GB)
Fused saving
216 k 512 k 1M 2M 4.9 M
1.93 2.42 3.43 5.09 11.06
1.86 2.35 3.14 5.13 10.39
1.89 2.41 3.25 5.00 10.07
0.04 GB 0.01 GB 0.18 GB 0.09 GB 0.99 GB
These snapshots grow approximately linearly with nelem for all paths, as expected from the O(nelem ) matrix-free memory model. The fused path achieves a modest end-of-run saving versus the three-stage FP64 baseline by eliminating the two 24-component per-element work arrays; in the deposited 4.9 M run the observed difference is about 1 GB. The total footprint is not dominated by any single small array: persistent buffers include the global DOF state, the elementto-DOF table, the density field, the sparse density-filter operators, and solver/preconditioner workspace, with the exact split depending on path and allocator reuse. The RTX 4090’s 24 GB budget is consistent with the deposited 8 M-element FEA-only allocator snapshot (10.56 GB used at end-of-run), but that point should be read as a capped memory-allocation stress test rather than as a peak-memory guarantee or a fully converged solve.
27
4.11
Correctness Verification
Numerical correctness of the fused FP32 kernel is assessed by three cross-checks. 1. Compliance parity. Across the reported tables, the fused FP32 selected-summary compliance remains close to the FP64 reference, with sub-0.2% deviations through 1 M elements, a 2.46% deviation at 2 M, and a 0.12% deviation at 4.9 M. The high-cap selectediterate validation currently covers only the 216 k and 512 k cantilever rows (Section 4.9). 2. Selected-field grayness. The reported render metadata for the selected fused density fields shown in Figures 8, 9, and 11 indicate essentially discrete designs (g ≈ 0), indicating that the fused solver reaches the same low-grayness regime as the reference paths. 3. Residual monitoring. Figure 13 shows a representative 216 k cold-start linear solve at uniform density ρ = 0.5. The FP64, FP32, and fused traces are nearly indistinguishable and terminate in 510–512 CG iterations, indicating that the fused operator introduces no visible additional error in this representative case. The supplementary trace file stores the exact per-iteration residual histories for this representative solve.
CG residual decay
cantilever 216k (120x60x30)
100
rk / f
10 1 10 2 10 3 10 4
FP64 three-stage baseline FP32 three-stage baseline Fused FP32 kernel
10 5 0
100
200
300
CG iteration
400
500
Figure 13: Representative CG residual decay for a cold-start linear solve at uniform density ρ = 0.5 on the 120 × 60 × 30 = 216,000-element cantilever benchmark. The plotted quantity is the unpreconditioned relative residual ∥rk ∥/∥f ∥. All three solver paths (FP64, FP32, fused FP32) follow nearly identical trajectories.
5
Discussion
5.1
Interpretation of Speedup Trends
The fused gather–GEMM–scatter kernel achieves 6.0–6.6× per-matvec speedup (synthetic hotpath microbenchmark) and end-to-end SIMP-120 wall-time speedup of 4.6–7.3× on cantilever plus 4.4× on the 499,125-element torsion benchmark over the FP64 three-stage baseline. Against the same-precision FP32 three-stage path, the fused solver delivers 2.3–4.6× on cantilever and 2.8× on torsion. CUDA-event measurements on the actual cantilever operator (Table 3) yield isolated per-call speedups of 8.9–13.8×—higher than the synthetic benchmark. That gap is consistent with the extra Python/CuPy dispatch overhead avoided when three API calls are collapsed into one launch; the synthetic benchmark reports pure GPU kernel time and thus 28
captures only the hardware-level benefit. Both measurements indicate that the operator is memory-bandwidth bound: the fused kernel achieves 61–179 GB/s effective bandwidth (6–18% of the 1,008 GB/s peak), versus 20–43 GB/s (2–4% of peak) for the FP64 baseline—both well below the roofline ridge point of ≈82 FLOP/B. The per-matvec speedup is consistent with two mechanisms that compound. The primary mechanism is DRAM traffic elimination: the fused kernel avoids writing the intermediate per-element work arrays (uelem and felem ) to DRAM between stages, reducing the effective DRAM load according to the profiling-guided traffic accounting in Section 4.2. The secondary mechanism is kernel launch overhead elimination: the baseline requires three separate GPU-resident stages (gather, batched GEMM, reduction scatter), each incurring dispatch and synchronization overhead that accumulates over hundreds of CG iterations per SIMP step. At 64 k elements, where the kernel execution time is short (∼28–168 µs per stage), this overhead is proportionally large, and its elimination accounts for a non-trivial fraction of the total speedup. The higher relative speedup at 2 M versus 1 M elements (7.3× vs. 5.0×) may appear counterintuitive since both sizes share the same kernel architecture. One plausible explanation is cache effects: the Kunit matrix (2.25 kB in the FP32 fused path; 2.0 kB for the BF16 padded e shared-memory tile in the WMMA path) is broadcast from L2 cache to all 2 M thread blocks at runtime; at very large sizes the sheer number of active thread blocks keeps the L2 warm, while at moderate sizes (∼1 M) there is a transition region where the cache-miss penalty from loading Kunit is visible but not yet amortized over enough active blocks. The paper does not include e cache-counter measurements, so this interpretation should be read as a plausible explanation rather than as a directly profiled mechanism, and we therefore treat it as an open performance question rather than as a validated cache account. The FP32 three-stage path achieves about 1.9–2.4× speedup over FP64 for large sizes (512 k– 4.9 M). This is consistent with a bandwidth-bound regime in which FP32 halves the per-entry data footprint relative to FP64, thereby reducing effective DRAM traffic, rather than with the RTX 4090’s much larger FP32:FP64 peak-FLOP ratio. At 216 k, the FP32 path shows essentially no difference relative to FP64 (80.5 s vs. 80.4 s) because the non-FEA SIMP overhead is still large enough to mask that matvec-level FP32 advantage in the full end-to-end timing. The bridge hard-problem stress test shows the same FP32-versus-FP64 three-stage trend under one distributed-load case: 1.27× speedup for the cold-start FEA solve and 1.51× for the 60-iteration SIMP run, with compliance parity staying within 0.5% and 0.2%, respectively. That result does not directly benchmark the fused kernel on the bridge case, but it is consistent with the same FP32-versus-FP64 three-stage trend on the tested bridge case; direct fused-kernel bridge benchmarking remains future work.
5.2
BF16 Tensor Cores: Promise and Precision Barrier
The reported profiling study includes a separate 14.3× BF16-versus-FP64 GEMM proxy timing at 512 k elements (80.6 µs vs. 1154.3 µs for FP64 batched DGEMM), yet the full-pipeline speedup saturates near 6.7×—nearly the same as the fused FP32 kernel. This result reveals an important architectural reality of the gather–GEMM–scatter pattern: the gather and scatter stages, which involve globally irregular memory accesses (driven by the DOF index table) and atomic reductions, are not accelerated by tensor cores. These stages impose an irreducible memory-bandwidth cost that sets the throughput ceiling regardless of how fast the GEMM stage runs. In the current fused kernel architecture, gather + scatter account for ∼196–1775 µs of the FP64 full-pipeline time across the reported 64 k–1 M profiling rows, while the GEMM stage accounts for 141–2281 µs. BF16 WMMA accelerates only the GEMM portion, so the total pipeline speedup is bounded by the non-GEMM share. Using the reported 512 k accounting argument from Section 4.2, even a theoretically zero-cost GEMM would still be capped near 8.6× for that fused-BF16 row.
29
The BF16 convergence failure is a more fundamental concern. Direct power-iteration estimates reported in Section 4.8 measure κ(K) ≈ 6.1 × 105 at 64 k, ≈ 1.3 × 106 at 216 k, and ≈ 2.3 × 106 at 512 k, for the uniform-density test state (ρ = 0.5), giving εBF16 · κ ≈ 2.4 × 103 , 5.2 × 103 , and 9.1 × 103 respectively—well above the 1/εBF16 = 256 sufficient threshold implied by the standard IR bound εBF16 κ(K) < 1 from Carson and Higham [20]. This places the tested systems firmly and deeply in the εBF16 ·κ(K) ≫ 1 regime where the IR guarantee does not apply. Iterative refinement (IR) does not resolve this because IR’s outer residual correction step itself requires a matrix–vector product with K, and when εBF16 · κ(K) ≫ 1, each inner BF16 solve introduces an error that cannot be fully corrected by the outer FP32 residual-correction step. In that regime, the BF16 solve is simply too inaccurate to reduce the outer residual below the BF16 noise floor for the tested systems. This result qualifies the optimistic projection of Henry et al. [22] for the specific tested TO systems. Henry et al. argue that BF16-IR converges “over a large range of condition numbers” when FP32 is used for the outer refinement; our results show that this claim does not extend to the reported 3D structural-elasticity experiments once the effective condition number rises above the BF16 IR threshold. The qualification is not a criticism of BF16-IR in general—for well-conditioned or properly preconditioned systems (e.g., the lattice QCD operator studied by Clark et al. [42], where mixed-precision Krylov solvers with reliable updates succeed in practice)—but it identifies the tested SIMP TO stiffness matrices as a regime where direct BF16-IR fails.
5.3
Contextual Comparison to Prior Implementations
The closest published competitor is Träff et al. [12], whose Futhark/OpenMP-C matrix-free solver achieves 65.5 M elements on an NVIDIA A100 (80 GB HBM2e) with SSOR V-cycle multigrid preconditioning. Our solver achieves 2 M elements on an RTX 4090 (24 GB GDDR6X) with Jacobi preconditioning with a few hundred CG iterations per cantilever SIMP step. The two implementations are not directly comparable for three reasons. First, the hardware differs: the A100 has 80 GB HBM2e at up to 1.94 TB/s versus the 4090’s 24 GB GDDR6X at about 1 TB/s [16, 61]. The VRAM gap contributes materially to the size difference, but it does not by itself explain the full problem-size gap; the 4090’s 24 GB limits the reported single-GPU study to ∼8 M elements in the current FEA-only stress test. Second, the preconditioner gap is responsible for the per-step efficiency difference. Träff et al. report much lower per-step times on comparable A100-scale problems; our Jacobi-PCG requires hundreds of iterations (per-step time 0.15–8 s depending on size), and the torsion stress test reaches the current 1,000-iteration limit in 108–110 of its 120 SIMP steps. This paper should therefore be read as an implementation-focused single-GPU study, not as a preconditionermatched head-to-head benchmark. Third, accessibility differs: Träff et al.’s solver requires the Futhark compiler or an additional build toolchain; our solver is Python-native once CuPy is installed. This accessibility advantage matters for the broad engineering community of TO practitioners who work in Python ecosystems. Wang et al. [14] achieve 128 M elements on a 64 GB CPU workstation using geometric multigrid with non-dyadic Galerkin coarsening. That result underscores the same point as Träff et al.’s work: the main remaining gap in the present implementation is preconditioning rather than matrix-free operator throughput.
5.4
Limitations
The present work has four principal limitations. √ Jacobi preconditioning. The Jacobi preconditioner limits CG convergence to O( κ) iterations, which dominates wall time for problems with high condition numbers. For the 3D MBB 30
beam problem with standard pin-roller boundary conditions, the near-rigid-body mode leads to very poor conditioning, and the accompanying hard-problem MBB rows remain cap-limited under the implemented 1,000-iteration Jacobi-PCG setting. This means the current solver does not yet handle all standard TO benchmarks. Geometric multigrid preconditioning would substantially reduce the per-step CG count and is the most impactful single improvement available; it is the primary direction of ongoing work. Benchmark scope. The quantitative main-text evidence now spans cantilever scaling, a bridge hard-problem stress test, and a torsion stress test in which 108–110 of 120 SIMP steps hit the current 1,000-iteration CG cap. The paper therefore still does not establish comparable behavior on the broader family of standard 3D TO benchmarks, and MBB remains an unresolved cap-limited case under the current solver configuration. A companion energy measurement workflow polls GPU power draw via the standard NVIDIA command-line telemetry utility at 100 ms intervals and integrates the trace to yield per-run energy in Joules; the workflow records measured energy alongside the TDP × twall upper bound. These energy rows come from a separate instrumented workflow and should be compared only within the matched FP64/fused runs of that workflow, not against the main timing table row-by-row. The 100 ms cadence cannot resolve sub-100 ms power transients, so the reported Joule values should be read as board-level energy estimates rather than as cycle-accurate power integrals. Only the 216 k and 1 M cantilever cases were instrumented in that workflow. Measured board-power traces (SIMP120 cantilever, RTX 4090) give 0.648 Wh for the fused 216 k run (mean 131.5 W, 29% of rated TDP) versus 2.098 Wh for the FP64 baseline—a 3.24× energy reduction. At 1 M elements the fused kernel uses 2.979 Wh (mean 231.7 W) versus 14.670 Wh for FP64, a 4.92× energy reduction. The per-path TDP × twall upper bounds overestimate actual measured energy by 1.9–3.4× across the four tested configurations, with the fused 216 k bound (2.22 Wh vs. actual 0.65 Wh; 3.4×) being the largest and the fp64 1 M bound (28.33 Wh vs. actual 14.67 Wh; 1.9×) the smallest, because the GPU operates at only 29–52% of its rated 450 W TDP across the instrumented FP64 and fused runs. For the 358 s fused 2 M run the TDP upper bound is ≈ 45 Wh, but extrapolating the measured 231.7 W mean power gives ≈ 23 Wh. BF16 arithmetic inapplicable to the CG solve. As demonstrated in Section 4.7, BF16 arithmetic stagnates in the CG linear solve due to the high condition numbers inherent to SIMP-penalized stiffness matrices. The BF16 WMMA kernel is therefore benchmarked for permatvec throughput only, not as a usable component of a production SIMP solver. Integration as a multigrid smoother—where the coarse-grid correction bounds the effective κ seen by the smoother—is identified as the viable path to extracting BF16 throughput benefit. Single-GPU scope. The present solver operates on the tested single RTX 4090 consumer GPU. The largest reported stress test is the 8 M-element FEA-only solve, which uses 10.56 GB on the RTX 4090. Full SIMP runs require additional storage for density fields, sensitivities, and filter workspace, so practical single-GPU SIMP limits on the 24 GB card are lower than a simple FEA-only extrapolation would suggest; larger problems ultimately require either a larger single GPU or a distributed-memory multi-GPU approach. Extending the fused kernel to multi-GPU via NCCL-based halo exchanges and domain decomposition is left for future work.
5.5
Path to Geometric Multigrid Integration
Geometric multigrid (GMG) is the highest-priority solver follow-on to the present operator paper. The present fused kernel already provides the fine-grid smoother kernel: each V-cycle pre/post-smoothing step is a fixed number of matvec-plus-update operations with the same gather–GEMM–scatter structure. The additional components required are: (i) an inter-grid restriction operator (R: inject or full-weighting) that maps fine-grid residuals to the coarser grid, (ii) a coarse-grid prolongation operator (P = R⊤ ), and (iii) a coarse-grid direct solver (cuSolver or cuSPARSE for the coarsest level). For the Cartesian structured mesh used in SIMP TO, all three components have closed-form stencil representations and can be implemented as 31
additional CuPy runtime-compiled kernel launchers without leaving the Python ecosystem. Following the classical multigrid treatments of Briggs et al. and Trottenberg et al. [62, 63] together with the TO-specific guidance of Peetz and Elbanna [52], a hybrid GMG–AMG strategy is planned: GMG for the first 40–60 SIMP iterations (when the topology is diffuse and the hierarchy is geometrically regular), transitioning to AMG via PyAMG after the topology has condensed (when AMG’s robustness to irregular connectivity becomes advantageous). This strategy is intended to substantially reduce the CG workload per SIMP step. Under the simple assumption that the current 2 M fused run is dominated by CG work, that would move the end-to-end wall time materially lower, into a much more favorable regime than the present Jacobi-PCG path. Additionally, once GMG bounds the effective condition number at each level below the BF16 IR threshold, the BF16 WMMA kernel becomes viable as the V-cycle smoother on the fine grid, potentially recovering the theoretical tensor-core throughput advantage in a production setting. This two-path roadmap—GMG for preconditioner improvement and BF16 as a smoother within the resulting bounded-κ subproblem—defines the next phase of this research program.
6
Conclusion
This paper presented four contributions for matrix-free 3D SIMP topology optimization on the tested RTX 4090 consumer GPU. First, we introduced a fused gather–GEMM–scatter CUDA kernel implemented through CuPy’s runtime kernel-compilation interface, preserving a Python-native workflow while eliminating the intermediate per-element DRAM round-trips of the three-stage baseline. In the reported study, that fused path achieves 6.0–6.6× per-matvec speedup (synthetic hot-path microbenchmark) and 8.9–13.8× in isolated CUDA-event measurements on the actual operator (Table 3), the latter also reflecting the Python/CuPy dispatch overhead avoided by collapsing three API calls into one launch; end-to-end SIMP-120 wall-time speedup is problem-sizedependent: 4.6–7.3× on cantilever and 4.4× on the 499,125-element torsion benchmark. Against the same-precision FP32 three-stage path, the fused solver yields 2.3–4.6× on cantilever and 2.8× on torsion. Measured board-power traces show that the fused path also delivers 3.2–4.9× energy reduction relative to matched FP64 instrumented runs (0.65–2.98 Wh vs. 2.10–14.67 Wh at 216 k and 1 M elements), with the GPU operating at 29–52% of its 450 W rated TDP across those instrumented runs. Second, we reported end-to-end SIMP scaling on the cantilever benchmark together with a harder torsion stress test and a bridge hard-problem stress test, making explicit that the paper’s SIMP tables use selected valid iterates rather than last-iterate metrics and that the large FEA-only points are capped stress-test solves rather than fully converged ones. The bridge family is consistent with the same FP32 three-stage gain persisting on a distributed-load case, although that bridge study is not a direct fused-kernel benchmark. MBB remains the unresolved cap-limited case within the same hard-problem family. Third, we implemented a BF16 WMMA tensor-core variant of the fused kernel and showed that the reported 512 k profiling study includes a 14.3× BF16-versus-FP64 GEMM proxy timing from a separate PyTorch BF16 GEMM benchmark, even though the full fused-BF16 matvec remains bounded by the gather/scatter share. Fourth, we characterized the mixed-precision failure mode of BF16 in the present Jacobipreconditioned CG setting. Direct power-iteration estimates of κ(K) reported in Section 4.8 give κ ≈ 6.1 × 105 at 64 k, ≈ 1.3 × 106 at 216 k, and ≈ 2.3 × 106 at 512 k elements, placing εBF16 · κ ≈ 2.4 × 103 –9.1 × 103 across the three tested sizes—more than an order of magnitude above the convergence threshold. These direct measurements show why BF16-IR stagnates: the systems are firmly in the non-convergence zone, not at the boundary (Section 4.8). That result points toward BF16 as a more plausible multigrid smoother (where coarse-grid correction 32
bounds the effective κ below 256) than as a drop-in CG inner solve. Taken together, these results improve on the three-stage Python/CuPy FP64 baseline used as the reference in this study while also clarifying the present scope. The fused FP32 path is directly supported on the reported cantilever and torsion studies; the bridge evidence in this paper is currently FP32/FP64 three-stage only. The current Jacobi-preconditioned solver still does not establish comparable performance on harder benchmark families. The implementation is intended to serve as the foundation for a separate geometric-multigrid solver study needed to close the remaining iteration-count gap.
CRediT authorship contribution statement Shaoliang Yang: Conceptualization, Methodology, Software, Investigation, Formal analysis, Visualization, Writing – original draft. Jun Wang: Supervision, Conceptualization, Methodology, Writing – review & editing, Funding acquisition. Yunsheng Wang: Validation, Investigation, Writing – review & editing.
Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Data availability The full code-and-experiments repository for the present implementation will be made public at https://github.com/nbbllxx0/Fused-Gather-GEMM-Scatter-Kernels once the arXiv version is online. Until that release, the present preprint records the methodological details, measurement conventions, and reporting rules needed to interpret the reported experiments.
A
Reproducibility Notes
This appendix records the reader-facing settings needed to interpret the reported experiments without exposing repository-internal bookkeeping. Software and hardware context. All GPU experiments for the present implementation use the single RTX 4090 configuration stated in Section 4.1. The local PyTopo3D comparison in Section 4.3 is a separate CPU baseline rerun and is not part of that GPU solver path. The reported manuscript values were generated under Microsoft Windows 10.0.26200.8037 with NVIDIA driver 595.71, Python 3.11, CuPy 13.6.0, PyTorch 2.5.1 built against CUDA 12.1, NumPy 2.2.6, SciPy 1.15.3, Matplotlib 3.10.7, scikit-image 0.25.2, Pandas 2.3.3, and PyVista 0.46.3; the CuPy runtime reports CUDA runtime version 12090. The GPU used the default 450 W power limit. This CUDA 12.1/12.9 split reflects PyTorch’s bundled CUDA 12.1 user-space runtime alongside the installed CUDA 12.9 runtime used by CuPy; that mixed configuration is supported by the installed NVIDIA driver on this host. For the contextual PyTopo3D comparison, the deposited CPU rerun log records 20 CPU cores and 69.50 GB available RAM on the same Windows host, and the accompanying environment manifest records PyTopo3D 0.1.0, PyPardiso 0.4.7, and default host scheduling with no explicit thread pinning; Figure 7 should therefore be read as contextual timing evidence rather than as a version-controlled benchmark. The synthetic hot-path microbenchmark uses the adaptive repeat rule max(100, min(2000, int(5e7 // n_elem))) to target about 50 ms of total measurement per size. What is and is not public at the preprint stage. The full code-and-experiments repository for the present implementation will be made public at https://github.com/nbbllxx0/Fused33
Gather-GEMM-Scatter-Kernels once the arXiv version is online. At preprint time, readers should therefore treat this manuscript and Appendix B as the authoritative reader-facing record of the experimental settings, measurement definitions, and reporting conventions used for the reported results. Key interpretation rules. Table 7 reports single cold-start linear solves at uniform density ρ = 0.5 and p = 3; they are not full SIMP runs. Table 6 mixes repeated cold-start FEA means (five timed calls in the benchmark script’s default mode) with single representative SIMP-60 runs. The 2 M, 4.9 M, and 8 M FEA-only points are capped-at-1,000-iteration stress-test solves rather than fully converged linear solves. The FEA scaling-ladder summaries average three timed solves for 216 k, 512 k, and 1 M elements, and two timed solves for 2 M, 4.9 M, and 8 M elements. The SIMP scaling rows in Table 4 and Figure 5 are single representative SIMP-120 runs from the dedicated scaling workflow, one deposited row per (size, path) pair; repeat variability is reported separately in Table 8. The condition-number workflow uses fixed seeds (42 for power iteration and 123 for inverse iteration) and inner SciPy CG solves at 10−8 relative and absolute tolerance with a 2,000-iteration cap. The energy numbers come from a separate instrumented workflow and should be compared only within that workflow’s matched FP64/fused runs. That workflow polls board power through nvidia-smi at an intended 100 ms interval, starts sampling immediately before the SIMP call, stops immediately after solver return, and computes energy by trapezoid integration over the recorded timestamped samples. This cadence cannot resolve sub-100 ms power transients, so the resulting Joule values should be read as board-level energy estimates rather than as cycle-accurate power integrals.
Figure 14: Centered-patch cantilever topology from a separate fused-FP32 SIMP rerun outside the quantitative benchmark set. The selected design field uses 160 × 80 × 80 = 1,024,000 elements on a 2 × 1 × 1 domain with Vf = 0.10, filter radius rmin = 3.0, and 100 design iterations. The rendered surface is a smoothed marching-cubes isosurface at ρ = 0.5 from a single fixed view (80-step Laplacian smoothing, relaxation factor 0.08); the red face indicates the clamped boundary and the green marker indicates the centered loaded-face patch on the right face. Selected compliance = 1.830, selected grayness = 1.35×10−5 , and the selected design comes from iteration 70 of the 100-step run.
34
B
Metric Conventions and Artifact Provenance
This appendix states the reporting conventions and the local workflows and result records used while preparing the main tables and figures in the preprint. At preprint time, it serves as the reader-facing provenance record supporting the main paper. Selected versus final values. The main paper reports selected compliance values for the SIMP studies. Here, “selected” means the lowest-compliance iterate among the iterates that pass the run’s validity checks. In the present implementation, those checks require p ≥ 3.0 and grayness g < 0.25. These selected values are not guaranteed to be identical to the last iterate of the reported continuation schedule (120 steps in the main cantilever/torsion tables and 60 steps in Table 6). Figure interpretation. Figures 8, 9, 11, and 12 render selected design fields as smoothed marching-cubes isosurfaces at ρ = 0.5, using the render script’s fixed-view layouts and 80-step Laplacian smoothing (niter = 80, relaxation factor 0.08). Figure 12 is a qualitative companion and is not the source of Table 6. Figure 7 uses completed 64 k and 216 k PyTopo3D rows. Scope of the validation studies. The condition-number tabulation in Section 4.8 covers only uniform-density states (ρ = 0.5) at p ∈ {3.0, 5.0} for 64 k, 216 k, and 512 k elements. A separate uniform-density BF16 extension at p = 4.5 supports the brief remark about higher-penalization stagnation in Section 4.8, but it is not part of that condition-number table. The selected-iterate high-cap validation in Section 4.9 covers 216 k and 512 k for the FP64 and fused paths. The 2 M and 4.9 M cantilever scaling rows are summary-only representative runs and are not part of that high-cap validation set. Repository release status. The full code-and-experiments repository for the present implementation will be made public at https://github.com/nbbllxx0/Fused-Gather-GEMM-ScatterKernels once the arXiv version is online. The list below is therefore the preprint-side provenance record of the local scripts and result files used during manuscript preparation. Preprint workflow/result mapping. • Profiling table/figure. Generators: the synthetic hot-path profiling workflow and the paper plotting workflow. Local result records: the profiling summary CSV and structured JSON export. Notes: synthetic hot-path microbenchmark matching operator tensor shapes but using seeded index-pattern probes rather than structured cantilever connectivity. • Cantilever FEA scaling figure. Generators: the cantilever scaling workflow and the paper plotting workflow. Local result records: the mid-range, large, and 8 M FEA scaling summary CSVs. Notes: cold-start FEA rows at uniform density ρ = 0.5; the 2 M, 4.9 M, and 8 M rows are cap-limited stresstest solves. • Cantilever SIMP scaling table/figures. Generators: the cantilever scaling workflow and the paper plotting workflow. Local result records: the mid-range, 1 M, and 2 M/4.9 M SIMP scaling summary CSVs. Notes: representative single-run SIMP-120 scaling rows; reported compliance is the selected best-valid compliance, not the last iterate. • External CPU comparison figure. Generators: the external CPU rerun comparison workflow and the paper plotting workflow. Local result records: the PyTopo3D-versus-GPU comparison CSV, the deposited PyTopo3D CPU log, the dedicated 64 k fused rerun CSV, and the mid-range SIMP scaling CSV. Notes: contextual timing comparison only; the 64 k fused bar comes from a separate deposited 64 k rerun, and the PyTopo3D record is paired with the environment manifest that records PyTopo3D 0.1.0, PyPardiso 0.4.7, and default host scheduling with no explicit thread pinning. • Torsion table/trajectory figure. Generators: the torsion benchmark workflow and the paper plotting workflow. Local result records: the torsion summary CSV and the FP64, FP32, and fused iteration-history JSON files. Notes: 108–110 of the 120 SIMP steps hit the current 1,000-iteration CG cap across the three deposited histories. • MBB/bridge hard-problem table. Generator: the hard-problem benchmark workflow. Local result record: the hard-problem summary CSV. Notes: the MBB and bridge rows cover FP64 and FP32 three-stage paths only; fused is not benchmarked in this table.
35
• BF16 convergence table. Generator: the BF16 iterative-refinement smoke-test workflow. Local result records: the BF16 smoke-test CSV and JSON outputs. Notes: single cold-start linear solves at uniform density ρ = 0.5 and p = 3; not full SIMP runs. • Condition-number study. Generator: the condition-number estimation workflow. Local result records: the condition-number CSV and JSON outputs. Notes: uniform-density states only (ρ = 0.5; p ∈ {3, 5}); late-SIMP conditioning statements remain explicitly inferential. • Repeat, determinism, high-cap, and energy studies. Generators: the repeat-study, determinism, high-cap validation, and energy benchmark workflows. Local result records: the repeatstudy CSV/JSON pair, determinism CSV/JSON pair, fully converged CSV/JSON pair, and energy CSV/JSON pair. Notes: supporting studies for variability, atomic-scatter determinism, selectediterate validity, and power-trace energy. • Residual figure. Generators: the residual-capture workflow and the paper plotting workflow. Local result record: the CG residual-history JSON file. Notes: representative cold-start cantilever solve at uniform density ρ = 0.5. • Qualitative render panels. Generator: the topology-render workflow. Local result records: the selected-density arrays and metadata sidecars for the 216 k cantilever, 1 M cantilever, 499,125-element torsion, 216 k bridge, and 1 M low-volume-fraction corner-load examples. Notes: all qualitative panels use selected-design density fields rendered as smoothed marching-cubes isosurfaces at ρ = 0.5. The metadata sidecars record the selected-run/source metadata, while the fixed camera layouts and panel compositions are defined in the render script.
References [1] Martin P. Bendsøe. Optimal shape design as a material distribution problem. Structural Optimization, 1(4):193–202, 1989. doi: 10.1007/BF01650949. [2] George I. N. Rozvany, Ming Zhou, and Tom Birker. Generalized shape optimization without homogenization. Structural Optimization, 4(3–4):250–252, 1992. doi: 10.1007/BF01742754. [3] Martin P. Bendsøe and Ole Sigmund. Topology Optimization: Theory, Methods, and Applications. Springer Berlin Heidelberg, 2004. doi: 10.1007/978-3-662-05086-6. [4] Niels Aage, Erik Andreassen, and Boyan S. Lazarov. Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework. Structural and Multidisciplinary Optimization, 51(3):565–572, 2015. doi: 10.1007/s00158-014-1157-0. [5] Niels Aage, Erik Andreassen, Boyan S. Lazarov, and Ole Sigmund. Giga-voxel computational morphogenesis for structural design. Nature, 550(7674):84–86, 2017. doi: 10.1038/nature23911. [6] Eddie Wadbro and Martin Berggren. Megapixel topology optimization on a graphics processing unit. SIAM Review, 51(4):707–721, 2009. doi: 10.1137/070699822. [7] Vivien J. Challis, Anthony P. Roberts, and Joseph F. Grotowski. High resolution topology optimization using graphics processing units (GPUs). Structural and Multidisciplinary Optimization, 49(2):315–325, 2014. doi: 10.1007/s00158-013-0980-z. [8] Jesús Martínez-Frutos and David Herrero-Pérez. Large-scale robust topology optimization using multi-GPU systems. Computer Methods in Applied Mechanics and Engineering, 311: 393–414, 2016. doi: 10.1016/j.cma.2016.08.016. [9] David Herrero-Pérez and Pedro J. Martínez Castejón. Multi-GPU acceleration of largescale density-based topology optimization. Advances in Engineering Software, 157–158: 103006, 2021. doi: 10.1016/j.advengsoft.2021.103006. 36
[10] Jiangnan Hou, Jiajie Li, Shengfeng Zhu, Xindi Hu, and Zeyang Yu. Parallel computing on GPU with CuPy and vectorized SpMV for large-scale topology optimization. Finite Elements in Analysis and Design, 250:104388, 2025. doi: 10.1016/j.finel.2025.104388. [11] Tianyuan Qi, Junpeng Zhao, and Chunjie Wang. An efficient GPU solver for 3d topology optimization of continuous fiber-reinforced composite structures. Computer Methods in Applied Mechanics and Engineering, 435:117675, 2025. doi: 10.1016/j.cma.2024.117675. [12] Erik A. Träff, Anton Rydahl, Sven Karlsson, Ole Sigmund, and Niels Aage. Simple and efficient GPU accelerated topology optimisation: Codes and applications. Computer Methods in Applied Mechanics and Engineering, 410:116043, 2023. doi: 10.1016/j.cma.2023.116043. [13] Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohei Hido, and Crissman Loomis. CuPy: A NumPy-compatible library for NVIDIA GPU calculations. In ML Systems Workshop (LearningSys) at NIPS, 2017. URL https://tech.preferred.jp/en/publications/cupya-numpy-compatible-library-for-nvidia-gpu-calculations/. [14] Junpeng Wang, Niels Aage, Jun Wu, Ole Sigmund, and Rüdiger Westermann. Efficient large-scale 3D topology optimization with matrix-free MATLAB code. Structural and Multidisciplinary Optimization, 68(9):174, 2025. doi: 10.1007/s00158-025-04127-3. [15] Samuel Williams, Andrew Waterman, and David Patterson. Roofline: an insightful visual performance model for multicore architectures. Communications of the ACM, 52(4):65–76, 2009. doi: 10.1145/1498765.1498785. [16] NVIDIA Corporation. NVIDIA Ada GPU Architecture. https://images.nvidia.com/aemdam/Solutions/Data-Center/l4/nvidia-ada-gpu-architecture-whitepaper-V2.02.pdf, 2023. Whitepaper, version 2.02; Appendix A includes GeForce RTX 4090 peak BF16, FP32, and memory-bandwidth specifications. [17] Stefano Markidis, Steven Wei Der Chien, Erwin Laure, Ivy Bo Peng, and Jeffrey S. Vetter. NVIDIA tensor core programmability, performance & precision. In 2018 IEEE Int. Parallel and Distributed Processing Symp. Workshops (IPDPSW), pages 522–531, 2018. doi: 10. 1109/IPDPSW.2018.00091. [18] Massimiliano Fasi, Nicholas J. Higham, Mantas Mikaitis, and Srikara Pranesh. Numerical behavior of NVIDIA tensor cores. PeerJ Computer Science, 7:e330, 2021. doi: 10.7717/ peerj-cs.330. [19] Azzam Haidar, Stanimire Tomov, Jack Dongarra, and Nicholas J. Higham. Harnessing GPU tensor cores for fast FP16 arithmetic to speed up mixed-precision iterative refinement solvers. In SC18: Int. Conf. for High Performance Computing, Networking, Storage and Analysis, pages 603–613, 2018. doi: 10.1109/SC.2018.00050. [20] Erin Carson and Nicholas J. Higham. Accelerating the solution of linear systems by iterative refinement in three precisions. SIAM Journal on Scientific Computing, 40(2):A817–A847, 2018. doi: 10.1137/17M1140819. [21] Nicholas J. Higham and Theo Mary. Mixed precision algorithms in numerical linear algebra. Acta Numerica, 31:347–414, 2022. doi: 10.1017/S0962492922000022. [22] Greg Henry, Ping Tak Peter Tang, and Alexander Heinecke. Leveraging the BFloat16 AI datatype for higher-precision computations. In Proc. 26th IEEE Symposium on Computer Arithmetic (ARITH), pages 69–76, 2019. doi: 10.1109/ARITH.2019.00019.
37
[23] Thomas Borrvall and Joakim Petersson. Large-scale topology optimization in 3D using parallel computing. Computer Methods in Applied Mechanics and Engineering, 190(46– 47):6201–6229, 2001. doi: 10.1016/S0045-7825(01)00216-X. [24] Kai Liu and Andrés Tovar. An efficient 3D topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 50(6):1175–1196, 2014. doi: 10.1007/s00158014-1107-x. [25] Ole Sigmund. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 21(2):120–127, 2001. doi: 10.1007/s001580050176. [26] Erik Andreassen, Anders Clausen, Mattias Schevenels, Boyan S. Lazarov, and Ole Sigmund. Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization, 43(1):1–16, 2011. doi: 10.1007/s00158-010-0594-7. [27] Thomas J. R. Hughes, Itzhak Levit, and James Winget. An element-by-element solution algorithm for problems of structural and solid mechanics. Computer Methods in Applied Mechanics and Engineering, 36(2):241–254, 1983. doi: 10.1016/0045-7825(83)90115-9. [28] Graham F. Carey and Bo-nan Jiang. Element-by-element linear and nonlinear solution schemes. Communications in Applied Numerical Methods, 2(2):145–153, 1986. doi: 10. 1002/cnm.1630020205. [29] Stephan Schmidt and Volker Schulz. A 2589 line topology optimization code written for the graphics card. Computing and Visualization in Science, 14(6):249–256, 2011. doi: 10.1007/s00791-012-0180-1. [30] Jun Wu, Christian Dick, and Rüdiger Westermann. A system for high-resolution topology optimization. IEEE Transactions on Visualization and Computer Graphics, 22(3):1195– 1208, 2016. doi: 10.1109/TVCG.2015.2502588. [31] Martin Kronbichler and Katharina Kormann. A generic interface for parallel cell-based finite element operator application. Computers & Fluids, 63:135–147, 2012. doi: 10.1016/ j.compfluid.2012.04.012. [32] Martin Kronbichler and Katharina Kormann. Fast matrix-free evaluation of discontinuous Galerkin finite element operators. ACM Transactions on Mathematical Software, 45(3): 29:1–29:40, 2019. doi: 10.1145/3325864. [33] Tzanio Kolev, Paul Fischer, Misun Min, Jack Dongarra, Jed Brown, Veselin Dobrev, Tim Warburton, Stanimire Tomov, Mark S. Shephard, Ahmad Abdelfattah, Valeria Barra, Natalie Beams, Jean-Sylvain Camier, Noel Chalmers, Yohann Dudouit, Ali Karakus, Ian Karlin, Stefan Kerkemeier, Yu-Hsiang Lan, David Medina, Elia Merzari, Aleksandr Obabko, Will Pazner, Thilina Rathnayake, Cameron W. Smith, Lukas Spies, Kasia Swirydowicz, Jeremy Thompson, Ananias Tomboulides, and Vladimir Tomov. Efficient exascale discretizations: high-order finite element methods. International Journal of High Performance Computing Applications, 35(6):527–552, 2021. doi: 10.1177/10943420211020803. [34] Julian Andrej, Nabil Atallah, Jan-Phillip Bäcker, Jean-Sylvain Camier, Dylan Copeland, Veselin Dobrev, Yohann Dudouit, Tobias Duswald, Brendan Keith, Dohyun Kim, Tzanio Kolev, Boyan Lazarov, Ketan Mittal, Will Pazner, Socratis Petrides, Syun’ichi Shiraiwa, Mark Stowell, and Vladimir Tomov. High-performance finite elements with MFEM. International Journal of High Performance Computing Applications, 38(5):447–467, 2024. doi: 10.1177/10943420241261981.
38
[35] Zijian Cao, Qiao Sun, Tiangong Zhang, and Huiyuan Li. Towards a higher roofline for matrix-vector multiplication in matrix-free HOSFEM, 2025. [36] Jiri Filipovič, Matouš Madzin, Jan Fousek, and Luděk Matyska. Optimizing CUDA code by kernel fusion: application on BLAS. Journal of Supercomputing, 71(10):3934–3957, 2015. doi: 10.1007/s11227-015-1483-z. [37] Mohamed Wahib and Naoya Maruyama. Scalable kernel fusion for memory-bound GPU applications. In SC’14: Int. Conf. for High Performance Computing, Networking, Storage and Analysis, pages 191–202, 2014. doi: 10.1109/SC.2014.21. [38] Mohamed Wahib and Naoya Maruyama. Automated GPU kernel transformations in largescale production stencil applications. In Proc. 24th Int. Symposium on High-Performance Parallel and Distributed Computing (HPDC), pages 259–270, 2015. doi: 10.1145/2749246. 2749255. [39] Jonathan Ragan-Kelley, Connelly Barnes, Andrew Adams, Sylvain Paris, Frédo Durand, and Saman Amarasinghe. Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines. In Proc. 34th ACM SIGPLAN Conf. on Programming Language Design and Implementation (PLDI), pages 519–530, 2013. doi: 10.1145/2499370.2462176. [40] Patrick Amestoy, Alfredo Buttari, Nicholas J. Higham, Jean-Yves L’Excellent, Theo Mary, and Bastien Vieublé. Five-precision GMRES-based iterative refinement. SIAM Journal on Matrix Analysis and Applications, 45(1):529–552, 2024. doi: 10.1137/23M1549079. [41] Dhiraj Kalamkar, Dheevatsa Mudigere, Naveen Mellempudi, Dipankar Das, Kunal Banerjee, Sasikanth Avancha, Dharma Teja Vooturi, Nataraj Jammalamadaka, Jianyu Huang, Hector Yuen, Jiyan Yang, Jongsoo Park, Alexander Heinecke, Evangelos Georganas, Sudarshan Srinivasan, Abhisek Kundu, Misha Smelyanskiy, Bharat Kaul, and Pradeep Dubey. A study of BFLOAT16 for deep learning training, 2019. [42] M. A. Clark, R. Babich, K. Barros, R. Brower, and C. Rebbi. Solving lattice QCD systems of equations using mixed precision solvers on GPUs. Computer Physics Communications, 181(9):1517–1528, 2010. doi: 10.1016/j.cpc.2010.05.002. [43] Steve F. McCormick, Joseph Benzaken, and Rasmus Tamstorf. Algebraic error analysis for mixed-precision multigrid solvers. SIAM Journal on Scientific Computing, 43(5):S392– S419, 2021. doi: 10.1137/20M1348571. [44] Dechuang Yang, Yuxuan Zhao, Yiduo Niu, Weile Jia, En Shao, Weifeng Liu, Guangming Tan, and Zhou Jin. Mille-feuille: a tile-grained mixed-precision single-kernel conjugate gradient solver on GPUs. In SC24: International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–16, 2024. doi: 10.1109/SC41406.2024.00064. [45] Abdul Dakkak, Cheng Li, Jinjun Xiong, Isaac Gelado, and Wen-mei Hwu. Accelerating reduction and scan using tensor core units. In Proc. ACM Int. Conf. on Supercomputing (ICS), pages 46–57, 2019. doi: 10.1145/3330345.3331057. [46] Hiroyuki Ootomo and Rio Yokota. Recovering single-precision accuracy from tensor cores while surpassing the FP32 theoretical peak performance. International Journal of High Performance Computing Applications, 36(4):475–491, 2022. doi: 10.1177/10943420221090256. [47] Hiroyuki Ootomo, Katsuhisa Ozaki, and Rio Yokota. DGEMM on integer matrix multiplication unit. International Journal of High Performance Computing Applications, 38(4): 297–313, 2024. doi: 10.1177/10943420241239588. 39
[48] Yuetao Chen, Kun Li, Yuhao Wang, Donglin Bai, Lei Wang, Lingxiao Ma, Liang Yuan, Yunquan Zhang, Ting Cao, and Mao Yang. ConvStencil: Transform stencil computation to matrix multiplication on tensor cores. In Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming, pages 333–347, 2024. doi: 10.1145/3627535.3638476. [49] Yiwei Zhang, Kun Li, Liang Yuan, Jiawen Cheng, Yunquan Zhang, Ting Cao, and Mao Yang. LoRAStencil: Low-rank adaptation of stencil computation on tensor cores. In SC24: International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–17, 2024. doi: 10.1109/SC41406.2024.00059. [50] Cu Cui. Acceleration of tensor-product operations with tensor cores. ACM Transactions on Parallel Computing, 11(4):15:1–15:24, 2024. doi: 10.1145/3695466. [51] Oded Amir, Niels Aage, and Boyan S. Lazarov. On multigrid-CG for efficient topology optimization. Structural and Multidisciplinary Optimization, 49(5):815–829, 2014. doi: 10.1007/s00158-013-1015-5. [52] Darin Peetz and Ahmed Elbanna. On the use of multigrid preconditioners for topology optimization. Structural and Multidisciplinary Optimization, 63(2):835–853, 2021. doi: 10.1007/s00158-020-02750-w. [53] Blaise Bourdin. Filters in topology optimization. International Journal for Numerical Methods in Engineering, 50(9):2143–2158, 2001. doi: 10.1002/nme.116. [54] Fengwen Wang, Boyan Stefanov Lazarov, and Ole Sigmund. On projection methods, convergence and robust formulations in topology optimization. Struct. Multidiscipl. Optim., 43(6):767–784, 2011. doi: 10.1007/s00158-010-0602-y. [55] Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6):409–436, 1952. doi: 10.6028/jres.049.044. [56] Yousef Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2nd edition, 2003. doi: 10.1137/1.9780898718003. [57] Nan Ding and Samuel Williams. An instruction roofline model for GPUs. In 2019 IEEE/ACM Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS), pages 7–18, 2019. doi: 10.1109/PMBS49563.2019.00007. [58] Jihoon Kim and Namwoo Kang. PyTopo3D: A Python framework for 3D SIMP-based topology optimization, 2025. arXiv preprint. [59] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins University Press, 4th edition, 2013. ISBN 978-1-4214-0794-4. doi: 10.56021/9781421407944. [60] Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM, 2nd edition, 2002. doi: 10.1137/1.9780898718027. [61] NVIDIA Corporation. NVIDIA A100 Tensor Core GPU. https://www.nvidia.com/en-us/ data-center/a100/, 2020. Product page; specifications table lists 80 GB HBM2e and up to 2.039 TB/s bandwidth on the A100 80GB SXM configuration; accessed 2026-04-17. [62] William L. Briggs, Van Emden Henson, and Steve F. McCormick. A Multigrid Tutorial: Second Edition. SIAM, 2nd edition, 2000. doi: 10.1137/1.9780898719505. [63] Ulrich Trottenberg, Cornelius W. Oosterlee, and Anton Schüller. Multigrid. Academic Press, 2001. ISBN 978-0-12-701070-0. 40