arXiv CSopen access

FusionRCG: Orchestrating Recursive Computation Graphs across GPU Memory Hierarchies

FusionRCG: Orchestrating Recursive Computation Graphs across GPU Memory Hierarchies Xinran Wei†

Yihong Zhang†

Junshi Chen

arXiv:2605.10312v1 [physics.comp-ph] 11 May 2026

University of Science and Technology of China Zhongguancun Academy University of Science and Technology of China Hefei, Anhui, China Beijing, China Hefei, Anhui, China [email protected] [email protected] [email protected]

Fusong Ju∗

Wei Hu∗

Jinlong Yang

Zhongguancun Academy Beijing, China [email protected]

University of Science and Technology of China Hefei, Anhui, China [email protected]

University of Science and Technology of China Hefei, Anhui, China [email protected]

Huanhuan Xia∗ Zhongguancun Academy Beijing, China [email protected]

Abstract—Evaluating high-dimensional integrals via deep hierarchical recurrences is a dominant cost in quantum chemistry. While CPUs manage these efficiently, GPUs suffer a critical mismatch: limited per-thread memory is quickly overwhelmed by an explosion of simultaneously live intermediate variables. As recurrence scales, this forces massive data spilling to global memory, collapsing performance into a severe memory-bound regime. We present FusionRCG, a framework that jointly optimizes computation graph structure and GPU memory mapping. Exploiting the inherent topological flexibility of recurrence graphs, using electron repulsion integrals as an example, we contribute: (1) liveness-aware graph orchestration to minimize peak live intermediates; (2) algebraic dimensionality reduction via stepwise Cartesian-to-spherical fusion, shrinking intermediate footprints by up to 7.7×; and (3) an adaptive multi-tier kernel architecture routing graphs across the memory hierarchy. Evaluated on NVIDIA A100 GPUs, FusionRCG achieves up to 3.09× end-toend SCF speedup over GPU4PySCF and maintains 75% parallel efficiency at 64 GPUs, successfully rescuing these workloads from memory-bound limits. Index Terms—GPU, memory hierarchy, recursive computation graph scheduling, register pressure, electron repulsion integrals

I. I NTRODUCTION In modern large-scale quantum chemistry simulations, evaluating high-dimensional tensor integrals via deep hierarchical recurrences constitutes a critical computational bottleneck [1]– [3]. From an algorithmic and computer architecture perspective, such workloads can be naturally modeled as Recursive Computation Graphs (RCGs), where vertices represent intermediate quantities and edges encode recurrence dependencies. As the polynomial degree or the dimensionality of the † These authors contributed equally to this work.

Corresponding authors: Fusong Ju ([email protected]), Wei Hu ([email protected]), Leo Xia ([email protected]).

problem increases, the number of intermediate variables that must be simultaneously maintained within this computation graph undergoes a combinatorial explosion. While CPUs can efficiently manage this super-polynomial growth by utilizing sophisticated out-of-order execution engines and deep, multilevel cache hierarchies to absorb the exploding intermediate state [4], [5], porting such recursive workloads to modern GPUs exposes a severe architectural mismatch. GPU performance relies on massive thread-level parallelism, which strictly constrains per-thread on-chip storage [6]. When the working set of simultaneously live intermediates breaches the hard ceiling of physical registers, the compiler is forced to spill massive amounts of data to the high-latency, off-chip global memory [7]. This pure memory-hierarchy crisis drastically collapses hardware thread occupancy, submerging the execution in memory access latency and degrading an inherently compute-bound kernel into a memory-bound bottleneck [8]. A representative and highly demanding instantiation of this memory-bound RCG is the calculation of electron repulsion integrals (ERIs). ERIs are 4-index tensor integrals that describe particle interactions, serving as a foundational operator in many high-fidelity scientific simulations and typically consuming 80% of total runtime [1]–[3], [9]. Accelerating ERI calculation has therefore been a long-standing objective within the computational chemistry community, and existing GPU engines have achieved notable success [2], [9]. Most successful GPU implementations are based on Rys quadrature or McMurchie–Davidson (MD) formulations [9]–[11]. These formulations are well aligned with SIMT execution because they expose regular batched parallelism and maintain a comparatively controlled intermediate footprint.

tightly coupled dimensions relevant to GPU execution: • Computation Graph Structure and Execution Order. The HGP recurrence does not induce a unique computation graph. Different legal reduction paths generate structurally distinct graphs, and different topological orders expose different live ranges even for the same graph. Together, these factors determine the peak number of simultaneously live intermediates. • Intermediate State Size Optimization. The cost of HGP on GPUs is governed not only by arithmetic count, but also by the physical size of the intermediate state carried between recurrence stages and contraction. The basis and representation in which these intermediates are maintained directly affect both storage demand and data movement. • Memory Mapping across GPU Hierarchy. As angular momentum increases, the working set of the recurrence inevitably outgrows the smallest on-chip storage. Performance therefore depends on how intermediate state is placed across registers, shared memory, and global memory, rather than on arithmetic optimization alone. Based on these insights, we present FusionRCG, a framework that addresses the memory wall for hierarchical recursive computation graphs by systematically orchestrating their evaluation. By co-designing the graph structure, memory mapping, and algebraic boundaries, FusionRCG preserves the exact mathematical arithmetic of HGP while eliminating its historical GPU bottlenecks. Our contributions are threefold, directly addressing the dimensions above:

Fig. 1: Recursive computation graph structure of HGP and its associated GPU bottlenecks. The first recurrence phase generates a rapidly expanding dependency graph with a large peaklive frontier, while the second recurrence phase introduces long-lived accumulators; these phases are formalized later as VRR and HRR, respectively, in Section II-A. Their overlap within the same register file creates the fundamental memory bottleneck.

The Head-Gordon-Pople (HGP) algorithm, by contrast, represents a markedly different ERI formulation. HGP and its variants remain highly influential in CPU engines due to their favorable contraction-aware arithmetic structure and their substantially lower floating-point operation count than Rys quadrature and McMurchie–Davidson formulations [1], [12], [13]. However, executing HGP efficiently on GPUs is notoriously difficult, and mainstream adoption has been severely limited. This difficulty is fundamentally driven by HGP’s mathematical formulation: it evaluates an ERI through deep, chained recurrences over Cartesian intermediates. When mapped onto SIMT hardware, a target integral expands into a massive dependency graph where many intermediates are produced early but consumed much later. This forces an enormous working set to remain simultaneously live, triggering the exact memory-hierarchy crisis described above. Because of this severe register pressure, massive data spilling, and excessive global-memory traffic, many optimization efforts for HGP on GPUs inevitably hit a performance wall. Figure 1 illustrates this recursive computation graph structure and the memory bottlenecks during two HGP phases (formalized in Section II-A). Bridging the gap between HGP’s arithmetic advantages and its limited practical viability on GPUs requires more than accelerating a fixed kernel. The key insight of this work is that a strictly defined mathematical result does not imply a unique computation graph. This reveals a previously under-exploited optimization space in which HGP can be analyzed along three

