CUTh-Solver: GPU-Accelerated Sparse Matrix Solver for High-Resolution Thermal Simulation of 3D ICs

Abstract: Coarse-grained thermal simulation tends to underestimate localized thermal issues, potentially missing critical hotspots. Accurate analysis, therefore, demands fine-grained information, which dramatically increases grid resolution and thus computational workload. Fortunately, the coefficient matrices are often sparse with regular sparsity patterns, offering optimization opportunities. However, existing general-purpose matrix solvers on GPUs rarely exploit these domain-specific properties, thereby encountering bottlenecks in data storage, memory access, parallelism, computational efficiency, and hardware utilization. Therefore, we propose CUTh-Solver, a co-designed GPU-accelerated Preconditioned Conjugate Gradient (PCG)-based sparse solver framework for Symmetric Positive Definite (SPD) systems arising from high-resolution steady-state and transient 3D IC thermal simulation. For data storage, CUTh-Solver condenses the Diagonal (DIA) storage format to remove redundancy. To optimize the memory access, CUTh-Solver employs diagonal-wise SpMV to achieve coalesced memory access. We further observe a critical conflict between parallelism and preconditioning quality and thus adopt a high-parallelism preconditioning strategy. To improve computational efficiency and hardware utilization, we employ an adaptive fine-grained mixed-precision strategy that leverages diverse floating-point units to avoid resource contention, enhancing throughput without compromising numerical stability. Experimental results show that CUTh-Solver achieves up to 25.8x speedup over GPU-accelerated COMSOL Multiphysics 6.4 and over 3x speedup over NVIDIA's native general-purpose libraries (AmgX, cuSPARSE, cuDSS). Ablation studies validate the individual contribution of each optimization. The code is available at: https://github.com/Chenghan-Wang/CUTh-Solver

• The paper introduces CUTh-Solver, a domain-specific GPU-accelerated PCG framework that leverages structured thermal matrices to achieve significant speedup over conventional solvers.
• It employs a novel Quad-Diag storage format and diagonal-wise, coalesced SpMV kernel to optimize GPU memory bandwidth and reduce redundancy in high-resolution thermal simulations.
• Adaptive mixed-precision and Chebyshev polynomial preconditioning strategies combine with optimized storage to deliver up to 103x acceleration while maintaining high simulation accuracy.

CUTh-Solver: A GPU-Accelerated Sparse Solver for High-Resolution Thermal Simulation of 3D ICs

Motivation and Problem Statement

Emerging 3D integrated circuits (3D ICs) necessitate high-resolution thermal simulation due to localized hotspots that coarse-grained simulators systematically underestimate. State-of-the-art tools require fine mesh sizes (down to 1 μm), resulting in extremely large sparse linear systems governing heat transfer. Solving these systems is bottlenecked by computational throughput, memory footprint, and GPU hardware utilization. Existing general-purpose GPU solvers, such as cuDSS, AmgX, and cuSPARSE, do not exploit the specific structure and thermal physics inherent to 3D ICs, especially the regular banded and SPD nature of the system matrix, leading to significant inefficiencies.

Figure 1: Fine-grained ($1\ \mu$m) power mapping identifies peak power densities up to $6.3\times$ higher than coarse-grained ($20\ \mu$m) mapping, underscoring the essentiality of high-resolution simulation.

Methodological Innovations

The paper introduces CUTh-Solver, a domain-specific, GPU-accelerated preconditioned conjugate gradient (PCG) framework optimized for high-resolution (millions of cells) steady-state and transient thermal analysis in 3D ICs. The architecture pivots on four key co-designed innovations:

1. Redundancy-Free Storage Format: CUTh-Solver proposes Quad-Diag, which leverages the SPD and seven-diagonal banded structure of thermal matrices to eliminate the inherent redundancy and padding in conventional DIA storage. Quad-Diag stores only the upper triangle, exploiting diagonal unalignment and data locality.
2. Diagonal-wise, Coalesced SpMV Kernel: By reordering sparse matrix-vector multiplication (SpMV) to be diagonal-aligned, CUTh-Solver achieves fully coalesced DRAM accesses, maximizing GPU memory bandwidth utilization and addressing the memory wall.
3. High-Parallelism Polynomial Preconditioning: The solver adopts a Chebyshev polynomial preconditioner instead of classic incomplete Cholesky or algebraic multigrid schemes. This preconditioner is comprised entirely of BLAS-2 routines, circumventing global reductions and serial dependencies, and coupling seamlessly with the Quad-Diag SpMV.
4. Adaptive Fine-Grained Mixed-Precision Strategy: Capitalizing on the error tolerance of thermal simulations, the solver dynamically assigns floating-point precision per phase and per component, exploiting the prevalence of lower-precision arithmetic units on recent consumer GPUs without sacrificing convergence or stability.

Figure 3: CUTh-Solver identifies and exploits four domain-specific opportunities—storage, memory access, parallelism, and mixed precision—missing in general-purpose solvers.

