Matrix-Free 3D SIMP Topology Optimization with Fused Gather-GEMM-Scatter Kernels

The Thesis

Topology optimization — the computational process of finding the most efficient material layout for a structure, like a lightweight aircraft bracket or bridge truss — is notoriously slow because it requires solving huge systems of equations repeatedly. This paper attacks one specific bottleneck: the wasted memory bandwidth caused by running three separate GPU operations (gather, matrix multiply, scatter) in sequence, each reading and writing results to slow DRAM. The authors fuse all three into a single CUDA kernel, meaning data stays in fast on-chip memory for the entire operation. The result is a 4.6–7.3x end-to-end speedup on a single RTX 4090 consumer GPU for problems up to 4.9 million finite elements — with a 3–5x reduction in energy per solve. The catch is significant: this only works on problems that fit in 24 GB of GPU memory, and a lower-precision variant (BF16) fails numerically due to poor conditioning.

Catalyst

Consumer GPU memory has crossed the 24 GB threshold with cards like the RTX 4090, making it practical to run medium-scale 3D topology optimization problems entirely on-device without distributed computing. Simultaneously, JIT (just-in-time) compilation interfaces like CuPy's runtime compilation now let researchers write and deploy custom CUDA kernels without a full C++ toolchain, lowering the barrier to fused-kernel experimentation. These two changes together make single-GPU, fused-kernel topology optimization newly tractable.

What's New

Prior matrix-free topology optimization solvers — including implementations built on PETSc and GPU-accelerated finite element frameworks — separated the gather, matrix multiply (GEMM), and scatter operations into distinct kernel calls, each requiring a round-trip to DRAM (the GPU's main memory). This paper replaces those three round-trips with a single fused CUDA kernel that keeps intermediate results in the GPU's L2 cache and registers throughout. The claimed advantage over the same-precision (FP32) three-stage baseline is still 2.3–4.6x, meaning the fusion itself — not just the switch from FP64 to FP32 — contributes meaningfully to the speedup.

The Counter

The speedups reported here are real, but the comparison is carefully chosen. The headline 4.6–7.3x figure is against a mixed-precision baseline (FP32 fused vs. FP64 three-stage), which conflates two separate advantages — lower precision and better memory access — in a way that flatters the result. The apples-to-apples FP32-versus-FP32 comparison yields a more modest 2.3–4.6x, and that range has wide variance across problem sizes. More importantly, the paper solves a problem that the topology optimization community has already been attacking from multiple angles: GPU-accelerated solvers from academic groups and commercial vendors already exploit matrix-free methods and mixed precision. The question is whether this specific fused-kernel trick generalizes beyond the two benchmark geometries tested, and the paper provides no evidence that it does. Industrial users also rarely run on a single RTX 4090 — they run on A100 or H100 clusters where the memory hierarchy, cache sizes, and interconnects are completely different, and those benchmarks are absent. Finally, the failure of BF16 is a real setback: modern GPU hardware is optimized for Tensor Core operations at BF16 or FP8, and a method that cannot use those units safely will fall further behind as hardware evolves.

The Paper