• Liveness-aware graph orchestration for registercontained execution (§III-A). By jointly optimizing the graph topology and the evaluation schedule, FusionRCG drastically minimizes the peak working set of live intermediates. This structural optimization effectively eliminates register spilling for low-to-moderate degree recurrences, keeping the dense execution entirely on-chip. • Algebraic dimensionality reduction to minimize globalmemory traffic (§III-B). We step-wise fuse the standard Cartesian-to-pure-spherical shell transformation [5], [14] directly into the final stages of the recurrence. This shrinks the contraction-bound intermediate tensors and the associated buffer writes / atomicAdd traffic by up to 7.7× without introducing any mathematical approximation. • An adaptive multi-tier execution strategy across the GPU memory hierarchy (§III-C). We demonstrate that memory pressure arises from distinct sources as the problem scales. To address this, the framework automatically routes each computation graph to the appropriate memory level—emitting fundamentally different kernel architectures (register-only or shared-memory-buffered) based on computational complexity.

2

Evaluated on NVIDIA A100 GPUs against the state-ofthe-art GPU4PySCF [9], FusionRCG achieves up to 3.09× end-to-end SCF speedup on representative molecular systems with an average of 2.4× at the cc-pVQZ level, and scales to 64 GPUs with 75% parallel efficiency, successfully rescuing the execution from the memory-bound regime. Figure 2 provides an overview of the code generation and memory mapping pipeline.

which requires maintaining the architectural state of thousands of in-flight threads simultaneously. To support this, NVIDIA GPUs provide a massive but strictly partitioned register file. Each Streaming Multiprocessor (SM) contains 65,536 32-bit registers shared among all resident threads, imposing a hard per-thread architectural ceiling of 255 registers. When a compute kernel’s working set exceeds this 255register limit, the compiler is forced to spill excess intermediate values to per-thread local memory, which is physically backed by off-chip High Bandwidth Memory (HBM). These spills are devastating to performance because the bandwidth gap across the memory hierarchy is extreme. The SM register file delivers an aggregate bandwidth of ∼19 TB/s, whereas global HBM provides only ∼2 TB/s—a disparity of nearly 10,000× at the per-SM level. Between these two extremes, the shared memory/L1 cache (19 TB/s, 164 KB/SM) and L2 cache (5 TB/s, 40 MB) offer vital intermediate capacity-bandwidth trade-offs. FusionRCG utilizes this exact hierarchy, mapping intermediates strictly to registers (Tier 1) or buffering them through shared memory and global memory (Tier 2) based on the computation’s exact spatial footprint.

II. BACKGROUND AND M OTIVATION A. The HGP Recurrence for ERIs To understand the origin of the memory bottleneck, we first examine the underlying mathematical formulation. An electron repulsion integral over four Gaussian basis functions a, b, c, d is defined as: ZZ ϕa (r1 )ϕb (r1 ) ϕc (r2 )ϕd (r2 ) (ab|cd) = dr1 dr2 (1) |r1 − r2 | where ϕa is a Cartesian Gaussian with angular momentum a = (ax , ay , az ), and total angular momentum |a| = la . In the notation below, (ab|cd) denotes the target contracted ERI, whereas the bracketed quantities [· · · | · · · ](m) denote the HGP auxiliary integral family indexed by the Boys order m. The HGP algorithm [1] evaluates Eq. (1) through a two-phase hierarchical recurrence:

C. From Recurrence to Register Wall When the HGP mathematical formulation is mapped onto the SIMT execution model, it manifests as a dense computation graph whose size grows combinatorially with the angular-momentum quartet l = (la , lb , lc , ld ). Because each node in this graph represents a double-precision intermediate occupying two 32-bit registers, even modest quantum systems generate massive register demand. P For example, at l = (1, 1, 1, 1) ( li = 4), the graph contains over 700 nodes. If all intermediates were simultaneously live, the naive register requirement would exceed the 255-register hardware limit by nearly 6×. By l = (2, 2, 2, 2) P ( li = 8), the graph grows to thousands of nodes, surpassing the limit by more than 28×. Fortunately, the hardware does not need to store the entire graph simultaneously. A register is occupied only while the corresponding intermediate remains needed by a future consumer; once all uses of that value are complete, its register is reclaimed. Therefore, the actual bottleneck within the recurrence graph is the peak simultaneously live set. Let G = (V, E) denote the computation graph, where each vertex v ∈ V is an intermediate and each edge encodes a data dependency. Given a topological ordering π of G, a vertex v is live at step t if it has been computed but at least one of its consumers has not yet been scheduled. The peak liveness is defined as:

Phase 1: Vertical Recurrence Relation (VRR). Starting from Boys-function base values Fm (T ) ≡ [0 0|0 0](m) , VRR builds up angular momentum one index component at a time: [a + 1i 0|c 0](m) = (Pi − Ai ) [a 0|c 0](m) + (Wi − Pi ) [a 0|c 0](m+1)
ρ ai [a − 1i 0|c 0](m) − [· · · ](m+1) + 2ζ ζ ci + [a 0|c − 1i 0](m+1) (2) 2(ζ + η) where 1i is the unit vector along axis i, and ζ, η, ρ, P, W are Gaussian-product parameters. Each node depends on 2 to 5 parents. Crucially, the choice of which index i to increment is not strictly fixed, providing a structural degree of freedom that we will exploit. Phase 2: Horizontal Recurrence Relation (HRR). After contraction, HRR transfers angular momentum from center-pair products to individual centers via a simpler binary recurrence: (a b + 1i |c d) = (a + 1i b|c d) + (Ai − Bi ) (a b|c d) (3) Each HRR node depends on exactly 2 parents. While the operation count of this two-phase approach is highly efficient, executing its data dependencies on massively parallel hardware exposes severe physical constraints.

PeakLive(π) = max {v ∈ V | v live at step t} t

B. GPU Memory Hierarchy and the Register Wall

(4)

This PeakLive metric is highly sensitive to both the graph’s structural topology and its evaluation order—the two degrees of freedom that FusionRCG optimizes.

Modern GPUs derive their massive throughput from the Single Instruction Multiple Threads (SIMT) execution model,

3

Fig. 2: FusionRCG overview. Given an angular-momentum quartet, the generator first optimizes the Phase 1 VRR graph for low peak liveness, then selects a register-only Tier 1 backend or a buffered Tier 2 backend based on decomposed register demand, and finally fuses the standard Cartesian-to-pure-spherical shell transformation into the Phase 2 HRR stage before generating specialized CUDA kernels. Algorithm 1 Axis-priority rule (consumer-first, z → y → x). 1: for i = 0, 1, 2 do 2: if c[i] > 0 then 3: axis ← i + 3 4: else if a[i] > 0 then 5: axis ← i 6: end if 7: end for

Rtotal ≈ 2 PeakLiveVRR + Nacc

(2li + 1)

(6)

Ultimately, the register wall is not caused by isolated bottlenecks, but by the overlap of a massive short-lived VRR frontier colliding with a persistent accumulator bank within the same finite register file. Figure 1 visually summarizes this sharedpressure bottleneck and how the combined footprint dictates the feasible execution regime.

However, peak liveness captures only the transient portion of the on-chip footprint. In a full HGP kernel, these shortlived VRR intermediates must coexist with a second, longlived term: the Phase 2 accumulators. These accumulators are allocated before the primitive loops begin and must remain resident across all primitive combinations. For a shell quartet l = (l1 , l2 , l3 , l4 ), their count is: 4 Y