Numerical Schemes and Kernel Design

The linear systems stem from standard discretizations of the heat equation using a 7-point finite difference stencil, yielding an SPD banded system. By formulating SpMV as diagonal-wise (versus classical row-wise), CUTh-Solver enables contiguous memory operations for each diagonal vector, paramount for high arithmetic throughput on modern GPU architectures.

Figure 2: Structure of the 3D IC mesh, 7-point stencil, and corresponding banded SPD matrix with seven diagonals.

Figure 4: Schematic of the Quad-Diag storage format, diagonal-wise SpMV algorithm, and resulting coalesced memory transaction patterns.

Importantly, the Chebyshev preconditioning routine is composed nearly exclusively of SpMV and AXPY routines, allowing the full benefit of the SpMV kernel optimizations. Furthermore, the mixed-precision scheduling exploits the physical structure of the matrix: outer-plane and in-plane conductances are assigned different precisions commensurate with their numerical significance and sensitivity.

Experimental Results

The authors benchmark CUTh-Solver against contemporary solvers on structured, real-world thermal matrices from high-resolution 3D IC designs. Three levels of comparison are emphasized: (1) end-to-end simulation time and accuracy versus established tools (COMSOL, 3D-ICE, HotSpot), (2) time-to-solution versus native GPU libraries, and (3) ablation studies on storage and computational kernel design.

Comparison Against Established Thermal Simulators

CUTh-Solver achieves up to $25.8\times$ speedup over GPU-accelerated COMSOL Multiphysics 6.4 and resolves the largest (up to $512^2\times10$) 3D ICs without out-of-memory failures, consistently outperforming 3D-ICE and HotSpot. Maximum and mean absolute errors with respect to COMSOL are below $0.2^\circ$C and $0.06^\circ$C, respectively.

Figure 5: End-to-end runtime comparison illustrating the performance inversion of classical preconditioners (IC-PCG) on GPUs, highlighting parallelism bottlenecks.

Figure 6: Benchmarked floorplans exhibit significant thermal heterogeneity requiring fine-grained resolution (SPARC, EV6, and multi-core floorplans).

Figure 7: Thermal fields computed by CUTh-Solver, COMSOL, and other baselines agree closely, with CUTh-Solver's error field confirming sign-off fidelity.

Comparison Against General-Purpose GPU Solvers

Against NVIDIA's cuSPARSE, cuDSS, and AmgX, CUTh-Solver attains $52.1\sim70.9\times$ speedup in steady-state cases and $53.8\sim103.4\times$ in transient thermal simulation, while always maintaining infinity-norm errors $<3.3\times 10^{-4}$ with respect to reference solutions. The framework does not encounter OOM limitations, unlike direct solvers and baseline iterative solvers.

Kernel Ablation and Synergy

Quad-Diag and coalesced SpMV contribute an isolated $6.3\times$0 kernel-level and $6.3\times$1 wall-clock speedup over classical CSR+iSpMV implementations. Adaptive mixed-precision scheduling (with diagonal-wise, component-aware and stage-aware scheduling) delivers $6.3\times$2 acceleration over fixed-precision solvers, with FP16-only approaches invariably diverging and FP32 incurring significant drift. Polynomial (Chebyshev) preconditioning confers a greater than $6.3\times$3 synergistic acceleration with these optimizations, as the SpMV-bound preconditioning amplifies the impact of storage and memory layout improvements.

Implications and Future Directions

The CUTh-Solver framework demonstrates that leveraging matrix symmetry, banded structure, and physical model awareness unlocks unprecedented runtime performance for thermal analysis of 3D ICs. The research empirically demonstrates that classical CPU-centric preconditioning (IC, AMG) induces parallelism bottlenecks and is fundamentally suboptimal on modern GPUs. The adoption of polynomial, SpMV-centric preconditioning, combined with storage innovations, is key to realizing the throughput potential of current and next-generation simulators.

Practically, the work makes sign-off accuracy thermal simulation feasible for chip-scale (tens of millions of unknowns) 3D ICs within seconds, overcoming previous memory and computational limitations. Theoretically, it exposes the necessity of holistic co-design—algorithm, storage, precision management, and hardware—for effective scientific computing on GPU accelerators.

The primary limitation is applicability to highly regular, structured mesh matrices; extensions to fully unstructured matrices (e.g., arising in FEM with tetrahedral meshing) are non-trivial and require hybrid or alternative storage formats (e.g., hybrid DIA-ELL/COO).

Conclusion

CUTh-Solver establishes a new state of the art in GPU-based sparse linear solvers for high-resolution thermal simulation in 3D ICs. By harmonizing storage format, memory access, preconditioner design, and mixed-precision execution, the framework delivers order-of-magnitude performance improvements over both legacy and contemporary alternatives, without trading off numerical accuracy. The paradigm advanced in this work is extensible to other scientific and engineering simulation problems where structured, banded, SPD linear systems arise.