// then left side

superposition:

Fig. 3: Axis selection reshapes the VRR computation graph before scheduling. Starting from the same recurrence state, consumer-first (left) and producer-first (right) generate graphs with identical arithmetic but different frontier width. The consumer-first graph is narrower, so fewer intermediates must remain live simultaneously.

Nacc =

// right side first

III. F USION RCG The three components of FusionRCG operate at distinct but coupled levels. At the mathematical level, HGP exposes flexibility in recurrence path selection and intermediate representation. At the algorithmic level, FusionRCG converts this flexibility into graph construction, scheduling, and staged transformation rules during code generation. At the architectural level, these decisions reduce register pressure, shrink data movement, and determine whether execution remains registerresident or must be routed through higher memory tiers.

(5)

i=1

The total shared register demand is therefore dictated by their

4

Dependency Scheduling Algorithm (DSA). Given a fixed computation graph G = (V, E) from the axis-priority construction, we then optimize only the evaluation order on that fixed computation graph. This subproblem is closely related to classical minimum-register sequencing and register-pressureaware instruction scheduling on trees, computation graphs, and GPU scheduling regions [15]–[20]. In contrast to Figure 3, which changes the graph structure itself, DSA keeps the nodes and edges unchanged and seeks a topological ordering π minimizing PeakLive(π). DSA is therefore not a new scheduling formulation; it is a lightweight domain-specific heuristic tailored to our offline recurrence code generation. Our greedy heuristic maintains a ready set R and at each step selects: X
v ∗ = arg max out(w) − in(w) (7)

FusionRCG integrates three co-designed optimizations within a single code-generation framework: liveness-aware recurrence orchestration (§III-A), algebraic dimensionality reduction of the contraction-bound intermediate space (§III-B), and adaptive multi-tier mapping across the GPU memory hierarchy (§III-C). All decisions in this section are made at code-generation time—the generated GPU kernel contains a fixed, fully unrolled instruction sequence with no runtime branching. A. Liveness-Aware Recurrence Orchestration Mathematically, HGP admits multiple equivalent recurrence paths for the same target integral. Algorithmically, we exploit this freedom through graph construction and topological scheduling at code-generation time. Architecturally, the objective is to minimize the peak live frontier so that Phase 1 remains as register-resident as possible. The HGP recurrence (Eq. 2) increments one component of the multi-index (a, c) per step. With a = (ax , ay , az ) and c = (cx , cy , cz ), there are up to six valid reduction axes at each node: three on the left side (ax , ay , az ) and three on the right side (cx , cy , cz ). Choosing different axes does not change the total node count or the total arithmetic—but it fundamentally reshapes the graph topology. A “wide” expansion (reducing high-degree axes first) creates broad wavefronts where many intermediates are simultaneously live. A “deep” expansion (reducing loworder axes first) creates narrow channels that limit concurrent liveness. Figure 3 isolates this graph-construction effect: two axis-selection strategies on the same starting node produce computation graphs with identical node count but 1.5–2× different PeakLive. Figure 2(b) then isolates the execution strategy after the computation-graph structure has been fixed. This observation is the foundation of our approach: before optimizing the evaluation order of a fixed graph (a compiler problem), we optimize the graph structure itself (a domainspecific algebraic problem).

v∈R

w∈Reach(v)

The intuition behind this scoring function is that it captures the “register release potential” of scheduling v next. Nodes reachable from v with high fan-out (out ≫ in) are distributors—scheduling them early enables their consumers to execute promptly and release registers. Nodes with high fan-in (in ≫ out) are accumulators—deferring them avoids premature liveness inflation. DSA thus prioritizes paths that lead to rapid register turnover (Figure 2(b)). The algorithm runs in O(|V |2 ) time, acceptable for offline code generation: the full schedule for l = (2, 2, 2, 2) (|V | = 3,545) completes in seconds on one CPU core. B. Algebraic Dimensionality Reduction via Spherical Fusion Mathematically, the Phase 2 result can be represented either in the Cartesian shell basis or, after an exact shell transformation, in the pure spherical basis. Algorithmically, FusionRCG fuses this basis change into the recurrence boundary instead of materializing the full Cartesian tensor first. Architecturally, this reduces the size of contraction-bound intermediates and the associated global-memory traffic. Phase 2 naturally produces contraction-bound intermediates resolved in the Cartesian Gaussian shell basis. In most electronic-structure implementations, however, the final shell representation is the pure spherical harmonic Gaussian basis, so an additional Cartesian-to-spherical transformation is required [5], [14]. For a shell of angular momentum l, this transformation has the standard form X (l) (sph) (cart) Ilm = tijk,lm Iijk , (8)

Axis-Priority Rule. FusionRCG uses a deterministic axispriority rule that selects the reduction axis at each VRR node (Algorithm 1). This ordering is motivated by two observations. First, reducing ket-side indices before bra-side indices (consumerfirst) generates subgraphs whose outputs are consumed sooner by the downstream HRR phase, thereby shortening bra-side accumulator lifetimes. Second, within each side, the z → y → x traversal order exploits the fact that z-components have the smallest fan-out per recursion step in the Cartesian basis ordering, creating narrower wavefronts that minimize simultaneously live nodes. On l = (1, 1, 1, 1), this priority reduces PeakLive by 42% versus the worst axis ordering (brafirst, x → y → z), with zero change to total arithmetic or node count.

i+j+k=l (l)

where the coefficients tijk,lm depend only on the angular momentum and can therefore be pretabulated or generated analytically [5], [14]. Accordingly, each transformed index shrinks from ncart (l) = (l + 1)(l + 2)/2 Cartesian components to nsph (l) = 2l + 1 pure spherical components.

5

The compression ratio grows rapidly with l: at l = 4, each index shrinks from 15 to 9 components, yielding 154 /94 = 7.7× fewer contraction-bound intermediates for a four-index tensor. This exact shell transformation directly reduces Phase 2 node count, atomicAdd traffic to global memory, and downstream contraction work.

Fig. 4: Cartesian-to-spherical transformation and stepwise fusion. (a) At l = 2, the standard shell transformation maps 6 Cartesian Gaussian components to 5 pure spherical components per index via the fixed coefficient matrix t(l) . (b) FusionRCG applies this transformation one index at a time under a 4-bit mask (0000 → 1000 → 1100 → 1110 → 1111), shrinking the intermediate tensor progressively instead of materializing the full Cartesian intermediate tensor first.

Stepwise Fusion into the Recurrence. A naive implementation would first materialize the full Cartesian intermediate tensor and only then apply Eq. (8). FusionRCG instead fuses this standard shell transformation into Phase 2, applying it incrementally at HRR boundaries. The transformation state is tracked by a 4-bit mask, so each step retires a Cartesian index immediately after it has been consumed and replaces it with its smaller spherical counterpart. The result is a progressively shrinking intermediate tensor rather than a materialize-thentransform pipeline. The correctness of this rearrangement follows from the fact that HRR operators are linear and act on disjoint index pairs. Let the contracted VRR output be the four-center Cartesian tensor
Hαβγδ := χA,cart χB,cart χC,cart χD,cart , (9) α γ β δ

T(la ) ⊗T(lb ) acts only on (α, β), the two operations commute:

Rcd T(la ) ⊗T(lb ) = T(la ) ⊗T(lb ) Rcd . (14) Advancing the bra transformation ahead of the ket HRR therefore does not change the mathematical result. Finally, the ket indices are transformed: X (l ) (l ) (3) (4) Gµνλσ := Tλγc′ Tσδd′ Gµνγ ′ δ′ . (15) γ ′ δ′

where α, β, γ, δ index Cartesian components on centers A, B, C, D. For a shell of angular momentum l, the corresponding pure spherical basis functions satisfy X (l) cart ϕsph Tµα χα . (10) µ :=

Chaining these steps gives the fused operator emitted by FusionRCG,

Gfinal = T(lc ) ⊗T(ld ) Rcd T(la ) ⊗T(lb ) Rab H. (16) Using Eq. (14), this is algebraically identical to the conventional materialize-then-transform pipeline,
T(la ) ⊗T(lb ) ⊗T(lc ) ⊗T(ld ) Rcd Rab H. (17)

α (la )

(lb )

(lc )

(ld )

The conventional scheme applies T ⊗T ⊗T ⊗T after both HRR sweeps. FusionRCG interleaves these transformations with the bra–ket HRR sequence. For compact notation, the following derivation groups the two bra transformations and the two ket transformations; the implementation realizes the same algebra through the per-index mask shown in Figure 4. First, the bra-side HRR operator Rab (·; Oab ), with Oab = A − B, produces   (1) Gαβγδ := Rab (H; Oab ) αβγδ . (11) The bra indices can then be transformed immediately: X (2) (la ) (lb ) (1) Gµνγδ := Tµα Tνβ Gαβγδ ,

Stepwise spherical fusion is therefore an implementation rearrangement, not a numerical approximation: it preserves the HGP result while avoiding the full Cartesian footprint during Phase 2. C. Adaptive Multi-Tier Kernel Architecture Mathematically and algorithmically, the first two components reduce the transient live set and the size of Phase 2 intermediates, but they do not remove the strong dependence of total footprint on angular momentum. Architecturally, the remaining problem is therefore one of memory placement: different quartets must be mapped to different levels of the GPU hierarchy according to their decomposed register demand. While the liveness-aware orchestration (§III-A) and algebraic dimensionality reduction (§III-B) presented above significantly reduce on-chip memory pressure, the residual register footprint still varies dramatically across angularmomentum quartets. A natural design choice would be to adopt a single buffered execution pipeline for all cases, since buffering through shared memory and global memory can accommodate arbitrarily large working sets. However, such a

(12)

αβ

or compactly G(2) = T(la ) ⊗T(lb ) G(1) . This step shrinks the bra half of the tensor from ncart (la )ncart (lb ) to nsph (la )nsph (lb ) components before the ket-side sweep begins. The ket-side HRR operator Rcd (·; Ocd ), with Ocd = C−D, then acts only on the ket indices:   (3) Gµνγ ′ δ′ := Rcd (G(2) ; Ocd ) µνγ ′ δ′ . (13) Because Rcd is linear and touches only (γ, δ), while

6

Algorithm 2 Tier 1 Register-Only ERI Kernel (Optimized Stepwise Spherical Fusion)

uniform strategy introduces non-negligible fixed overheads— global buffer allocation, a layout-conversion kernel, and two separate kernel launches with synchronization barriers. For low-to-moderate degree quartets where total register spill remains below ∼1 KB, these fixed costs exceed the spill penalty itself, making a lightweight register-only path consistently faster (validated in §V-B). Conversely, for high-degree quartets whose persistent Phase 2 accumulator bank alone surpasses the hardware register budget, no scheduling improvement can fit the data on chip, and buffered execution becomes structurally mandatory. This asymmetry motivates our adaptive two-tier kernel architecture. The tier decision is derived from the decomposed register pressure model in Eq. (6): the same register file is jointly occupied by the transient VRR live frontier and the persistent accumulator bank (Figure 1), so the feasible execution regime depends on their combined footprint relative to the 255-register architectural limit. Concretely, this yields two hardware regimes:

1: Input: Shell quartet (a, b, c, d) with primitives {(αi , ci )},

density D 2: H[0 . . . Nhrr − 1] ← 0

▷ contraction accumulators (registers) 3: for each primitive quartet (αa , αb , αc , αd ) do 4: ζ, η, P, Q, T ← G AUSSIAN PAIRS(αa , αb , αc , αd ) 5: F0 , . . . , Fm ← B OYS(T, m) 6: V ← VRRπ∗ (F0 . . . Fm , ζ, η, P, Q) ▷ DSA-ordered schedule π ∗ 7: H += ca cb cc cd · V 8: end for 9: Gcart ← HRR bra (H, Oab ) ▷ stepwise spherical fusion Gsph ← (T(la ) ⊗ T(lb ) ) · Gcart Gcart ← HRRket (Gsph , Ocd ) Gsph ← (T(lc ) ⊗ T(ld ) ) · Gcart 10: J, K += S CATTER(Gsph , D)

overhead-free path outperforms the two-kernel Tier 2 pipeline despite the minor spill. P 2) Tier 2: Buffered Two-Kernel Pipeline ( li ≥ 7): When Nacc alone exceeds the register budget (Regime II), accumulator variables cannot reside in registers regardless of VRR scheduling. This structural constraint—not an implementation choice—necessitates a fundamentally different architecture with two cooperating kernels (Figure 2 (a)).

P Regime I ( li ≤ 6). The combined footprint of transient VRR intermediates and persistent accumulators still fits in registers. Nacc reaches at most 225 (at l = (2, 2, 1, 1): 5 × 5 × 3 × 3); after DSA, PeakLiveVRR ≲ 50. Shared register occupancy therefore remains within the hardware budget with at most negligible spill (<1 KB), so a register-only fused kernel is emitted (§III-C1). P Regime II ( li ≥ 7). The persistent accumulator bank already dominates the shared footprint. At (2, 2, 2, 1): 5 × 5 × 5 × 3 = 375 doubles = 750 registers—3× the hardware limit before any transient VRR term is added. In this regime, keeping accumulators on chip is structurally impossible, so a buffered two-kernel pipeline is emitted (§III-C2) while VRR intermediates are still kept register-resident as much as possible. The tier boundaries are not tuned hyperparameters: they follow directly from Eq. (6) and the 255-register limit. P 1) Tier 1: Register-Only Kernels ( li ≤ 6): When the shared footprint fits in registers (Regime I), FusionRCG generates a single fused kernel that keeps all intermediates in the register file. Phase 1 (VRR) executes inside four nested primitive loops using DSA-scheduled register variables. Results accumulate into Phase 2 (HRR) input variables via +=, exploiting the well-known accumulate-then-transform restructuring [1] that avoids n4p redundant HRR evaluations. Phase 2 executes once after the loops, and the output contraction (J/K matrix update via atomicAdd) is fused inline—no intermediate writes to global memory at any point. Algorithm 2 illustrates the generated kernel structure. This path has zero fixed overhead: no global buffer allocation, no transpose, no extra kernel launch. For quartets where residual spill is below ∼1 KB (e.g., 484 B for (1, 1, 1, 1)), this

Kernel 1 (Evaluation). A warp-sized batch of 32 threads cooperatively evaluates their respective shell quartets. VRR intermediates remain in registers (bounded by DSA); HRR accumulators are written to an interleaved global buffer layout, enabling coalesced memory access across the warp. A two-pass lifetime allocator (forward dry-run for birth/death timestamps, reverse-time best-fit packing) enables slot reuse: for l = (2, 2, 2, 2), this compresses 3,545 slots to 1,865—a 1.9× reduction. Kernel 2 (Contraction). Each thread block processes one task with 64 threads. Both the density matrix tile D and the intermediate buffer are loaded into shared memory (≤ 48 KB), delivering ∼19 TB/s bandwidth—10× over HBM. The conP traction Jcd += ab Dab · I[a, b, c, d] exploits 8-fold ERI permutation symmetry via scale factors. IV. I MPLEMENTATION FusionRCG is implemented as an offline code generator rather than a hand-written collection of CUDA kernels. This design is necessary because each angular-momentum quartet induces a different recurrence graph, liveness profile, tier decision, and lowering strategy; encoding these combinations manually would be both error-prone and difficult to maintain as the optimization rules evolve. For each parameter tuple l = (l1 , l2 , l3 , l4 ), the generator therefore executes a nine-stage

7

Throughput and Speedup from DSA Scheduling

E7KURXJKSXWVSHHGXSIURPD[LVSULRULW\

 ×

 

×

×

×

× × ×

         _ _ _ _ _ _ _ _ _         

)XVLRQ5&*

$[LVGLVDEOHG

0

0

× × × × ×

×

×

FusionRCG

[

×

[ [

         _ _ _ _ _ _ _ _ _         

)XVLRQ5&* $[LVGLVDEOHG

DSA disabled

Speedup

[

6SHHGXS

Fig. 5: Ablation P of axis-priority graph construction across six quartets ( li = 2–10). (a) PeakLive reduction factors. (b) Throughput gain from axis priority.

3.0

1000.0M 100.0M

2.2x

2.1x

1.9x

2.5x

2.5

Speedup (x-fold)

×

[

Shell quartets / second



×

6SHHGXS

×

6KHOOTXDUWHWVVHFRQG

3HDN/LYHLQWHUPHGLDWHV

D3HDN/LYHUHGXFWLRQIURPD[LVSULRULW\ 

2.0 1.5

1.0x

1.0x

10.0M

1.0x

1.0x

1.0 0.5

)

|1,0

(1,0

)

|2,0

(2,0

)

|2,0

(3,0

)

|3,0

(3,0

)

|2,2

(3,2

)

|3,2

(3,2

)

|3,3

(4,1

)

|4,2

(4,2

0.0

P Fig. 6: Ablation of DSA scheduling across six quartets ( li = 4–10). Throughput and speedup from DSA.

pipeline: (1) select tier from Eq. (6); (2) construct Phase 1 computation graph using the axis-priority rule; (3) schedule with DSA and insert barriers for Tier 2; (4) [Tier 2] dry-run lifetime analysis and buffer slot packing; (5) build Phase 2 computation graph with stepwise Cartesian-to-spherical fusion; (6) schedule Phase 2; (7) generate the accumulate-thentransform loop structure; (8) emit specialized CUDA kernel (.cu); (9) compile to device object with nvcc.

conventional la →lc , x→y→z ordering), No DSA (a safe but non-optimized topological order), and No spherical (spherical fusion disabled). We also compare Tier 1-only vs. Tier 2only execution to motivate adaptive routing. Each variant is compiled as a complete library so that compiler reports and hardware counters reflect the true fused code path. 1) Effect of Axis-Priority Graph Construction: Figure 5 isolates the axis-priority rule (Algorithm 1) by comparing the full pipeline against the conventional ordering (la →lc , x→y→z), labeled “Axis disabled” in the figure. Both variants use the same DSA schedule; all differences arise solely from the graph topology choice. The benefit scales consistently with angular momentum P (Figure 5). At low li both topologies fit within the register budget and the throughput difference is marginal. Once register P pressure becomes significant ( li ≥ 7), axis priority compresses PeakLive by up to 4.2× and improves throughput by up to +1.47%, confirming that the gain concentrates precisely where spill cost is highest. 2) Effect of DSA Scheduling: Figure 6 isolates DSA by comparing the full pipeline against a variant using a safe but non-optimized topological order on the same graph structure. P For light quartets ( li ≤ 6), any reasonable topological order fits within the register file and DSA adds little (Figure 6). P Once spilling sets in ( li ≥ 8), DSA delivers up to 2.3× throughput gain and reduces dynamic local-memory traffic by up to 142×, confirming the causal chain: optimized evaluation order → fewer compiler spills → reduced off-chip traffic → higher throughput. 3) Effect of Spherical Fusion: Figure 7 compares the full pipeline (with stepwise spherical fusion) against a variant operating entirely in the Cartesian basis (“No spherical”). At low angular momentum the Cartesian-to-spherical size gap is small and the effect is negligible (Figure 7). As l increases, spherical fusion eliminates a growing fraction of floating-point operations—up to 12.81% FLOPs reduction at the highest tested quartet—translating to a peak throughput speedup of 1.91×. The stepwise strategy is key: each index is

Kernel specialization. For Tier 1, the final contraction is fused inline with zero intermediate global writes. For Tier 2, the dedicated _jk_transform kernel tiles the density matrix D into shared memory and parallelizes the contraction across 64 threads per block, exploiting 8-fold ERI permutation symmetry via compile-time scale factors.
Build system. The parameter space has lmax4 +4 unique angular-momentum configurations (120 for lmax = 4). Including first-derivative variants and both tier architectures, the full library comprises 565 kernel files. Generation is embarrassingly parallel: each l-configuration is independent, so the pipeline is distributed across MPI ranks. On a 32-core node, the complete library (generation + nvcc compilation) builds in several hours. V. E XPERIMENTAL E VALUATION A. Experimental Setup Platform. All experiments run on NVIDIA A100-80GB GPUs (sm 80, CUDA 12.4). Single-GPU evaluations use one A100; multi-GPU scaling tests use up to 64 GPUs across 8 nodes. Baseline. We compare against GPU4PySCF [9] (GPUaccelerated PySCF), the current state-of-the-art GPU quantum chemistry framework for Hartree–Fock and DFT calculations. GPU4PySCF employs Rys-quadrature-based CUDA kernels and serves as the strongest publicly available GPU baseline for end-to-end SCF performance. All reported timing results are averaged over three independent runs. B. Ablation Study To validate each optimization component in isolation, we generate variant kernel libraries from the same code generator by toggling one component at a time: Common-way axis (the

8

Throughput Speedup from Spherical Fusion 2.68x

100.0M 1.28x

1.70x

1.46x

1.59x

2.5x

2.08x

1.95x

1.80x

3.0x

Speedup (x-fold)

Shell quartets / second

1000.0M

2.0x 1.5x

10.0M

1.0x 0.5x

1.0M

(2,0

|2,0

)

(3,0

|2,0

)

)

|3,0

(3,0

FusionRCG

)

|2,2

(2,2

)

|2,2

(3,2

)

|3,2

(3,2

No spherical

)

|3,3

(4,1

)

|4,2

(4,2

0.0x

Fig. 9: Time per SCF iteration: GPU4PySCF vs. FusionRCG on a single A100 GPU. (a) Varying system size with ccpVDZ (lmax = 2): five molecular systems from 60 to 272 atoms. (b) Higher angular-momentum basis sets: cc-pVTZ (lmax = 3) and cc-pVQZ (lmax = 4) on three systems each. The red line (right axis) shows the speedup of FusionRCG over GPU4PySCF.

Speedup

Fig. P 7: Ablation of stepwise spherical fusion across six quartets ( li = 4–10). Throughput speedup over the no-spherical variant.

routes each quartet to the execution regime where it performs best. These components are complementary: their combined effect in the full FusionRCG pipeline exceeds any single optimization in isolation. C. End-to-End SCF Performance We compare the time per SCF iteration against GPU4PySCF on a single A100. Since GPU4PySCF tightly integrates ERI evaluation within its SCF solver, isolated kernel extraction is infeasible; we therefore compare at the application level, where both systems execute the identical SCF calculation on each molecular system.

P

Fig. 8: Tier 1-only vs. Tier 2-only throughput over li = 0– P 16. At low l Tier 1 wins due to zero pipeline overhead; i P at high li Tier 2 winsP by eliminating register spills. Tier 1 fails to compile beyond li = 12.

Varying system size at fixed basis (cc-pVDZ). Figure 9(a) evaluates five molecular systems of increasing size—from the water cluster (H2 O)20 to the mini-protein Trp-cage (272 atoms)—with the cc-pVDZ basis set (lmax = 2). FusionRCG consistently outperforms GPU4PySCF, achieving 1.33–1.62× speedup across the entire range. The advantage is largest on small systems and gradually narrows as system size grows, because cc-pVDZ involves primarily low angular-momentum shells and non-ERI SCF components (e.g., diagonalization, density fitting) constitute a growing fraction of each iteration.

contracted and projected immediately rather than materializing the full Cartesian tensor, keeping the peak intermediate at the current-stage size and eliminating redundant arithmetic at the algebraic level. 4) Effect of Adaptive Multi-Tier Routing: Figure 8 compares Tier 1-only (register-only for all quartets) against Tier 2P only (buffered pipeline for all quartets) across li = 0–16. P At low angular momentum ( li ≤ 6), Tier 1 outperforms Tier 2 because the entire working set fits in registers, making Tier 2’s additional pipeline stages (buffer allocation, layout P conversion, dual kernel launch) unnecessary overhead. As li grows, Tier 1’s register working set overflows the 255-register limit, causing massive spilling; Tier 2 avoids this because its buffered two-kernel design keeps the per-thread memory footprint far lower, absorbing the overflow through shared and P global memory rather than register spills. Beyond li = 12, Tier 1 fails to compile entirely. This crossover motivates P FusionRCG’s adaptive routing: Tier 1 for li ≤ 6 and Tier 2 P li ≥ 7, a boundary that follows directly from Eq. (6). for Taken together, the four ablations confirm that axis-priority graph construction and DSA scheduling jointly control register pressure and spill traffic, spherical fusion reduces both arithmetic and memory footprint at high l, and adaptive tiering

Higher angular-momentum basis sets. Figure 9(b) demonstrates FusionRCG’s increasing advantage as the basis set incorporates higher angular-momentum functions. With ccpVTZ (lmax = 3), speedups range from 1.65× to 2.33×. Moving to cc-pVQZ (lmax = 4), where f - and g-type shell quartets become prevalent, FusionRCG achieves its peak advantage of 3.09× on the drug cluster Aspirin21 , with all tested systems exceeding 1.98×. This trend directly reflects the scaling of FusionRCG’s optimizations: liveness-aware orchestration and spherical fusion provide progressively larger compression as angular momentum increases, while adaptive tiering becomes structurally essential for the g-function quartets that dominate cc-pVQZ workloads.

9

Fig. 10: Average SCF speedup of FusionRCG over GPU4PySCF across basis sets with progressively higher angular momentum. The advantage scales from 1.5× (cc-pVDZ, lmax = 2) to 2.0× (cc-pVTZ, lmax = 3) to 2.4× (cc-pVQZ, lmax = 4).

Scaling trend across basis hierarchy. Figure 10 summarizes the aggregate trend: the average speedup over GPU4PySCF increases monotonically from 1.5× at cc-pVDZ to 2.0× at cc-pVTZ and 2.4× at cc-pVQZ. This consistent scaling validates our central thesis—that co-optimizing graph structure, intermediate representation, and memory mapping yields compounding benefits as the computational complexity of the recurrence grows. The 2.4× average improvement at cc-pVQZ, where high-l integrals dominate the SCF cost, demonstrates that FusionRCG effectively rescues these workloads from the memory-bound regime. To demonstrate that kernel-level register optimizations translate to system-level HPC scalability, we evaluate FusionRCG’s strong scaling on a multi-node A100 cluster for Ubiquitin (PDB: 1UBQ [21], 1,231 atoms, 76 residues) with the cc-pVTZ basis set. As shown in Figure 11, FusionRCG scales to 5.97× speedup at 64 GPUs (8 nodes) relative to the 8-GPU baseline, maintaining 75% parallel efficiency despite inter-node communication overhead and load imbalance among shell-quartet batches. The sub-linear scaling is primarily attributable to the final J/K-matrix reduction across nodes and to workload heterogeneity among shell-quartet classes, which causes faster GPUs to idle while waiting for the slowest batch to complete. This demonstrates that per-kernel register optimizations compose well with distributed-memory parallelism.

Fig. 11: Multi-GPU strong scaling on Ubiquitin (PDB: 1UBQ, 1,231 atoms)/cc-pVTZ. (a) Speedup relative to 8 GPUs; the dashed line indicates ideal linear scaling. (b) Parallel efficiency, defined as speedup/(GPU count/8).

performance for complex basis sets. LibintX [22] ports the widely used recursive ERI evaluation to GPUs, representing the most mature publicly available GPU-native ERI library for molecular integrals. Mako [23] reformulates ERIs as matrix multiplications to leverage Tensor Cores via CUTLASS, demonstrating impressive speedups. Though its closed-source nature precludes direct kernel-level comparison, it represents an orthogonal algorithmic approach to our thread-local recurrence method. In contrast, FusionRCG addresses the fundamental register-limit bottleneck of the recursive per-thread model by dynamically co-optimizing the algebraic recurrence paths and instruction scheduling for arbitrary l, enabling automated, high-performance kernel generation without manual tuning. A. Multi-GPU Scalability Graph scheduling and register pressure management. General-purpose minimum-register instruction scheduling is inherently NP-hard [24]. Classical compiler heuristics [16]– [18], [20] and GPU-specific optimizations, such as occupancyaware instruction reordering [6], [19] or kernel fission [25], operate on a fixed dependence graph or scheduling region. Similar assumptions appear in dependency-rich GPU kernels such as sparse triangular solve, where an analysis phase first builds a graph from the matrix sparsity pattern and then executes that graph level by level [26]. These methods improve execution order, synchronization, or resource usage on a prescribed DAG, but they do not change the graph itself. HGP differs in two essential respects. First, the same target ERI admits multiple legal recurrence paths, so the computation

VI. R ELATED W ORK GPU-accelerated quantum chemistry. TeraChem [2] pioneered GPU ERI evaluation. While early implementations relied heavily on hand-tuned kernels for low angular momenta, subsequent developments have incorporated metaprogramming to support higher l; nevertheless, managing the exponential growth of register pressure remains a fundamental bottleneck. GPU4PySCF [9] provides an efficient CUDAaccelerated backend for PySCF and serves as our primary end-to-end baseline (§V-C). Although it generates kernels for higher angular momenta, it relies heavily on local memory and experiences significant register spills, leading to degraded

10

VII. C ONCLUSION

graph is not a fixed input to scheduling: different axis choices induce genuinely different graph topologies and therefore different live ranges. Second, the dominant resource bottleneck is not merely dependency satisfaction on a fixed graph, but the overlap between a transient VRR live frontier and a persistent HRR accumulator bank under a hard per-thread register budget. General schedulers lack the domain information needed to exploit this algebraic degree of freedom. FusionRCG therefore couples graph construction with register-aware scheduling, reducing pressure before the generated kernel reaches NVCC rather than only reordering a predetermined DAG.

We have presented FusionRCG, a framework that addresses the GPU memory wall for hierarchical recursive computation graphs through three co-designed optimizations. Our key insight—that the recurrence graph’s structure is a degree of freedom, not a fixed constraint—enables liveness-aware graph orchestration that reduces peak simultaneously live intermediates by up to 3.2× and boosts throughput by up to 44.7%. Combined with stepwise spherical fusion that cuts FLOPs by up to 12.81% and delivers up to 1.91× kernel speedup, and an adaptive two-tier kernel architecture derived from registerpressure decomposition, FusionRCG achieves zero spill for all P angular-momentum classes up to li = 6 and extends feasiP ble GPU execution to the full g-function range ( li = 16). In end-to-end SCF comparison with GPU4PySCF on NVIDIA A100, FusionRCG achieves up to 3.09× speedup, with the average advantage scaling from 1.5× at cc-pVDZ to 2.4× at cc-pVQZ—confirming that the benefit grows with angularmomentum complexity. Furthermore, FusionRCG scales to 64 GPUs with 75% parallel efficiency, demonstrating that per-kernel register optimizations compose effectively with distributed-memory parallelism. Future directions include extending the framework to derivative integrals required for geometry optimization and molecular dynamics, and exploring Tensor Core co-execution for hybrid recurrence–matmul workloads.

Recursive and irregular GPU algorithms. A separate line of work studies recursive or highly irregular GPU workloads such as tree traversals, divide-and-conquer algorithms, and nested parallelism. For recursive tree traversals, prior work focuses on mitigating control and memory divergence across traversals by regrouping threads or by using hybrid CPU–GPU scheduling [27]. For divide-and-conquer algorithms, GPU implementations such as QuickHull reformulate recursively generated subproblems as segmented data-parallel operations over contiguous array ranges [28]. More generally, nestedparallel and dynamic-parallel approaches use child-kernel launches and workload consolidation to cope with runtimegenerated work, while identifying kernel-launch overhead and GPU underutilization as the main difficulty [29], [30]. These techniques do not transfer directly to HGP. For a given shell quartet, the full HGP recurrence can be constructed statically at code-generation time; the core challenge is therefore not runtime task discovery, warp regrouping, or recursive kernel-launch granularity. Instead, HGP suffers from a compile-time, per-thread live-state explosion inside a small but dense algebraic recurrence DAG. Methods designed for segmented subproblem management, dynamic load balancing, or nested kernel consolidation do not reduce the VRR live frontier, the Phase 2 accumulator bank, or their overlap in the register file. FusionRCG targets this gap by treating HGP as a domain-specific recursive computation graph whose topology, schedule, representation, and memory placement must be cooptimized jointly.

ACKNOWLEDGMENT This project was supported by the Zhongguancun Academy under the Internal Research Grant No. C20250501. R EFERENCES [1] M. Head-Gordon and J. A. Pople, “A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations,” The Journal of Chemical Physics, vol. 89, no. 9, pp. 5777–5786, 1988. [Online]. Available: https://doi.org/10.1063/1.455553 [2] I. S. Ufimtsev and T. J. Martı́nez, “Quantum chemistry on graphical processing units. 1. strategies for two-electron integral evaluation,” Journal of Chemical Theory and Computation, vol. 4, no. 2, pp. 222–231, 2008. [Online]. Available: https://doi.org/10.1021/ct700268q [3] H. Laqua, J. Kussmann, and C. Ochsenfeld, “Accelerating seminumerical fock-exchange calculations using mixed single- and double-precision arithmetic,” The Journal of Chemical Physics, vol. 154, no. 21, p. 214116, 2021. [Online]. Available: https://doi.org/10.1063/5.0045084 [4] E. F. Valeev, “Libint: Machine-generated library for efficient evaluation of molecular integrals over Gaussians,” 2025. [Online]. Available: https://evaleev.github.io/libint/ [5] F. Neese, “The SHARK integral generation and digestion system,” Journal of Computational Chemistry, vol. 44, no. 3, pp. 381–396, 2023. [Online]. Available: https://doi.org/10.1002/jcc.26942 [6] V. Volkov, “Better performance at lower occupancy,” 2010, GPU Technology Conference. [Online]. Available: https://dmacssite.github. io/materials/volkov10-GTC.pdf [7] A. W. Appel and L. George, “Optimal spilling for CISC machines with few registers,” in Proceedings of the ACM SIGPLAN 2001 Conference on Programming Language Design and Implementation, 2001, pp. 243–253. [Online]. Available: https://doi.org/10.1145/378795.378854

Generation of Domain-specific code. Frameworks like ATLAS/FFTW [31], [32] (for linear algebra and FFTs), Halide [33], and TVM [34] (for deep learning) have demonstrated the immense power of auto-tuning and domain-specific compilation. While these tools excel at optimizing loop nests, memory tiling, and tensor operations, they are not designed to handle the combinatorial explosion of intermediate variables inherent in deeply recursive quantum chemistry integrals. FusionRCG uniquely expands the scope of domain-specific compilation by tightly coupling mathematical recurrence generation, algorithmic path finding, and hardware-aware register management into a unified ERI code generator.

11

[25] A. Sethia and S. A. Mahlke, “Equalizer: Dynamic tuning of GPU resources for efficient execution,” in Proceedings of the 47th Annual IEEE/ACM International Symposium on Microarchitecture, 2014, pp. 647–658. [Online]. Available: https://doi.org/10.1109/MICRO.2014.16 [26] M. Naumov, “Parallel solution of sparse triangular linear systems in the preconditioned iterative methods on the GPU,” NVIDIA Corporation, Tech. Rep. NVR-2011-001, 2011. [Online]. Available: https://research.nvidia.com/publication/2011-06 parallel-solution-sparse-triangular-linear-systems-preconditioned-iterative [27] W. Liu, B. Schmidt, and W. W. Hwu, “Hybrid CPU-GPU scheduling and execution of tree traversals,” in Proceedings of the 21st ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, 2016, pp. 330–340. [Online]. Available: https://doi.org/ 10.1145/2851141.2851174 [28] G. Mei, N. Xu, and L. Xu, “A sample implementation for parallelizing divide-and-conquer algorithms on the GPU,” Heliyon, vol. 4, no. 3, p. e00512, 2018. [Online]. Available: https://doi.org/10.1016/j.heliyon. 2018.e00512 [29] W. Li, G. Shobaki, and T. El-Ghazawi, “Nested parallelism on GPU: Exploring parallelization templates for irregular loops and recursive computations,” in 2015 44th International Conference on Parallel Processing, 2015, pp. 595–604. [Online]. Available: https://doi.org/10.1109/ICPP.2015.107 [30] G. Shobaki, W. Li, and T. El-Ghazawi, “Compiler-assisted workload consolidation for efficient dynamic parallelism on GPU,” in 2016 IEEE International Parallel and Distributed Processing Symposium, 2016, pp. 534–543. [Online]. Available: https://doi.org/10.1109/IPDPS.2016.110 [31] R. C. Whaley and J. J. Dongarra, “Automatically tuned linear algebra software,” in Proceedings of the ACM/IEEE Conference on Supercomputing, 1998. [Online]. Available: https://doi.org/10.1109/SC. 1998.10004 [32] M. Frigo and S. G. Johnson, “FFTW: An adaptive software architecture for the FFT,” in Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, vol. 3, 1998, pp. 1381– 1384. [Online]. Available: https://doi.org/10.1109/ICASSP.1998.681704 [33] J. Ragan-Kelley, C. Barnes, A. Adams, S. Paris, F. Durand, and S. Amarasinghe, “Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines,” in Proceedings of the 34th ACM SIGPLAN Conference on Programming Language Design and Implementation, 2013, pp. 519–530. [Online]. Available: https://doi.org/10.1145/2499370.2462176 [34] T. Chen, T. Moreau, Z. Jiang et al., “TVM: An automated endto-end optimizing compiler for deep learning,” in 13th USENIX Symposium on Operating Systems Design and Implementation, 2018, pp. 578–594. [Online]. Available: https://www.usenix.org/conference/ osdi18/presentation/chen

[8] S. Williams, A. Waterman, and D. Patterson, “Roofline: An insightful visual performance model for multicore architectures,” Communications of the ACM, vol. 52, no. 4, pp. 65–76, 2009. [Online]. Available: https://doi.org/10.1145/1498765.1498785 [9] Q. Sun, T. Zhu, N. S. Blunt et al., “Introducing GPU acceleration into the Python-based simulations of chemistry framework,” The Journal of Physical Chemistry A, vol. 129, no. 5, pp. 1459–1468, 2025. [Online]. Available: https://doi.org/10.1021/acs.jpca.4c05876 [10] A. Asadchev, V. Allada, J. Felder, B. M. Bode, M. S. Gordon, and T. L. Windus, “Uncontracted Rys quadrature implementation of up to G functions on graphical processing units,” Journal of Chemical Theory and Computation, vol. 6, no. 3, pp. 696–704, 2010. [Online]. Available: https://doi.org/10.1021/ct9005079 [11] E. Palethorpe and G. M. J. Barca, “High-performance, high-angularmomentum J engine on graphics processing units,” Journal of Chemical Theory and Computation, vol. 21, no. 19, pp. 9388–9403, 2025. [Online]. Available: https://doi.org/10.1021/acs.jctc.5c00775 [12] G. Samu and M. Kállay, “Efficient evaluation of three-center coulomb integrals,” The Journal of Chemical Physics, vol. 146, no. 20, p. 204101, 2017. [Online]. Available: https://doi.org/10.1063/1.4983393 [13] A. Asadchev and E. F. Valeev, “Implementation of McMurchie– Davidson algorithm for Gaussian AO integrals suited for SIMD processors,” The Journal of Physical Chemistry A, vol. 129, no. 42, pp. 9788–9797, 2025. [Online]. Available: https://doi.org/10.1021/acs.jpca. 5c04136 [14] H. B. Schlegel and M. J. Frisch, “Transformation between cartesian and pure spherical harmonic Gaussians,” International Journal of Quantum Chemistry, vol. 54, no. 2, pp. 83–87, 1995. [Online]. Available: https://doi.org/10.1002/qua.560540202 [15] R. Sethi and J. D. Ullman, “The generation of optimal code for arithmetic expressions,” Journal of the ACM, vol. 17, no. 4, pp. 715–728, 1970. [Online]. Available: https://doi.org/10.1145/321607.321620 [16] R. Castañeda Lozano and C. Schulte, “Survey on combinatorial register allocation and instruction scheduling,” ACM Computing Surveys, vol. 52, no. 3, pp. 62:1–62:50, 2019. [Online]. Available: https://doi.org/10.1145/3340313 [17] S. Park, A. Nicolau, A. Shrivastava, Y. Paek, N. Dutt, and E. Earlie, “Bypass aware instruction scheduling for register file power reduction,” ACM SIGPLAN Notices, vol. 41, no. 7, pp. 173–181, 2006. [Online]. Available: https://doi.org/10.1145/1159974.1134675 [18] C. W. Kessler, “Scheduling expression DAGs for minimal register need,” Computer Languages, vol. 24, no. 1, pp. 33–53, 1998. [Online]. Available: https://doi.org/10.1016/S0096-0551(98)00002-2 [19] G. Shobaki, A. Kerbow, and S. Mekhanoshin, “Optimizing occupancy and ILP on the GPU using a combinatorial approach,” in Proceedings of the 18th ACM/IEEE International Symposium on Code Generation and Optimization, 2020, pp. 133–144. [Online]. Available: https: //doi.org/10.1145/3368826.3377918 [20] G. Shobaki, J. Bassett, M. Heffernan, and A. Kerbow, “Graph transformations for register-pressure-aware instruction scheduling,” in Proceedings of the 31st ACM SIGPLAN International Conference on Compiler Construction, 2022, pp. 41–53. [Online]. Available: https://doi.org/10.1145/3497776.3517771 [21] S. Vijay-Kumar, C. E. Bugg, and W. J. Cook, “Structure of ubiquitin refined at 1.8 Å resolution,” Journal of Molecular Biology, vol. 194, no. 3, pp. 531–544, 1987. [Online]. Available: https://doi.org/10.1016/0022-2836(87)90679-6 [22] A. Asadchev and E. F. Valeev, “LibintX: High-performance library for scalable molecular integral evaluation,” 2023. [Online]. Available: https://github.com/ValeevGroup/LibintX [23] H. Han, K. Li, F. Ju, Q. Li, H. An, Y. Chen, Y. Zhang, T. Cao, and M. Yang, “Matrix is all you need: Rearchitecting quantum chemistry to scale on AI accelerators,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2025, pp. 2126–2142. [Online]. Available: https://doi.org/10.1145/3712285.3759829 [24] R. Sethi, “Complete register allocation problems,” SIAM Journal on Computing, vol. 4, no. 3, pp. 226–248, 1975. [Online]. Available: https://doi.org/10.1137/0204020

12