arXiv:2609.06449v1 [cs.DS] 6 Sep 2026
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution Shina Guo
Weiguo Gao
Yuan Tang∗
College of Computer Science and Artificial Intelligence, Fudan University Shanghai, P. R. China [email protected]
School of Mathematical Sciences, School of Data Science, Fudan University Shanghai, P. R. China [email protected]
College of Computer Science and Artificial Intelligence, Fudan University Shanghai, P. R. China [email protected]
Abstract
1
Fast Fourier transforms (FFTs) on parallel machines incur two distinct forms of data movement: processor-local transfers through the memory hierarchy and global redistribution between processors. Existing four-step FFT organizations can switch the active transform dimension with one global transpose-like exchange, but this property alone does not yield cache-efficient local computation. Conversely, cache-oblivious FFTs obtain asymptotically optimal local memory traffic through recursively applied layout transformations, whose direct materialization can require additional datarearrangement passes and global exchanges. We present PACO, a fully cache-oblivious parallel FFT framework that reconciles these objectives. PACO executes
The fast Fourier transform (FFT) is a fundamental primitive in scientific computing, signal processing, and data analysis. The Cooley– Tukey algorithm reduces the arithmetic cost of an 𝑁 -point discrete Fourier transform (DFT) from Θ(𝑁 2 ) to Θ(𝑁 log 𝑁 ) operations [4]. On modern parallel machines, however, arithmetic is only one component of the execution cost. An FFT must move data through each processor’s local memory hierarchy and redistribute data between processors when the active transform dimension changes. Thus, an efficient parallel FFT should simultaneously control processor-local cache traffic and interprocessor communication. Existing FFT organizations typically optimize one of these forms of data movement more directly than the other. Four-step and sixstep FFT organizations reshape a one-dimensional transform into matrix views and use transpose-like operations to switch the active transform dimension [2, 13]. With compatible local kernels and data distributions, a two-factor four-step organization can use one global slab redistribution while producing the canonical DFT coefficient order in a factor-swapped output view. This algebraic one-exchange structure, however, does not by itself guarantee fully cache-oblivious local execution. Conversely, cache-oblivious FFTs obtain asymptotically optimal ideal-cache complexity through recursively applied transpose-like layout transformations [8]. Directly materializing these transformations in a distributed setting may introduce extra local data passes and, when an affected digit field is distributed, additional global redistributions. Other lowcommunication FFTs reduce the number of global all-to-all operations through alternative factorizations, approximation, oversampling, or additional neighbor-level communication [7, 18]; their numerical and distributional contracts differ from the exact FFT setting considered here. Problem setting: We consider an exact bounded-radix Cooley– Tukey FFT under an explicit source and target ownership contract. Let
LocalFFT −→ OneGlobalPermutation −→ LocalFFT. Its local stages recursively partition both transform and batch dimensions without knowledge of the cache parameters. Rather than materializing the transpose-like layouts induced by this recursion, PACO defers them. We show that the deferred layouts compose into a base-𝑏 digit-reversal permutation, which PACO fuses with the redistribution already required to change the local transform dimension. The resulting middle stage is a perfectly balanced parallel cache-oblivious digit-reversal permutation: each source– destination processor pair exchanges exactly 𝑁 /𝑝 2 elements (the 𝑝 self-pairs incur no network traffic). Under an exact base-𝑏 slab decomposition in a hybrid idealcache/BSP model, PACO computes an 𝑁 -point DFT exactly with maximum per-processor work 𝑁 Θ log 𝑁 𝑝 and maximum per-processor cache complexity 𝑁 Θ 1 + log𝑀 𝑁 . 𝑝𝐵 It uses exactly one global redistribution round. PACO returns canonical logical DFT coefficients under a factor-swapped target slab ownership. The single redistribution is necessary under the stated no-replication ownership model, while its communication volume is optimal for the prescribed fused permutation and source–target slab distributions. ∗ Corresponding Author; Also affiliated with Shanghai Key Lab. of Intelligent Informa-
tion Processing.
Introduction
𝑁 = 𝑁1 𝑁2 . The input is viewed as an 𝑁 1 × 𝑁 2 row-major matrix distributed by 𝑁 2 -column slabs. The output is viewed as an 𝑁 2 × 𝑁 1 rowmajor matrix distributed by 𝑁 1 -column slabs. Its coordinate (𝑘 2, 𝑘 1 ) represents 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ]. Hence the output has canonical logical DFT coefficient indices, but its physical ownership is factor-swapped relative to the input distribution. Restoring the initial slab ownership is not part of the problem and may require a further redistribution.
Shina Guo, Weiguo Gao, and Yuan Tang
This paper presents PACO, a Processor-Aware Cache-Oblivious FFT framework that combines fully cache-oblivious local computation with exactly one global redistribution. PACO follows the fixed structure
Problem and ownership contract. Fix a constant radix 𝑏 = 2𝑡 = 𝑂 (1). We assume
Its local stages use a cache-oblivious batched FFT scheduler that recursively partitions both transform and batch dimensions without knowledge of the cache capacity 𝑀 or block size 𝐵. The central observation is that a transpose-free cache-oblivious local FFT does not leave its result in arbitrary disorder. Instead, the layout transformations deferred during local recursion compose into a structured base-𝑏 digit-reversal permutation. PACO materializes this complete permutation exactly once, at the redistribution already required between the two top-level Cooley–Tukey transform groups. Thus, the middle stage simultaneously corrects the first local stage’s deferred layout, swaps the two factor fields, prepares the layout consumed by the second local stage, and transfers each item to its factor-swapped target owner. PACO does not introduce a separate global communication phase merely to repair digit order. Contributions. This paper makes the following contributions: • We design a fully cache-oblivious batched Cooley–Tukey FFT scheduler that recursively partitions transform and batch dimensions, avoids recursion-level physical transposes, and achieves asymptotically optimal local work and cache complexity under the stated memory model. • We show that the lazy layouts induced by the local recursion compose into PARA-CO-DRP𝑏 , a base-𝑏 PARAllel Cache-Oblivious Digit-Reversal Permutation, and we fuse that permutation with the ownership change required between the two top-level transform groups. The resulting PARA-CO-DRP𝑏 middle stage uses exactly one global redistribution and routes exactly 𝑁 /𝑝 2 items from every source processor to every destination processor. • Under an exact base-𝑏 slab decomposition in a hybrid idealcache/BSP model, we prove maximum per-processor work 𝑁 Θ log 𝑁 𝑝 and maximum per-processor cache complexity 𝑁 Θ 1 + log𝑀 𝑁 . 𝑝𝐵 We further prove that one global redistribution round is necessary under the no-replication ownership model, and that PACO’s communication volume is optimal for the prescribed fused slab-to-slab permutation.
2
Model and Main Result
We study the canonical length-𝑁 discrete Fourier transform 2𝜋𝚤 𝑦 = Ω𝑁 𝑥, (Ω𝑁 )𝑘,𝑗 = 𝜔 𝑁−𝑘 𝑗 , 𝜔 𝑁 = exp . 𝑁 Let [𝑛] := {0, 1, . . . , 𝑛 − 1}. One complex scalar is one item in both the cache and communication models.
𝑁 1 = 𝑏𝑚 1 ,
𝑁 = 𝑁 1 𝑁 2,
LocalFFT −→ OneGlobalPermutation −→ LocalFFT.
𝑁 2 = 𝑏𝑚 2 ,
and write 𝑚 := 𝑚 1 + 𝑚 2 = log𝑏 𝑁 . There are 𝑝 ≥ 2 processors, indexed by [𝑝]. For the exact balance results, 𝑝, 𝑁 1 , and 𝑁 2 are powers of 𝑏, and 𝑝 | 𝑁 1,
𝑝 | 𝑁2 .
The input vector is viewed as an 𝑁 1 × 𝑁 2 row-major matrix 𝑋 [ 𝑗1, 𝑗2 ] = 𝑥 [ 𝑗1 𝑁 2 + 𝑗2 ],
𝑗1 ∈ [𝑁 1 ],
𝑗2 ∈ [𝑁 2 ],
distributed by equal column slabs: owner0 ( 𝑗1, 𝑗2 ) =
𝑗2 𝑝 . 𝑁2
(1)
Thus, processor 𝑞 initially owns the columns 𝑁2 𝑁2 ≤ 𝑗2 < (𝑞 + 1) . 𝑞 𝑝 𝑝 We factor the canonical frequency index as 𝑘 = 𝑘 1 + 𝑘 2 𝑁 1,
𝑘 1 ∈ [𝑁 1 ],
𝑘 2 ∈ [𝑁 2 ],
and view the output as an 𝑁 2 × 𝑁 1 row-major matrix 𝑌 [𝑘 2, 𝑘 1 ] = 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ]. The target ownership is the column-slab distribution 𝑘1𝑝 owner1 (𝑘 2, 𝑘 1 ) = . 𝑁1
(2)
(3)
Hence, PACO returns canonical DFT coefficients 𝑦 [𝑘], but under a factor-swapped target ownership. It does not restore the initial slab distribution. Machine and execution model. Each processor has a private fully associative ideal cache of 𝑀 items and block size 𝐵, with optimal replacement and initially empty cache. We assume the tall-cache condition 𝑀 = Ω(𝐵 2 ) (4) and the slab-width conditions 𝑁1 𝑁2 ≥ 𝐵, ≥ 𝐵. (5) 𝑝 𝑝 Interprocessor communication follows the BSP model. In one BSP ℎ-relation, every processor sends and receives at most ℎ remote items; its cost is 𝑔ℎ + 𝐿, where 𝑔 is the per-item gap and 𝐿 is the synchronization cost. A global redistribution round consists of local packing, one BSP ℎ-relation, and local unpacking. Packing and unpacking are charged to local work and cache misses, not to ℎ. We use an owner-computes model: input values are initially unreplicated, no processor initially owns the entire input, and a processor obtains a remotely owned value only through an explicit BSP communication superstep. Arithmetic vertices of the fixedradix Cooley–Tukey DAG are not recomputed. Let𝑊max and 𝑄 max denote, respectively, the maximum per-processor local work and cache misses. Let 𝑅 be the number of global redistribution rounds, ℎ the maximum per-processor remote volume in a round, and 𝜇 the aggregate number of remotely transmitted items.
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
Theorem 1 (Main result). Under assumptions (1)–(5), PACO computes the canonical length-𝑁 DFT under the input ownership (1) and factor-swapped target ownership (3). In particular, 𝑌 [𝑘 2, 𝑘 1 ] = (Ω𝑁 𝑥) [𝑘 1 + 𝑘 2 𝑁 1 ] for every 𝑘 1 ∈ [𝑁 1 ] and 𝑘 2 ∈ [𝑁 2 ]. Its maximum per-processor work and cache complexity are 𝑁 𝑊max = Θ log 𝑁 𝑝
Thus, Phase I performs one length-𝑁 1 FFT for every 𝑗2 , Phase II applies the corresponding twiddle factor and redistributes the intermediate values, and Phase III performs one length-𝑁 2 FFT for every 𝑘 1 . Under the input ownership (1), every Phase-I transform is local to one processor. Under the target ownership (3), every PhaseIII transform is local to one processor. The purpose of the middle stage is therefore not to create a new transform factorization, but to change ownership and layout between these two local transform groups.
and 𝑁 1 + log𝑀 𝑁 . 𝑝𝐵 PACO performs exactly one global redistribution round, 𝑄 max = Θ
3.2
For a permutation 𝜋 : [𝑛] −→ [𝑛],
𝑅 = 1,
let 𝑃𝜋 denote the physical-to-logical layout operator (𝑃𝜋 𝑧) [ 𝑗] = 𝑧 [𝜋 ( 𝑗)].
whose BSP volume and aggregate remote migration volume are 𝑁 1 1 ℎ= 1− , 𝜇 =𝑁 1− . 𝑝 𝑝 𝑝 Scope of optimality. The optimality claims in Theorem 1 have different scopes. The one-round bound is optimal for the complete DFT under the stated owner-computes and no-replication assumptions. The work bound is optimal for the fixed-radix Cooley–Tukey computation DAG without arithmetic recomputation. The cache bound matches the lower bound for the local FFT DAGs assigned by the PACO decomposition. Finally, the migration bound is optimal for the prescribed fused digit-reversal permutation between the source and target slab ownerships; it is not a communicationvolume lower bound for all distributed FFT schedules with different intermediate ownerships. The next section develops the two-factor decomposition and the layout algebra that allows PACO to fuse the complete digit reversal with its sole global redistribution.
3
PACO: One Redistribution from Lazy Layouts
PACO combines two transpose-free cache-oblivious local FFT phases with one global redistribution: LocalFFT −→ OneGlobalPermutation −→ LocalFFT. This section explains why the layout changes deferred by the local recursions can be materialized exactly at the ownership change already required between the two top-level Cooley–Tukey transform groups.
3.1
Lazy Layouts of the Local FFT
Two-Factor Cooley–Tukey Decomposition 𝑁 = 𝑁1 𝑁2 .
Writing 𝑘 = 𝑘 1 + 𝑘 2 𝑁 1,
𝑗1, 𝑘 1 ∈ [𝑁 1 ],
𝑗2, 𝑘 2 ∈ [𝑁 2 ],
then 𝜋ℓ (𝑢) = 𝑑 0 ++ · · · ++ 𝑑 ℓ −1 . Digit reversal is an involution. Lemma 2 (Lazy-layout FFT). Consider a transpose-free local radix-𝑏 FFT of length 𝑏 ℓ with natural-order constant-size leaves and twiddle factors evaluated from logical coordinates. Its output-layout form satisfies Aℓ = 𝑃𝜋ℓ Ω𝑏 ℓ . (7) Conversely, the same recursion can consume the digit-reversed input layout and produce a natural-order DFT output: Bℓ = Ω𝑏 ℓ 𝑃𝜋ℓ .
(8)
Proof sketch. A recursive Cooley–Tukey step consists of two child-transform groups, a diagonal twiddle operator, and a final swap of two transform digit fields. PACO performs the child transforms and twiddle multiplication, but does not physically materialize the final digit-field swap. Recursively deferred swaps compose into base-𝑏 digit reversal: for every split ℓ = ℓ1 + ℓ2, the physical digit fields 𝑈 ++ 𝑉 represent the logical fields 𝜋ℓ2 (𝑉 ) ++ 𝜋ℓ1 (𝑈 ). Twiddle factors remain correct because they are evaluated using the logical coordinates represented by each physical item (see the layout-aware twiddle operator, Appendix A.2); equivalently, the diagonal twiddle operator is conjugated by the current layout permutation. Since digit reversal is an involution, conjugating the output-layout form yields the input-layout form. The complete induction and digit-field view invariant are given in the appendix. □
Recall that
𝑗 = 𝑗1 𝑁 2 + 𝑗2,
For 0 ≤ ℓ ≤ 𝑚, let 𝜋ℓ denote base-𝑏 digit reversal on ℓ-digit strings. In particular, if 𝑢 = 𝑑 ℓ −1 ++ · · · ++ 𝑑 0,
where gives 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ] =
𝑁∑︁ 2 −1 𝑁 1 −1 ∑︁ 𝑗2 =0
! 𝑥 [ 𝑗1 𝑁 2 + 𝑗2 ]𝜔 𝑁−𝑘1 1 𝑗1
𝑗1 =0
· 𝜔 𝑁−𝑘1 𝑗2 𝜔 𝑁−𝑘2 2 𝑗2 .
(6)
Lemma 2 shows that a local PACO FFT does not leave its result in arbitrary disorder. It produces a precise and composable layout: base-𝑏 digit reversal. PACO exploits this structure rather than materializing a separate local transpose or digit-reversal pass after every recursive factorization.
Shina Guo, Weiguo Gao, and Yuan Tang
3.3
One Fused Middle Permutation
Phase I applies A𝑚1 to the locally owned length-𝑁 1 transforms indexed by 𝑗2 . After Phase I, let (𝑎, 𝑐) denote a physical array coordinate, where 𝑎 ∈ [𝑁 1 ] is the physical row coordinate and 𝑐 ∈ [𝑁 2 ] is the unchanged second-stage input coordinate. Thus 𝑐 = 𝑗2 , whereas 𝑎 is not the original coordinate 𝑗1 : it represents the logical first-stage frequency 𝑘 1 = 𝜋𝑚1 (𝑎). When the item at physical coordinate (𝑎, 𝑐) is packed for the global redistribution, PACO applies the top-level twiddle factor 𝜔 𝑁−𝑘1𝑐 using its represented logical coordinates 𝑘 1 = 𝜋𝑚1 (𝑎) and 𝑐 = 𝑗2 . PACO then materializes the complete 𝑚-digit reversal, rather than only the 𝑚 1 -digit layout exposed by Phase I. For digit fields 𝑎 ++ 𝑐,
|𝑎| = 𝑚 1,
−𝜋𝑚 (𝑎) 𝑐
𝐷 tw is the diagonal top-level twiddle (𝐷 tw ) (𝑎,𝑐 ) = 𝜔 𝑁 1 , evaluated from the logical coordinate 𝑘 1 = 𝜋𝑚1 (𝑎) carried at physical row 𝑎; 𝑃𝜋𝑚 : C𝑁1 ⊗ C𝑁2 → C𝑁2 ⊗ C𝑁1 is the base-𝑏 digit reversal (9), which simultaneously cancels 𝑃𝜋𝑚1 , swaps the two factors, and installs 𝑃𝜋𝑚2 ; and B𝑚2 ⊗ 𝐼 𝑁1 applies the input-layout length-𝑁 2 transform B𝑚2 = Ω𝑁2 𝑃𝜋𝑚2 to the row factor (batched over the 𝑁 1 columns). Theorem 3 (End-to-end correctness). Under the ownership assumptions of Sect. 2, PACO computes the canonical DFT 𝑦 = Ω𝑁 𝑥 . More precisely, its factor-swapped output view satisfies
|𝑐 | = 𝑚 2,
𝑌 [𝑘 2, 𝑘 1 ] = 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ]
the middle permutation is 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎).
for every (9)
Equation (9) performs four operations simultaneously: (1) it cancels the 𝜋𝑚1 output layout exposed by Phase I; (2) it swaps the 𝑁 1 - and 𝑁 2 -digit fields; (3) it establishes the 𝜋𝑚2 input layout consumed by Phase III; and (4) it changes ownership from owner0 to owner1 , making every complete length-𝑁 2 Phase-III transform processor local. More concretely, an item at physical Phase-I coordinate (𝑎, 𝑐) is routed to the factor-swapped coordinate 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) . Its destination processor is therefore determined by the second coordinate, 𝜋𝑚1 (𝑎)𝑝 . 𝑁1 The detailed packing, routing, local transposition, and unpacking procedure that realizes this mapping cache-obliviously is given in Sect. 4.2. Its exact source–destination balance is analyzed there as well. The ownership change is required by the data dependencies of the second transform group: for a fixed 𝑘 1 , a length-𝑁 2 FFT requires values from all 𝑗2 ∈ [𝑁 2 ], whereas the input slabs distribute those coordinates among the processors. PACO therefore uses the unavoidable communication boundary to materialize the full deferred layout transformation, rather than introducing a separate redistribution solely to repair digit order.
3.4
where A𝑚1 ⊗ 𝐼 𝑁2 applies the output-layout length-𝑁 1 transform A𝑚1 = 𝑃𝜋𝑚1 Ω𝑁1 to the row factor (batched over the 𝑁 2 columns);
End-to-End Correctness
PACO applies the output-layout local operator of Lemma 2 in Phase I, the fused permutation (9) together with the top-level twiddle diagonal in Phase II, and the input-layout local operator in Phase III. Identifying the input array with C𝑁1 ⊗ C𝑁2 (physical coordinate (𝑎, 𝑐), 𝑎 ∈ [𝑁 1 ], 𝑐 ∈ [𝑁 2 ]) and the output array with C𝑁2 ⊗ C𝑁1 (physical coordinate (𝑢, 𝑣), 𝑢 ∈ [𝑁 2 ], 𝑣 ∈ [𝑁 1 ]), PACO’s operator composition is F = B𝑚2 ⊗ 𝐼 𝑁1 𝑃𝜋𝑚 𝐷 tw A𝑚1 ⊗ 𝐼 𝑁2 , (10)
𝑘 1 ∈ [𝑁 1 ],
𝑘 2 ∈ [𝑁 2 ],
and this coefficient is owned by owner1 (𝑘 2, 𝑘 1 ) =
𝑘1𝑝 . 𝑁1
Proof. Write the input array as C𝑁1 ⊗ C𝑁2 and the output array as C𝑁2 ⊗ C𝑁1 , as in (10). By Lemma 2 (in batched form), Phase I realizes A𝑚1 ⊗ 𝐼 𝑁2 = (𝑃𝜋𝑚1 ⊗ 𝐼 𝑁2 )(Ω𝑁1 ⊗ 𝐼 𝑁2 ) and Phase III realizes B𝑚2 ⊗ 𝐼 𝑁1 = (Ω𝑁2 ⊗ 𝐼 𝑁1 )(𝑃𝜋𝑚2 ⊗ 𝐼 𝑁1 ). The middle stage factors as 𝑃𝜋𝑚 = (𝑃𝜋𝑚2 ⊗ 𝐼 𝑁1 ) Πswap (𝑃𝜋𝑚1 ⊗ 𝐼 𝑁2 ), where Πswap : C𝑁1 ⊗ C𝑁2 → C𝑁2 ⊗ C𝑁1 is the factor transpose (Π swap𝑧) [𝑢, 𝑣] = 𝑧 [𝑣, 𝑢]. Indeed, both sides send physical (𝑢, 𝑣) to 𝑧 [𝜋𝑚1 (𝑣), 𝜋𝑚2 (𝑢)], since reversing the 𝑚-digit string 𝑢 ++𝑣 yields 𝜋𝑚1 (𝑣)++𝜋𝑚2 (𝑢). Likewise the layout-aware twiddle diagonal satisfies e (𝑃𝜋𝑚 ⊗ 𝐼 𝑁2 ), 𝐷 tw = (𝑃𝜋𝑚1 ⊗ 𝐼 𝑁2 ) 𝐷 1
e(𝑘1 ,𝑗2 ) = 𝜔 −𝑘1 𝑗2 , 𝐷 𝑁
because conjugating a diagonal by the involution 𝑃𝜋𝑚1 ⊗ 𝐼 𝑁2 merely e(𝜋𝑚 (𝑎),𝑐 ) = 𝜔 −𝜋𝑚1 (𝑎)𝑐 . relabels its entries to 𝐷 𝑁 1 Substituting these identities into (10) and cancelling each adjacent pair (𝑃𝜋𝑚1 ⊗ 𝐼 𝑁2 ) 2 = 𝐼 and (𝑃𝜋𝑚2 ⊗ 𝐼 𝑁1 ) 2 = 𝐼 gives e (Ω𝑁1 ⊗ 𝐼 𝑁2 ). F = (Ω𝑁2 ⊗ 𝐼 𝑁1 ) Πswap 𝐷 This is the two-factor four-step operator. Evaluating it on the input 𝑥 reproduces the right-hand side of (6): (F 𝑥) [𝑘 2, 𝑘 1 ] =
𝑁∑︁ 2 −1 𝑁 1 −1 ∑︁
𝑥 [ 𝑗1 𝑁 2 +𝑗2 ] 𝜔 𝑁−𝑘1 1 𝑗1 𝜔 𝑁−𝑘1 𝑗2 𝜔 𝑁−𝑘2 2 𝑗2 = 𝑦 [𝑘 1 +𝑘 2 𝑁 1 ].
𝑗2 =0 𝑗1 =0
Hence 𝑌 [𝑘 2, 𝑘 1 ] = 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ] for all 𝑘 1 ∈ [𝑁 1 ], 𝑘 2 ∈ [𝑁 2 ]. Finally, the column factor 𝑣 is untouched by Phase III (which transforms only the row factor) and equals 𝜋𝑚1 (𝑎) = 𝑘 1 . Thus 𝑌 [𝑘 2, 𝑘 1 ] resides in physical column 𝑣 = 𝑘 1 of the output 𝑁 2 × 𝑁 1 slab, whose owner under (3) is ⌊𝑘 1 𝑝/𝑁 1 ⌋ = owner1 (𝑘 2, 𝑘 1 ). □
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
4
Cache-Oblivious Local FFT and Fused Redistribution
PACO is implemented by two cache-oblivious building blocks. The first is a batched local FFT that carries recursive layout changes lazily, without materializing recursion-level transposes. The second is a fused redistribution that applies the top-level twiddles, materializes the complete digit reversal, and changes ownership in one BSP communication round.
4.1
Cache-Oblivious Batched Local FFT
A local PACO subproblem consists of 𝑠 independent transforms of length 𝑛 = 𝑏 ℓ . The transforms are represented by recursive digitfield views rather than copied or transposed subarrays. One batch field remains the least significant physical address field; varying this field accesses a contiguous memory run. This distinguished field is preserved whenever the transform recursion creates new batch coordinates. The scheduler uses the relative sizes of 𝑛 and 𝑠, but does not inspect either 𝑀 or 𝐵. If 𝑠 ≥ 𝑛, it recursively splits one batch digit, thereby partitioning the transform instances into 𝑏 data-disjoint subbatches. It first consumes batch digits introduced by previous transform splits and splits the distinguished contiguous batch field only when no such outer batch digit remains. If 𝑠 < 𝑛, the scheduler performs a balanced digit-aligned Cooley–Tukey split. Specifically, it writes 𝑛 = 𝑛 1𝑛 2 ,
𝑛 1 = 𝑏 ⌈ℓ/2⌉ ,
𝑛 2 = 𝑏 ⌊ℓ/2⌋ ,
Proposition 1 (Cache-oblivious batched local FFT). For 𝑠 transforms of length 𝑛 = 𝑏 ℓ , the local PACO scheduler performs 𝑊local (𝑛, 𝑠) = 𝑂 (𝑛𝑠 (1 + log 𝑛))
(11)
local work and uses 𝑂 (log(𝑛𝑠)) recursion-stack words. Under the tall-cache assumption 𝑀 = Ω(𝐵 2 ), and for an initial contiguous batch width at least 𝐵, it incurs 𝑛𝑠 𝑄 local (𝑛, 𝑠) = 𝑂 1 + log𝑀 𝑛 𝐵 cache misses.
(12)
Proof sketch. Fix a sufficiently small constant 𝛼 > 0, let 𝑀 ★ = 𝛼𝑀, and consider the maximal recursive views satisfying 𝑛𝑠 ≤ 𝑀 ★ . This is only an analytical cache-fitting frontier: the algorithm itself never tests this condition. Frontier geometry. A transform split preserves the working-set size of each child group. Indeed, if
𝑛1 ≤ 𝑏. 𝑛2 Writing a logical transform coordinate as 1≤
𝑢 ∈ [𝑛 1 ],
Thus, the scheduler realizes the lazy-layout operators of Lemma 2 on every transform in the batch.
ℓ = log𝑏 𝑛,
so that
𝑥 = 𝑢𝑛 2 + 𝑣,
The top-level twiddle factors of each recursive Cooley–Tukey step are evaluated from the logical digit coordinates represented by the view. They are fused into the first accesses of the second child-transform group. At constant-size leaves, PACO invokes a natural-order local DFT kernel. In output-layout mode, the resulting batched operator is 𝐼𝑠 ⊗ 𝑃𝜋ℓ Ω𝑛 ; in input-layout mode, it is 𝐼𝑠 ⊗ Ω𝑛 𝑃𝜋ℓ .
𝑣 ∈ [𝑛 2 ],
and denoting an existing batch coordinate by 𝛽 ∈ [𝑡], the first child group consists of 𝑠𝑛 2 transforms of length 𝑛 1 . Its transform coordinate is 𝑢, while 𝑣 becomes a new outer batch coordinate: its batched transforms are indexed by
𝑛 = 𝑛 1𝑛 2 , then its two child groups have parameters (𝑛 1, 𝑠𝑛 2 )
and
(𝑛 2, 𝑠𝑛 1 ),
and both have volume
(𝑣, 𝛽).
𝑛 1 (𝑠𝑛 2 ) = 𝑛 2 (𝑠𝑛 1 ) = 𝑛𝑠.
Thus, the 𝑛 2 -digit field precedes the pre-existing batch digits in the child view, but the least-significant contiguous batch field within 𝛽 is unchanged. Symmetrically, the second child group consists of 𝑠𝑛 1 transforms of length 𝑛 2 , with transform coordinate 𝑣 and batch coordinate
Consequently, a non-cache-fitting view can first enter the cachefitting region only through a batch split. Every nonroot maximal frontier view therefore satisfies 𝑀★ < 𝑛𝑠 ≤ 𝑀 ★, (13) 𝑏 because its parent has the same transform length and 𝑏 times as many batch instances. We next bound the transform length at a frontier reached after at least one transform split. A transform split is executed only when 𝑠 < 𝑛. For a balanced digit-aligned factorization 𝑛1 ≤ 𝑏, 𝑛 = 𝑛 1𝑛 2 , 1≤ 𝑛2 the first child has parameters (𝑛 1, 𝑠𝑛 2 ) =: (𝑛 ′, 𝑠 ′ ). Since 𝑠 < 𝑛 1𝑛 2 ,
(𝑢, 𝛽). The associated Cooley–Tukey twiddle 𝜔𝑛−𝑢𝑣 is evaluated from these logical coordinates and is multiplied when the second child group first accesses the corresponding intermediate value. Hence, transform splitting changes only the recursive view: it creates an outer batch field but materializes neither a transpose nor a separate twiddle pass.
𝑠𝑛 2 < 𝑛 1𝑛 22 ≤ 𝑏 2𝑛 31 .
Shina Guo, Weiguo Gao, and Yuan Tang
Symmetrically, the second child satisfies
cache misses in total.
𝑠𝑛 1 < 𝑏 2𝑛 32 .
4.2
Thus, after every transform split, and hence after any later batch split (a batch split only decreases 𝑠 without changing 𝑛), the reachable state satisfies the following invariant: 𝑠 ≤ 𝑏 2𝑛 3 .
(14)
Combining (13) and (14) gives 𝑀★ < 𝑛𝑠 ≤ 𝑏 2𝑛 4, 𝑏 and therefore every frontier reached after at least one transform split has 𝑛 = Ω(𝑀 1/4 ). (15) Transform-induced frontier expansion. Let 𝐾 be the maximum number of transform splits on a root-to-frontier path. If the active transform length before one split is 𝑛 𝑗 = 𝑏𝑑 𝑗 , then the larger child has length at most √︁ 𝑏 ⌈𝑑 𝑗 /2⌉ ≤ 𝑏𝑛 𝑗 . Induction gives 𝐾
𝐾
𝐾
𝑛𝐾 ≤ 𝑏 1−1/2 𝑛 1/2 < 𝑏 𝑛 1/2 . For 𝐾 ≥ 1, combining this inequality with (15) yields 𝐾
𝑐𝑏 𝑀 1/4 ≤ 𝑏 𝑛 1/2
for a constant 𝑐𝑏 > 0 depending only on 𝑏 and 𝛼. Hence, 2𝐾 = 𝑂 1 + log𝑀 𝑛 . A batch split partitions working-set mass among its children, whereas a transform split creates two child groups, each of the parent’s mass. Thus, the total mass represented at the cache-fitting frontier is at most 𝑂 𝑛𝑠 2𝐾 = 𝑂 𝑛𝑠 1 + log𝑀 𝑛 . Spatial locality at the frontier. A separate argument (Lemma 17 in Appendix B.6) shows that whenever the distinguished contiguous field is split, the tall-cache assumption forces the resulting run width √ to satisfy 𝑐 = Ω( 𝑀) = Ω(𝐵); if it is never split it retains its initial width ≥ 𝐵 by (5). Hence a frontier view of mass 𝑛𝑠 is covered by 𝑛𝑠/𝑐 runs of length Ω(𝐵) and touches 𝑂 (𝑛𝑠/𝐵) blocks, which remain resident during its depth-first descendant execution. Summing over the frontier yields 𝑄 local (𝑛, 𝑠) = 𝑂 𝑛𝑠 𝐵 (1 + log𝑀 𝑛) . The work and stack bounds follow from the fixed-radix Cooley– Tukey recurrence and the digit-aligned recursion depth. More details in Appendix B. □ For PACO Phase I, Proposition 1 applies with 𝑁2 𝑛 = 𝑁 1, 𝑠= . 𝑝 For Phase III, it applies with 𝑁1 𝑛 = 𝑁 2, 𝑠= . 𝑝 The two local phases therefore incur, per processor, 𝑁 1 + log𝑀 𝑁 𝑂 𝑝𝐵
Fused Redistribution
The middle stage materializes the full digit reversal 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎) identified in Sect. 3.3. It consists of local packing, one balanced BSP exchange, and local unpacking. Local pack. Processor 𝑞 owns the Phase-I slab 𝑁2 e 𝑐], + 𝑐 loc . 𝐴[𝑎, 𝑐 =𝑞 𝑝 Because e 𝑐] = 𝐴[𝜋𝑚1 (𝑎), 𝑐], 𝐴[𝑎, a row macro-element permutation by 𝜋𝑚1 places the logical firststage value 𝐴[𝑣, 𝑐] in packed row 𝑣 = 𝜋𝑚1 (𝑎). During the same scan, PACO forms 𝑍 [𝑣, 𝑐] = 𝐴[𝑣, 𝑐]𝜔 𝑁−𝑣𝑐 . One BSP exchange. Decompose the packed logical row coordinate as 𝑁1 𝑁1 + 𝑣 loc, 𝑞 ′ ∈ [𝑝], 𝑣 loc ∈ . 𝑣 = 𝑞′ 𝑝 𝑝 The processor 𝑞 ′ is the destination owner of this item. After row packing, all rows with a common 𝑞 ′ form one contiguous source– destination block. For every source–destination pair (𝑞, 𝑞 ′ ), this block has shape 𝑁1 𝑁2 × 𝑝 𝑝 and contains exactly (𝑁 1 /𝑝) · (𝑁 2 /𝑝) = 𝑁 /𝑝 2 items. PACO sends these blocks in one balanced All-to-All BSP exchange. Local unpack. At destination 𝑞 ′ , the block received from source 𝑞 is indexed by (𝑣 loc, 𝑐 loc ). PACO transposes this local rectangle, producing the coordinate order (𝑐 loc, 𝑣 loc ). The source identifier 𝑞 and the first coordinate reconstruct the complete former column coordinate 𝑁2 + 𝑐 loc . 𝑐 =𝑞 𝑝 Finally, PACO applies 𝜋𝑚2 as a digit reversal of the 𝑁 2 row macroelements, each of size 𝑁 1 /𝑝. Writing 𝑢 = 𝜋𝑚2 (𝑐), the item is placed at the final local coordinate (𝑢, 𝑣) = 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) . Equivalently, the resulting local array satisfies 𝑍b[𝑢, 𝑣] = 𝑍 [𝑣, 𝜋𝑚2 (𝑢)],
(16)
(17)
which is precisely the input layout consumed by the Phase-III local FFT. Macro-element digit reversal. The row and column digit reversals in the fused stage are not elementwise random-access permutations. During packing, each of the 𝑁 1 rows is treated as one
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
Phase 1 pack rows
Phase 2 all-to-all ℎ = 𝑁 /𝑝 · (1 − 1/𝑝 ) P0 blocks routed around
𝑎 ↦→ 𝜋𝑚1 (𝑎)
Phase 3b unpack cols
Phase 3a transpose
𝑐 ↦→ 𝜋𝑚2 (𝑐 )
00,00
00,01
00,10
00,11
00,00
00,01
00,10
00,11
00,00
00,01
00,10
00,11
00,00
10,00
01,00
11,00
00,00
10,00
01,00
11,00
01,00
01,01
01,10
01,11
10,00
10,01
10,10
10,11
10,00
10,01
10,10
10,11
00,01
10,01
01,01
11,01
00,10
10,10
01,10
11,10
10,00
10,01
10,10
10,11
01,00
01,01
01,10
01,11
01,00
01,01
01,10
01,11
00,10
10,10
01,10
11,10
00,01
10,01
01,01
11,01
11,00
11,01
11,10
11,11
11,00
11,01
11,10
11,11
11,00
11,01
11,10
11,11
00,11
10,11
01,11
11,11
00,11
10,11
01,11
11,11
initial
packed rows 𝑎 → 𝜋𝑚 1 (𝑎)
𝑁 1 × 𝑝2 (col slab) 𝑁
transposed (𝜋𝑚 1 (𝑎) ↔𝑐 )
received (row slab)
Cell = ( 𝑎, 𝑐 ) in binary; fill colour = owning processor
P0
output
𝑁 2 × 𝑝1 (col slab) 𝑁
P1
Figure 1: PARA-CO-DRP𝑏 : tracking logical indices (𝑎, 𝑐), say (01, 00) in blue, through one global exchange 𝑁 1 =𝑁 2 =4, 𝑝 =2, 𝑏 =21 =2 (digit reversal on 2 bits: 1 ↔ 2, 0, 3 fixed) blue = P0 , orange = P1 contiguous macro-element of length 𝑁 2 /𝑝; during unpacking, each of the 𝑁 2 rows is treated as one contiguous macro-element of length 𝑁 1 /𝑝. A digit-reversal two-cycle is implemented by sequentially swapping its two macro-elements, using one temporary item. Coordinate trace. For an item initially stored at source processor 𝑞 with physical Phase-I coordinate 𝑁2 (𝑎, 𝑐), 𝑐 =𝑞 + 𝑐 loc, 𝑝 the fused middle stage performs the following coordinate trace: pack
(𝑎, 𝑐) ↦−−−→ (𝑣, 𝑐), route
𝑣 = 𝜋𝑚1 (𝑎), 𝑁1 𝑣 = 𝑞′ + 𝑣 loc, 𝑝
(𝑣, 𝑐) ↦−−−−→ 𝑞 ′ ; 𝑣 loc, 𝑐 loc ,
transpose 𝑞; 𝑣 loc, 𝑐 loc − ↦ −−−−−−→ 𝑞; 𝑐 loc, 𝑣 loc , unpack 𝑞; 𝑐 loc, 𝑣 loc ↦−−−−−→ (𝑢, 𝑣),
(18)
Hence, (𝑎, 𝑐) ↦−→ 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) = 𝜋𝑚 (𝑎 ++ 𝑐).
(19)
Theorem 4 (Correctness and complexity of the fused global permutation). Under assumptions (2) and (5), the composition of local packing, one BSP redistribution, and local unpacking materializes the full permutation 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎).
(20)
Applied to the Phase-I physical output, it produces 𝑍b[𝑢, 𝑣] = 𝑍 [𝑣, 𝜋𝑚2 (𝑢)], which is precisely the input layout required by Phase III. Each processor performs 𝑁 arith 𝑊𝑞,twiddle =Θ 𝑝
𝑁 𝑝
(27)
Excluding the mandatory local input and output arrays, the auxiliary workspace is 𝑁 𝑂 2 + log 𝑁 (28) 𝑝 items or words per processor. Proof. The coordinate trace (18) maps every Phase-I physical coordinate (𝑎, 𝑐) to (𝑢, 𝑣) = 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) . It therefore materializes 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎).
(21)
Since the Phase-I physical layout satisfies e 𝑐] = 𝐴[𝜋𝑚1 (𝑎), 𝑐], 𝐴[𝑎,
(22)
the packed row 𝑣 = 𝜋𝑚1 (𝑎) contains 𝐴[𝑣, 𝑐]. Fusing the top-level twiddle multiplication during this scan forms 𝑍 [𝑣, 𝑐] = 𝐴[𝑣, 𝑐]𝜔 𝑁−𝑣𝑐 .
twiddle arithmetic and rearr 𝑊𝑞,perm =𝑂
(24)
The permutation uses one global redistribution round with exactly 𝑁 (25) 𝑝2 items for every source–destination processor pair. Excluding the self block, every processor sends and receives 𝑁 1 ℎ= 1− (26) 𝑝 𝑝 remote items, and the aggregate remote migration volume is 1 𝜇 =𝑁 1− . 𝑝
at 𝑞 ′ (source = 𝑞) 𝑢 = 𝜋𝑚2 (𝑐).
local rearrangement work. Its local cache complexity is 𝑁 . 𝑄𝑞,perm = 𝑂 𝑝𝐵
The final coordinate mapping consequently gives
(23)
𝑍b[𝑢, 𝑣] = 𝑍 [𝑣, 𝜋𝑚2 (𝑢)],
Shina Guo, Weiguo Gao, and Yuan Tang
which is the Phase-III input layout. Packing applies macro-element digit reversal to 𝑁 1 rows of length 𝑁 2 /𝑝, and unpacking applies it to 𝑁 2 rows of length 𝑁 1 /𝑝. By the slab-width assumptions (5), both macro-element widths are at least 𝐵. Thus, together they require 𝑁 𝑂 𝑝 local rearrangement work and
𝑁 𝑂 𝑝𝐵
remote items sent and received by each processor. Its aggregate remote communication volume is 1 𝜇 =𝑁 1− . (32) 𝑝 Proof sketch. Consider an item initially stored at source coordinate (𝑎, 𝑐). The packing permutation sends it to logical first-stage frequency 𝑣 = 𝜋𝑚1 (𝑎),
cache misses; the fused twiddle scan contributes Θ(𝑁 /𝑝) arithmetic work without an additional cache pass. Each processor transposes 𝑝 blocks of size 𝑁 /𝑝 2 . The cacheoblivious transpose bound gives 𝑁 𝑂 𝑝+ 𝑝𝐵
thereby cancelling the output layout produced by Phase I; the same scan applies the twiddle 𝜔 𝑁−𝑣𝑐 . Writing 𝑣 = 𝑞′
routes the item to its unique Phase-III owner 𝑞 ′ . At that destination, transposing the received block exchanges (𝑣 loc, 𝑐 loc ) into (𝑐 loc, 𝑣 loc ), while the source identifier reconstructs
misses in total. Since
𝑐 =𝑞 𝑁 𝑁1 𝑁2 = ≥ 𝐵 2, 𝑝2 𝑝 𝑝
𝑁1 + 𝑣 loc 𝑝
𝑁2 + 𝑐 loc . 𝑝
The final row digit reversal places the item at
we have 𝑝 = 𝑂 (𝑁 /(𝑝𝐵)), so this cost is also 𝑂 (𝑁 /(𝑝𝐵)). For every source–destination pair, the block has
𝑢 = 𝜋𝑚2 (𝑐). Thus, the complete coordinate mapping is
𝑁1 𝑁2 𝑁 = 2 𝑝 𝑝 𝑝
(𝑎, 𝑐) ↦−→ 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) ,
items. Hence, excluding the self block, each processor sends and receives 𝑁 1 1− 𝑝 𝑝 remote items, and the aggregate remote volume is 1 𝑁 1− . 𝑝 The raw incoming regions are part of the mandatory output array; the only non-output array workspace is one reusable 𝑁 /𝑝 2 -item transpose buffer, plus 𝑂 (log 𝑁 ) recursion-stack words. □ Theorem 5 (Fused redistribution). The fused middle stage materializes the complete permutation 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎),
which is the full digit reversal 𝜋𝑚 (𝑎 ++ 𝑐). Both macro-element digit reversals and the local rectangular transpositions [8] are cache-oblivious linear-work, linear-cachemiss rearrangements: The row and column digit reversals are not element-wise random-access permutations. They are conducted on contiguous macro-elements of length 𝑁 2 /𝑝 and 𝑁 1 /𝑝, respectively. Given the slab-width conditions (5), these permutations take linear local work and 𝑂 (𝑁 /(𝑝𝐵)) cache misses. Since 𝜋𝑚1 is a bijection, each source has exactly 𝑁 1 /𝑝 packed rows for every destination 𝑞 ′ , and hence exactly 𝑁 /𝑝 2 items in every source–destination block. This gives the stated balanced communication volume and local complexity bounds. The complete in-place digit-reversal construction, cache analysis, and workspace details are deferred to the appendix. □
5
applies the top-level Cooley–Tukey twiddle factors, and produces the input layout consumed by Phase III. Per processor, it performs 𝑁 𝑊perm,local = Θ (29) 𝑝
Analysis and Optimality
This section combines the two local FFT phases and the fused middle redistribution, and then states the matching lower bounds under their respective scopes. Correctness was established in Theorem 3; here we analyze only local work, cache traffic, and interprocessor communication.
local work and incurs 𝑄 perm,local = Θ
𝑁 𝑝𝐵
(30)
cache misses. The unique communication round is a balanced BSP ℎ-relation with 𝑁 1 ℎ= 1− (31) 𝑝 𝑝
5.1
Upper Bounds
Table 1 summarizes the per-processor costs. In Phase I, every processor evaluates 𝑁 2 /𝑝 transforms of length 𝑁 1 . In Phase III, it evaluates 𝑁 1 /𝑝 transforms of length 𝑁 2 . The middle stage performs one fused redistribution: its top-level twiddles are included in local packing, and its local rearrangements consist of two macro-element digit reversals and rectangular block transpositions.
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
Table 1: Per-processor costs of the three PACO components.
5.2
Component
We now prove the lower bounds under the scopes stated after Theorem 1 and the fixed-radix scope 1 (see Remark 1 in Appendix D).
Local work Cache misses Communication 𝑁 𝑁 Local Phase I 𝑂 log 𝑁 1 𝑂 (1 + log𝑀 𝑁 1 ) 0 𝑝 𝑝𝐵 𝑁 𝑁 Fused middle stage 𝑂 𝑂 one ℎ-relation 𝑝 𝑝𝐵 𝑁 𝑁 Local Phase III 𝑂 log 𝑁 2 𝑂 (1 + log𝑀 𝑁 2 ) 0 𝑝 𝑝𝐵
Proposition 2 (PACO upper bounds). Under assumptions (2)– (5), PACO has maximum per-processor local work 𝑁 𝑊max = 𝑂 log 𝑁 (33) 𝑝 and maximum per-processor cache complexity 𝑁 1 + log𝑀 𝑁 . 𝑄 max = 𝑂 𝑝𝐵
(34)
𝑅 ≥ 1. Proof. Every DFT coefficient depends on every input item, since every entry of Ω𝑁 is nonzero. No processor initially owns the complete input. If 𝑅 = 0, then a processor’s local state and every output it returns are independent of some remotely owned input item, which contradicts correctness after changing only that input. □
𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎) between the source ownership (1) and target ownership (3). Any such algorithm must communicate at least 1 𝜇 ≥𝑁 1− (37) 𝑝
𝑅 = 1, with 1 𝑁 1− 𝑝 𝑝
Proposition 3 (At least one redistribution is necessary). Every correct parallel DFT algorithm satisfying the owner-computes and no-replication assumptions of Sect. 2 requires
Proposition 4 (Exact migration bound for the fused permutation). Consider materializing
It performs exactly one global redistribution round,
ℎ=
Lower Bounds and Exact Scope
(35)
remote items sent and received by each processor, and aggregate remote volume 1 . (36) 𝜇 =𝑁 1− 𝑝 Proof. Proposition 1 gives the local work and cache bounds for Phases I and III, while Theorem 5 gives the linear local cost of the middle stage. Since log 𝑁 1 + log 𝑁 2 = log 𝑁 , their local work bounds sum to (33); the same identity yields (34) after absorbing the constant number of linear scans. For every source–destination processor pair, the fused permutation routes exactly 𝑁 𝑝2 items. Excluding the self block, every processor therefore sends and receives 𝑁 𝑁 1 (𝑝 − 1) 2 = 1− 𝑝 𝑝 𝑝 remote items. Summing these balanced send volumes gives (36). □ The corresponding aggregate work and cache bounds are 𝑁 𝑂 (𝑁 log 𝑁 ) and 𝑂 1 + log𝑀 𝑁 , 𝐵 respectively. The local FFT recursion uses 𝑂 (log 𝑁 ) stack words, and the fused middle stage uses one reusable 𝑁 /𝑝 2 -item transpose scratch array; these workspace details are established in Section 4.2 and the appendix.
remote items. Moreover, 1 𝑁 max ℎ𝑞send, ℎ𝑞recv ≥ 1− . 𝑝 𝑝 𝑞 ∈ [𝑝 ]
(38)
Proof. For each fixed source column 𝑐, exactly 𝑁 1 /𝑝 row coordinates 𝑎 satisfy 𝜋𝑚1 (𝑎) 𝑐 = , 𝑁 1 /𝑝 𝑁 2 /𝑝 because 𝜋𝑚1 is a bijection. Hence exactly 𝑁2
𝑁1 𝑁 = 𝑝 𝑝
items retain their owner, while all remaining 𝑁−
𝑁 𝑝
items must cross a processor boundary at least once. Averaging the resulting aggregate send and receive volumes over the 𝑝 processors gives (38) since the maximum over 𝑝 processors is at least the average. □ Proposition 5 (Maximum Per-Processor work lower bound). For the fixed-base-𝑏 Cooley–Tukey DAG used by PACO, under the no-recomputation assumption, 𝑁 log 𝑁 . 𝑊max = Ω 𝑝 1 Let 𝑏 = 2𝑡 , 𝑡 = 𝑂 (1) . If every radix-𝑏 butterfly is expanded into its 𝑡 constituent
radix-2 layers, the resulting dependency graph is the standard radix-2 FFT DAG. Hence the maximum per-processor work/cache lower bounds are proved based on radix-2 DAG.
Shina Guo, Weiguo Gao, and Yuan Tang
Table 2: Representative FFT frameworks. “Communication” counts global exchanges; additional interprocessor communication is noted where relevant. Method
Numerical
Communication
Output ownership / cache guarantee
Four-/six-step FFT Frigo et al. CO-FFT Edelman et al. Tang et al. PACO
Exact Exact Tunable accuracy Tunable accuracy Exact
One or more transpose-like exchanges N/A (Sequential Alg.) One transpose + hierarchical communication One all-to-all + neighbor exchange One BSP redistribution only
Factor-swapped view possible; no cache-oblivious guarantee Natural-order output; optimal ideal-cache complexity Typically block ownership; no local ideal-cache theorem Typically block ownership; no local ideal-cache theorem Canonical output under factor-swapped slabs; cache-oblivious
Proof. The fixed-base FFT DAG contains Θ(𝑁 log 𝑁 ) arithmetic vertices. Without recomputation, these vertices must be executed somewhere, so at least one processor executes at least their average number, 𝑁 Ω log 𝑁 . 𝑝 □
compulsory misses, which is the required bound in this case. The detailed segment argument is deferred to Appendix. □ Together, Propositions 2– 6 establish Theorem 1 under the scopes stated following that theorem. As a byproduct, Appendix D gives a finite block-aware lower bound for arbitrary no-recomputation schedules of the radix-2 FFT DAG. In the stated asymptotic regime, it yields
Proposition 6 (Maximum Per-Processor cache lower bound). Under the exact slab decomposition and no arithmetic recomputation, the two local FFT phases satisfy 𝑁 𝑄 max = Ω 1 + log𝑀 𝑁 . (39) 𝑝𝐵
This stronger general statement is not needed for the PACO upperbound analysis and is therefore kept out of the main proof flow.
Consequently, the complete PACO execution satisfies the same lower bound.
5.3
Proof roadmap. Each processor evaluates a batch of 𝑁 2 /𝑝 length𝑁 1 FFT DAGs in Phase I and a batch of 𝑁 1 /𝑝 length-𝑁 2 FFT DAGs in Phase III. Applying an FFT-DAG input-boundary lower bound [15] to these data-disjoint batches yields 𝑁 Ω 1 + log𝑀 𝑁 1 𝑝𝐵 and
𝑁 1 + log𝑀 𝑁 2 𝑝𝐵 misses, respectively, up to the standard initial-cache term. Fix a sufficiently large constant 𝜖. If 𝑁 ≥ 𝜖𝑀, 𝑝 Ω
the Θ(𝑀/𝐵) boundary term is absorbed. Otherwise, 𝑁 < 𝜖𝑀. 𝑝 Because 𝑝 | 𝑁 1 and 𝑝 | 𝑁 2 , 𝑁2 ≥ 1, 𝑝 Hence 𝑁2 𝑁 𝑁1 ≤ 𝑁1 = , 𝑝 𝑝 Therefore,
𝑁1 ≥ 1. 𝑝 𝑁2 ≤ 𝑁2
𝑁1 𝑁 = . 𝑝 𝑝
2 𝑁 < 𝜖 2 𝑀 2, 𝑝 and thus 1 + log𝑀 𝑁 = 𝑂 (1). Since Phase I begins with an empty cache and accesses Θ(𝑁 /𝑝) distinct local input items, it incurs 𝑁 Ω 𝑝𝐵
𝑁 = 𝑁1 𝑁2 ≤
𝑄 miss ≥ (1 − 𝑜 (1))
𝑁 log𝑀 𝑁 . 𝐵
Related Work
Classical four-step and six-step FFTs decompose a transform into local FFT groups separated by twiddle factors and transpose-like changes of view [2, 13]. Exact distributed-memory FFTs have further explored the tradeoffs among data allocation, input/output ordering, and communication volume [10, 11]. These works show that a factor-swapped or otherwise reordered output can save a transpose. PACO adopts an explicit factor-swapped slab contract, but additionally proves that its local computation is fully cache oblivious and that its sole redistribution materializes the layout deferred by the local recursions. Cache-oblivious FFTs obtain ideal-cache locality through recursive decomposition and recursive layout transformations [8]. PACO uses the same cache-oblivious objective in a distributed setting, but does not materialize a transpose-like layout transformation at each recursive split. Instead, it represents these transformations lazily, shows that they compose into base-𝑏 digit reversal, and fuses that permutation with the ownership change between the two top-level FFT groups. Several lines of work reduce communication under different numerical, distributional, or practical contracts. Approximate distributed transforms use low-rank structure or oversampling to reduce global communication [7, 18]. Exact cyclic-to-cyclic multidimensional FFTs can also use a single all-to-all exchange [12]; unlike PACO, they retain cyclic ownership rather than returning a factorswapped slab distribution. Other multidimensional FFT frameworks optimize decomposition choices, transpose order, communication overlap, or MPI datatype-based redistribution [1, 5, 6, 14, 17]. They target practical multidimensional performance and do not establish PACO’s processor-local ideal-cache bound or its fused digit-reversalpermutation guarantee. Finally, Hong and Kung introduced the red-blue pebble framework for FFT I/O lower bounds [9], while Ranjan, Savage, and
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
Zubair established FFT-DAG boundary results used by our localcache lower bound [15]. BSP and distributed-memory lower-bound work also gives communication lower bounds for FFT computation under broader schedule classes [3, 16]. PACO’s claims are more specific: one redistribution is necessary under its owner-computes and no-replication assumptions (Proposition 3), and its migration bound is scoped as in the Scope of optimality paragraph of Sect. 2.
5.4
Discussion: Ownership, Permutations, and Composition
PACO’s factor-swapped output ownership is an explicit interface contract, rather than an omitted cost. The sole redistribution changes from input-column slabs indexed by the 𝑁 2 field to output-column slabs indexed by the 𝑁 1 field, exactly as required to make every second-phase length-𝑁 2 transform local. This contract can nevertheless be useful when subsequent computation can consume the factor-swapped frequency view directly. More importantly for the present result, it separates two issues that are often conflated: the global ownership change required by the Cooley–Tukey dependency structure, and the local layout changes introduced by cache-oblivious recursion. PACO shows that, for the stated slab decomposition, the latter need not cause an additional redistribution: their complete base-𝑏 digit reversal can be fused into the former. The resulting middle stage remains a balanced all-to-all, with every source–destination pair carrying exactly 𝑁 /𝑝 2 items.
5.5
Limitations
PACO assumes a fixed power-of-two radix and exact divisibility: 𝑁 = 𝑁 1 𝑁 2 = 𝑏𝑚 , with 𝑝, 𝑁 1 , and 𝑁 2 aligned to the base-𝑏 slab decomposition. Arbitrary sizes and processor counts require an imbalance analysis. The algorithm returns canonical coefficient indices under factor-swapped ownership; restoring the original ownership may require another redistribution. The present analysis is limited to one-dimensional slab decompositions, rather than pencil or block distributions. Finally, it assumes homogeneous idealcache and BSP parameters and excludes arithmetic recomputation from its work and cache optimality claims.
5.6
Conclusion
PACO shows that cache-oblivious local FFT execution and a oneredistribution parallel organization can coexist when recursive layout changes are represented lazily and fused with the unavoidable ownership change. Under the stated contract, it achieves optimal maximum per-processor work and cache complexity, uses the minimum number of global redistribution rounds, and attains the exact migration cost of its prescribed fused permutation. Extending these guarantees to more general ownership distributions and processor configurations remains open.
Shina Guo, Weiguo Gao, and Yuan Tang
References [1] Alan Ayala, Stanimire Tomov, Azzam Haidar, and Jack Dongarra. 2020. heFFTe: Highly Efficient FFT for Exascale. In Computational Science – ICCS 2020 (Lecture Notes in Computer Science, Vol. 12137). Springer, 262–275. doi:10.1007/978-3-03050371-0_19 [2] David H. Bailey. 1990. FFTs in External or Hierarchical Memory. The Journal of Supercomputing 4, 1 (1990), 23–35. doi:10.1007/BF00162341 [3] Gianfranco Bilardi, Michele Scquizzato, and Francesco Silvestri. 2018. A Lower Bound Technique for Communication in BSP. ACM Transactions on Parallel Computing 4, 3, Article 14 (2018), 27 pages. doi:10.1145/3181776 [4] James W. Cooley and John W. Tukey. 1965. An Algorithm for the Machine Calculation of Complex Fourier Series. Math. Comp. 19, 90 (1965), 297–301. doi:10.1090/S0025-5718-1965-0178586-1 [5] Lisandro Dalcin, Mikael Mortensen, and David E. Keyes. 2019. Fast Parallel Multidimensional FFT Using Advanced MPI. J. Parallel and Distrib. Comput. 128 (2019), 137–150. doi:10.1016/j.jpdc.2019.02.006 [6] Truong Vinh Truong Duy and Taisuke Ozaki. 2014. A Decomposition Method with Minimum Communication Amount for Parallelization of Multi-Dimensional FFTs. Computer Physics Communications 185, 1 (2014), 153–164. doi:10.1016/j. cpc.2013.08.028 [7] Alan Edelman, Peter McCorquodale, and Sivan Toledo. 1999. The Future Fast Fourier Transform? SIAM Journal on Scientific Computing 20, 3 (1999), 1094–1114. doi:10.1137/S1064827597316266 [8] Matteo Frigo, Charles E. Leiserson, Harald Prokop, and Sridhar Ramachandran. 1999. Cache-Oblivious Algorithms. In Proceedings of the 40th Annual Symposium on Foundations of Computer Science. IEEE Computer Society, 285–297. doi:10. 1109/SFFCS.1999.814600 [9] Jia-Wei Hong and H. T. Kung. 1981. I/O Complexity: The Red-Blue Pebble Game. In Proceedings of the 13th Annual ACM Symposium on Theory of Computing (STOC ’81). Association for Computing Machinery, 326–333. doi:10.1145/800076.802486
[10] Márcia A. Inda and Rob H. Bisseling. 2001. A Simple and Efficient Parallel FFT Algorithm Using the BSP Model. Parallel Comput. 27, 14 (2001), 1847–1878. doi:10.1016/S0167-8191(01)00118-1 [11] S. Lennart Johnsson, Michel Jacquemin, and Robert L. Krawitz. 1992. Communication Efficient Multi-Processor FFT. J. Comput. Phys. 102, 2 (1992), 381–397. doi:10.1016/0021-9991(92)90380-H [12] Thomas Koopman and Rob H. Bisseling. 2023. Minimizing Communication in the Multidimensional FFT. SIAM Journal on Scientific Computing 45, 6 (2023), C330–C347. doi:10.1137/22M1487242 [13] Charles F. Van Loan. 1992. Computational Frameworks for the Fast Fourier Transform. Frontiers in Applied Mathematics, Vol. 10. Society for Industrial and Applied Mathematics, Philadelphia, PA. doi:10.1137/1.9781611970999 [14] Dmitry Pekurovsky. 2012. P3DFFT: A Framework for Parallel Computations of Fourier Transforms in Three Dimensions. SIAM Journal on Scientific Computing 34, 4 (2012), C192–C209. doi:10.1137/11082748X [15] Desh Ranjan, John E. Savage, and Mohammad Zubair. 2011. Strong I/O Lower Bounds for Binomial and FFT Computation Graphs. In Computing and Combinatorics: 17th Annual International Conference, COCOON 2011 (Lecture Notes in Computer Science). Springer. doi:10.1007/978-3-642-22685-4_12 [16] Michele Scquizzato and Francesco Silvestri. 2014. Communication Lower Bounds for Distributed-Memory Computations. In 31st International Symposium on Theoretical Aspects of Computer Science (STACS 2014) (Leibniz International Proceedings in Informatics (LIPIcs), Vol. 25). Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 627–638. doi:10.4230/LIPIcs.STACS.2014.627 [17] Sukhyun Song and Jeffrey K. Hollingsworth. 2014. Designing and Auto-Tuning Parallel 3-D FFT for Computation-Communication Overlap. In Proceedings of the 19th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. ACM, 181–192. doi:10.1145/2692916.2555249 [18] Ping Tak Peter Tang, Jongsoo Park, Daehyun Kim, and Vladimir Petrov. 2012. A Framework for Low-Communication 1-D FFT. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’12). IEEE Computer Society, 1–11. doi:10.1109/SC.2012.5
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
A
Deferred Layout Algebra and Local FFT Correctness
This appendix formalizes the deferred-layout interpretation used by the local PACO FFT recursion. Its purpose is to prove that recursive transform splits do not materialize physical transposes, and that the resulting lazy layouts are precisely base-𝑏 digit reversals. In particular, we establish the two local operator identities used in Sect. 3: Aℓ = 𝑃𝜋ℓ Ω𝑏 ℓ and Bℓ = Ω𝑏 ℓ 𝑃𝜋ℓ . The global materialization and redistribution of these layouts are handled separately in Appendix B.
A.1
Coordinate, Digit-Field, and Layout Conventions
Fix the constant radix 𝑏 = 2𝑡 . An ℓ-digit base-𝑏 string is written from most significant to least significant digit, 𝑢 = 𝑑 ℓ −1 ++ 𝑑 ℓ −2 ++ · · · ++ 𝑑 0,
𝑑𝑟 ∈ [𝑏],
where ++ denotes digit-field concatenation: if 𝑈 and 𝑉 have lengths ℓ𝑈 and ℓ𝑉 , then 𝑈 ++ 𝑉 is the (ℓ𝑈 + ℓ𝑉 )-digit field whose digits in 𝑈 are more significant than those in 𝑉 . For 0 ≤ ℓ, the base-𝑏 digit-reversal permutation 𝜋ℓ : [𝑏 ℓ ] → [𝑏 ℓ ] is 𝜋ℓ (𝑑 ℓ −1 ++ · · · ++ 𝑑 0 ) = 𝑑 0 ++ · · · ++ 𝑑 ℓ −1, an involution: 𝜋ℓ−1 = 𝜋ℓ and 𝜋ℓ2 = id. For a permutation 𝜋 : [𝑛] → [𝑛], the physical-to-logical layout operator is (𝑃𝜋 𝑧) [𝑖] = 𝑧 [𝜋 (𝑖)]; thus, under layout 𝑃𝜋 , the value at physical position 𝑖 represents logical coordinate 𝜋 (𝑖). Since digit reversal is an involution, 𝑃𝜋2ℓ = 𝐼 .
(40)
We distinguish throughout: • Physical digit slots: fixed base-𝑏 positions in the physical address of an item. • Logical transform-input digits: the input coordinate of an active local transform. • Logical output-frequency digits: a transform frequency coordinate after a child transform has completed. • Batch digits: distinguish independent transform instances. • Fixed digits: selected by a recursive batch branch and hence holding one fixed value within that branch. A transform split may change the logical role of a physical digit slot (e.g. from a transform-input digit to an outer-batch digit); it never moves the item stored at that physical address.
A.2
Layout-Aware Twiddle Operator
Consider one local two-factor Cooley–Tukey split 𝑛 = 𝑛 1𝑛 2 , 𝑛 1 = 𝑏 ℓ1 , 𝑛 2 = 𝑏 ℓ2 , ℓ = ℓ1 + ℓ2 , with 𝑗 = 𝑗1𝑛 2 + 𝑗2 and 𝑘 = 𝑘 1 + 𝑘 2𝑛 1 . Suppose that after the first child-transform group the physical coordinate 𝑎 ∈ [𝑛 1 ] represents the logical first-stage frequency 𝑘 1 = 𝜋ℓ1 (𝑎). The twiddle multiplication must use this represented logical coordinate rather than 𝑎: −𝜋ℓ (𝑎) 𝑗2 view 𝐷 tw 𝑧 [𝑎, 𝑗2 ] = 𝜔𝑛 1 𝑧 [𝑎, 𝑗2 ]. (41)
Equivalently, the ordinary diagonal twiddle operator is conjugated by the current layout permutation and evaluated lazily from the view descriptor; no physical permutation is needed. The correctness of this lazy evaluation is established in Appendix A.7.
A.3
Admissible Digit-Field Views
Let 𝐷 = {𝑑 0, . . . , 𝑑𝐷 −1 } be the digit labels of the physical base-𝑏 slots of a root local array (slot 0 least significant). A recursive view is a tuple 𝑉 = (𝑇 , 𝑂, 𝐶, 𝐹, 𝑃, mode): • 𝑇 : ordered active transform digit field; • 𝑂: ordered list of active outer-batch digit fields; • 𝐶: distinguished contiguous batch field; • 𝐹 : partial assignment from fixed digit labels to values in [𝑏]; • 𝑃: bijection from current digit labels to physical digit slots; • mode ∈ {OutputLayout, InputLayout}. The active labels in𝑇 , 𝑂, 𝐶 and the fixed labels in dom(𝐹 ) are disjoint and cover the root digit labels. A label may be renamed when its logical role changes, but its physical slot under 𝑃 is unchanged. Let 𝑛 = 𝑏 |𝑇 | ,
𝑜 = 𝑏 |𝑂 | ,
𝑐 = 𝑏 |𝐶 | ,
𝑠 = 𝑜𝑐,
so the view represents 𝑠 transforms of length 𝑛; a logical item has coordinate (𝜏, 𝜂, 𝜒) ∈ [𝑛] × [𝑜] × [𝑐]. For a digit label 𝑑, define stride𝑃 (𝑑) = 𝑏 pos𝑃 (𝑑 ) . The active digits of 𝐶 occupy the least-significant contiguous physical interval, slots 0, . . . , |𝐶 | − 1 (fixed digits may occupy immediately higher slots without interrupting the remaining active interval). The address map is ∑︁ addr𝑉 (𝜏, 𝜂, 𝜒) = base0 + 𝐹 (𝑑) stride𝑃 (𝑑) 𝑑 ∈dom(𝐹 )
+
∑︁
digit𝑑 (𝜏) stride𝑃 (𝑑)
𝑑 ∈𝑇
+
∑︁
digit𝑑 (𝜂) stride𝑃 (𝑑) + 𝜒 .
(42)
𝑑 ∈𝑂
For fixed (𝜏, 𝜂, 𝐹 ), set base𝜏,𝜂,𝐹 = addr𝑉 (𝜏, 𝜂, 0); then addr𝑉 (𝜏, 𝜂, 𝜒) = base𝜏,𝜂,𝐹 + 𝜒, so varying 𝜒 traverses a contiguous run. Definition 6 (Admissible digit-field view). A view is admissible if: (1) its active and fixed labels partition the root digit slots; (2) 𝑃 assigns each current label a unique physical slot; (3) the active digits of 𝐶 form one least-significant contiguous interval; (4) its address map is (42); and (5) every logical coordinate maps to a distinct physical item. Lemma 7 (Active-run decomposition). An admissible view with 𝑠 transforms of length 𝑛 is the disjoint union of 𝑛𝑠/𝑐 contiguous physical runs, each of length 𝑐. Proof. There is one run per pair (𝜏, 𝜂) ∈ [𝑛] × [𝑜]; varying 𝜒 ∈ [𝑐] gives a consecutive interval of length 𝑐. There are 𝑛𝑜 = 𝑛𝑠/𝑐 such pairs, and injectivity of the address map makes the intervals disjoint and covering. □ The distinguished field 𝐶 is retained to expose spatial locality even when the active transform field is strided; this is used by the cache analysis in Appendix C.
Shina Guo, Weiguo Gao, and Yuan Tang
A.4
Legal View Operations
The scheduler uses outer-batch splits, contiguous-field splits, and transform splits. Batch splits partition batch instances into datadisjoint children; a transform split creates two sequential child groups reinterpreting the same parent working set. Outer-batch split. If 𝑑 ∈ 𝑂 is selected, the 𝑟 -th child (𝑟 ∈ [𝑏]) has 𝑇𝑟 = 𝑇 , 𝑂𝑟 = 𝑂 \ {𝑑 }, 𝐶𝑟 = 𝐶, 𝐹𝑟 = 𝐹 ∪ {𝑑 ↦→ 𝑟 }; 𝑃 and mode are unchanged. The 𝑏 children are data disjoint. Contiguous-field split. Used only when 𝑂 = ∅. Write 𝐶 = 𝑑 ++ 𝐶 ′ with 𝑑 the most significant active digit of 𝐶; the 𝑟 -th child has 𝑇𝑟 = 𝑇 , 𝑂𝑟 = ∅, 𝐶𝑟 = 𝐶 ′ , 𝐹𝑟 = 𝐹 ∪ {𝑑 ↦→ 𝑟 }. The child width is 𝑐/𝑏, and 𝐶 ′ still occupies the least-significant contiguous interval. Transform split. Decompose 𝑇 = 𝑇1 ++ 𝑇2 with ℓ1 = |𝑇1 |, ℓ2 = |𝑇2 |, 𝑛 1 = 𝑏 ℓ1 , 𝑛 2 = 𝑏 ℓ2 , 𝑛 = 𝑛 1𝑛 2 , and write 𝑗 = 𝑗1𝑛 2 + 𝑗2 . The first child group performs length-𝑛 1 transforms over 𝑗1 : 𝑇I = 𝑇1,
𝑂 I = 𝑇2 ++ 𝑂,
𝐶 I = 𝐶,
with 𝑠𝑛 2 transforms and inherited address map addr𝑉I 𝑗1, ( 𝑗2, 𝜂), 𝜒 = addr𝑉 𝑗1𝑛 2 + 𝑗2, 𝜂, 𝜒 .
(43) (44)
After the first child completes, the slots formerly holding 𝑇1 contain first-stage outputs; denote this physical field 𝐴. This is a semantic relabeling only: pos𝑃 (𝐴 𝑗 ) = pos𝑃 (𝑇1 ) 𝑗 . (45) In output-layout mode a physical coordinate 𝑎 ∈ [𝑛 1 ] of 𝐴 represents 𝑘 1 = 𝜋ℓ1 (𝑎). The second child group performs length-𝑛 2 transforms over 𝑇2 : 𝑇II = 𝑇2,
𝑂 II = 𝐴 ++ 𝑂,
𝐶 II = 𝐶,
with 𝑠𝑛 1 transforms and address map addr𝑉II 𝑗2, (𝑎, 𝜂), 𝜒 = addr𝑉 𝑎𝑛 2 + 𝑗2, 𝜂, 𝜒 .
(46) (47)
The field 𝑇2 may be strided past the physical field 𝐴; 𝑃 records this interpretation, and no exchange of the two fields is performed. The −𝜋ℓ (𝑎) 𝑗2 twiddle 𝜔𝑛 1 is fused into the first access of the corresponding second-child input, as in (41).
run into 𝑏 subruns. A transform split has child maps (44) and (47), each a bijective reindexing of the parent set; it does not change 𝐹 or the slots of 𝐶, and the relabeling 𝑇1 → 𝐴 preserves slots by (45). Thus every operation preserves Definition 6. □ The fixed assignment 𝐹 is necessary because batch and transform splits can interleave. For example, with physical digit order [𝑎 1 ] [𝑎 0 ] [𝑐 1 ] [𝑐 0 ] and 𝑇 = [𝑎 1, 𝑎 0 ], 𝑂 = ∅, 𝐶 = [𝑐 1, 𝑐 0 ], 𝐹 = ∅, a contiguous-field split fixing 𝑐 1 = 𝑑 gives 𝐶 = [𝑐 0 ], 𝐹 = {𝑐 1 ↦→ 𝑑 }; splitting [𝑎 1, 𝑎 0 ] = [𝑎 1 ] ++ [𝑎 0 ] makes the second child’s transform field [𝑎 0 ] with outer-batch field [𝐴1 ]; a later batch split fixing 𝐴1 = 𝑒 yields 𝑇 = [𝑎 0 ], 𝑂 = ∅, 𝐶 = [𝑐 0 ], 𝐹 = {𝑐 1 ↦→ 𝑑, 𝐴1 ↦→ 𝑒}. Although the fixed slots need not be adjacent, varying 𝑐 0 still traverses a contiguous run—exactly the information recorded by 𝐹 and 𝑃.
A.6
Theorem 9 (No physical transposes inside a local transform). No transform split in the local PACO recursion materializes a physical permutation. The transition from the first child group to the second changes only (1) the logical role of the physical 𝑇1 slots, renamed as 𝐴; (2) the active transform field, from 𝑇1 to 𝑇2 ; and (3) the descriptor 𝑃. All items remain in their original physical positions. Proof. The child maps (44) and (47) refer only to physical addresses of the parent. The first child overwrites each parent item with its first-stage result; relabeling 𝑇1 → 𝐴 preserves positions by (45); the second child accesses the inherited 𝑇2 slots through view-defined strides while treating 𝐴 as outer-batch. Hence no address-to-address copy, transpose, or recursive layout permutation occurs, and the argument applies recursively at every transformsplit node. □ The only physical digit permutation in the complete PACO execution is the fused middle permutation of Sect. 4.2 and Appendix C.5; the macro-element digit reversals used there are components of that single materialization, not recursive transposes.
A.7 A.5
Closure Under Arbitrary Split Interleavings
Theorem 8 (Closure of admissible views). Starting from an admissible root view, every finite sequence of legal outer-batch, contiguousfield, and transform splits produces admissible child views. Moreover: (1) active and fixed labels continue to partition the physical slots; (2) each label retains its physical slot under 𝑃, except for the semantic relabeling of a transform-input field as a transform-output field; (3) 𝐶 remains one least-significant contiguous interval; (4) the children of every batch split are data disjoint and partition the parent address set; (5) each transform child group is a bijective reindexing of the complete parent working set; and (6) invariant (42) is preserved under arbitrary interleavings. Proof. Induction on the number of operations; the base case is Definition 6. An outer-batch split moves one digit from 𝑂 to 𝐹 with a fixed value, leaving 𝑇 , 𝐶 and all physical slots unchanged, and its 𝑏 values partition the parent batch. A contiguous-field split moves the most significant digit of 𝐶 to 𝐹 ; the remaining active digits still form the least-significant interval, and the 𝑏 values partition every
Transform Splits Are View Transformations
Output-Layout Local FFT Correctness
Let a view contain 𝑠 data-disjoint transforms of length 𝑛 = 𝑏 ℓ = Í 𝑛 1𝑛 2 , and for batch instance 𝑞 let 𝑌 [𝑘, 𝑞] = 𝑗 𝑉 [ 𝑗, 𝑞] 𝜔𝑛−𝑘 𝑗 . Theorem 10 (Output-layout local FFT identity). For every ℓ ≥ 0 the output-layout traversal implements Aℓ = 𝑃𝜋ℓ Ω𝑏 ℓ ; equivalently, for a physical output coordinate 𝑟 , 𝑌 phys [𝑟, 𝑞] = 𝑌 [𝜋ℓ (𝑟 ), 𝑞]. Proof. Induction on ℓ. For ℓ = 0 there is no transform coordinate; for ℓ = 1 the routine performs one base-𝑏 DFT and 𝜋 1 = id. Let ℓ > 1 and split ℓ = ℓ1 + ℓ2 , 𝑗 = 𝑗1𝑛 2 + 𝑗2 , 𝑘 = 𝑘 1 + 𝑘 2𝑛 1 . The first child computes, for each ( 𝑗2, 𝑞), 𝐻 [𝑘 1, 𝑗2, 𝑞] =
𝑛∑︁ 1 −1
𝑉 [ 𝑗1𝑛 2 + 𝑗2, 𝑞] 𝜔𝑛−𝑘1 1 𝑗1 ,
(48)
𝑗1 =0
so by induction its physical coordinate 𝑎 represents 𝑘 1 = 𝜋ℓ1 (𝑎): 𝐻 phys [𝑎, 𝑗2, 𝑞] = 𝐻 [𝜋ℓ1 (𝑎), 𝑗2, 𝑞]. Applying the twiddle, 𝑍 [𝑘 1, 𝑗2, 𝑞] = 𝐻 [𝑘 1, 𝑗2, 𝑞] 𝜔𝑛−𝑘1 𝑗2 , gives 𝑍 phys [𝑎, 𝑗2, 𝑞] = 𝑍 [𝜋ℓ1 (𝑎), 𝑗2, 𝑞].
(49)
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
CO-Batched-FFT(𝑉 ) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
𝑛 = 𝑏 |𝑇 | 𝑠 = 𝑏 |𝑂 |+|𝐶 | if 𝑛 ≤ 𝑏 and 𝑠 ≤ 𝑏 // base case invoke the constant-size natural-order DFT kernels return if 𝑠 ≥ 𝑛 // batch split if 𝑂 ≠ ∅ choose an active outer-batch digit 𝑑 ∈ 𝑂 else choose the most significant active digit 𝑑 ∈ 𝐶 for 𝑟 = 0 to 𝑏 − 1 CO-Batched-FFT Batch-Child(𝑉 , 𝑑, 𝑟 ) return // transform split split 𝑇 = 𝑇1 ++ 𝑇2 with |𝑇1 | − |𝑇2 | ≤ 1 𝑉I = First-Transform-Child(𝑉 ,𝑇1,𝑇2 ) CO-Batched-FFT(𝑉I ) relabel the completed physical 𝑇1 -field as output field 𝐴 𝑉II = Second-Transform-Child(𝑉 ,𝑇1,𝑇2, 𝐴) fuse layout-aware twiddles into the first accesses of 𝑉II CO-Batched-FFT(𝑉II )
Proof. The input-layout traversal uses digit-reversed entrance and exit interpretations, so Bℓ = 𝑃𝜋ℓ Aℓ 𝑃𝜋ℓ . By Theorem 10 and (40), Bℓ = 𝑃𝜋ℓ (𝑃𝜋ℓ Ω𝑏 ℓ )𝑃𝜋ℓ = Ω𝑏 ℓ 𝑃𝜋ℓ , hence Bℓ 𝑥 phys = Ω𝑏 ℓ 𝑃𝜋ℓ 𝑃𝜋ℓ 𝑥 = Ω𝑏 ℓ 𝑥. The two occurrences of 𝑃𝜋ℓ are view interpretations, not runtime physical permutations. □
A.9
Local Correctness Summary
Corollary 12 (Batched local FFT operators). For 𝑠 independent transforms of length 𝑏 ℓ , the output- and input-layout traversals implement Aℓ,𝑠 = 𝐼𝑠 ⊗ (𝑃𝜋ℓ Ω𝑏 ℓ ) and Bℓ,𝑠 = 𝐼𝑠 ⊗ (Ω𝑏 ℓ 𝑃𝜋ℓ ), respectively. Proof. Theorem 8 preserves data disjointness across batch instances; applying Theorems 10 and 11 per instance gives the Kronecker forms. □ Consequently, Phase I invokes A𝑚1 (physical row 𝑎 represents 𝜋𝑚1 (𝑎)); the fused middle stage materializes 𝜋𝑚 (𝑎 ++𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎), converting the Phase-I output layout into the Phase-III input layout; and Phase III invokes B𝑚2 , producing natural-order coefficients as used in Theorem 3.
B
Cache-Oblivious Batched Local FFT
Figure 2: Cache-oblivious batched local FFT scheduler. The operations Batch-Child, First-Transform-Child, and Second-Transform-Child are the admissible view operations defined in Appendix A.
This appendix supplies the complete scheduler and the rigorous cache analysis underlying Proposition 1 of the main text. The scheduler operates on the admissible digit-field views defined in Appendix A. Consider an admissible view
The second child computes, for each (𝑘 1, 𝑞), 𝑌 [𝑘 1 + 𝑘 2𝑛 1, 𝑞] = Í −𝑘 2 𝑗2 . Substituting (48) and the twiddle, 𝑗2 𝑍 [𝑘 1 , 𝑗 2 , 𝑞] 𝜔𝑛 2
Let
𝑌 [𝑘 1 + 𝑘 2𝑛 1, 𝑞] =
𝑛∑︁ 2 −1 𝑛∑︁ 1 −1
𝑉 [ 𝑗1𝑛 2 + 𝑗2, 𝑞] 𝜔𝑛−𝑘1 1 𝑗1 𝜔𝑛−𝑘1 𝑗2 𝜔𝑛−𝑘2 2 𝑗2 .
𝑗2 =0 𝑗1 =0
𝑉 = (𝑇 , 𝑂, 𝐶, 𝐹, 𝑃, mode). 𝑛 = 𝑏 |𝑇 | , 𝑠 = 𝑏 |𝑂 |+|𝐶 | , 𝑐 = 𝑏 |𝐶 | , (52) where 𝑛 is the transform length, 𝑠 is the number of batched transforms, and 𝑐 is the width of the distinguished contiguous batch field. We call the product
Using 𝜔𝑛𝑛2 = 𝜔𝑛1 and 𝜔𝑛𝑛1 = 𝜔𝑛2 , the three factors combine into 𝜇 (𝑉 ) = 𝑛𝑠 𝜔𝑛− (𝑘1 𝑗1𝑛2 +𝑘1 𝑗2 +𝑘2𝑛1 𝑗2 ) . Adjoining the cross term 𝜔𝑛−𝑘2𝑛1 𝑗1𝑛2 = 𝜔𝑛−𝑘2 𝑗1 (𝑛1𝑛2 ) = 𝑛 −𝑘 𝑗 2 1 the mass of the view; it is the number of items in its working set. (𝜔𝑛 ) = 1 completes the exponent, so ∑︁ The scheduler is independent of the cache parameters 𝑀 and 𝐵. 𝑌 [𝑘 1 + 𝑘 2𝑛 1, 𝑞] = 𝑉 [ 𝑗1𝑛 2 + 𝑗2, 𝑞] 𝜔𝑛− (𝑘1 +𝑘2𝑛1 ) ( 𝑗1𝑛2 +𝑗2 ) , (50) 𝑗1 ,𝑗2
which, as 𝑗 = 𝑗1𝑛 2 + 𝑗2 ranges over [𝑛], is the length-𝑛 DFT at frequency 𝑘 1 + 𝑘 2𝑛 1 . By (49) the second child’s input for fixed physical 𝑎 is 𝑗2 ↦→ 𝑍 [𝜋ℓ1 (𝑎), 𝑗2, 𝑞]; by induction its physical output coordinate 𝑐 represents 𝑘 2 = 𝜋ℓ2 (𝑐), whence 𝑌 phys [𝑎, 𝑐, 𝑞] = 𝑌 𝜋ℓ1 (𝑎) + 𝜋ℓ2 (𝑐)𝑛 1, 𝑞 . (51) Finally, the parent physical coordinate 𝑟 = 𝑎𝑛 2 + 𝑐 has digit string (digits of 𝑎) ++ (digits of 𝑐), so 𝜋ℓ (𝑎𝑛 2 + 𝑐) = 𝜋ℓ1 (𝑎) + 𝜋ℓ2 (𝑐)𝑛 1 ; combining with (51) yields 𝑌 phys [𝑟, 𝑞] = 𝑌 [𝜋ℓ (𝑟 ), 𝑞]. □
A.8
Input-Layout Local FFT Correctness
Theorem 11 (Input-layout local FFT identity). For every ℓ ≥ 0 the input-layout traversal implements Bℓ = Ω𝑏 ℓ 𝑃𝜋ℓ . In particular, if 𝑥 phys = 𝑃𝜋ℓ 𝑥, then Bℓ 𝑥 phys = Ω𝑏 ℓ 𝑥.
B.1
Scheduler Pseudocode
Figure 2 realizes the scheduler described in the main text: it invokes a constant-size kernel at the leaves, performs a batch split when 𝑠 ≥ 𝑛 (splitting an outer-batch digit whenever one is available, and the distinguished contiguous field only otherwise), and performs a balanced digit-aligned transform split when 𝑠 < 𝑛. The transform split in Figure 2 is only a recursive view transformation: by the no-physical-transpose property of Appendix A, the first recursive call computes the first-stage transforms, the physical 𝑇1 -slots are then relabeled as the output field 𝐴, and the second call treats 𝐴 as an outer-batch field. The layout-aware twiddle is evaluated from the logical coordinates of the current view and fused into the second child group, so the scheduler introduces no separate recursive twiddle pass or physical transpose.
Shina Guo, Weiguo Gao, and Yuan Tang
B.2
Local Work Bound
Proposition 7 (Local work). For 𝑠 transforms of length 𝑛 = 𝑏 ℓ , the scheduler performs 𝑊local (𝑛, 𝑠) = 𝑂 𝑛𝑠 (1 + log 𝑛) (53) local operations. Proof. We separate the scheduler work into three parts: (i) constant-size leaf DFT work, (ii) fused twiddle and transform-split control work, and (iii) batch-split control work. Fix one of the 𝑠 logical transforms at the root. Ignore batch splits momentarily: they only partition the collection of independent transforms and never alter the Cooley–Tukey recursion followed by an individual transform. Whenever a transform subproblem of length 𝑛 ′ = 𝑏𝑑 is split into factors 𝑛 1 = 𝑏𝑑1 and 𝑛 2 = 𝑏𝑑2 , where 𝑑 1 + 𝑑 2 = 𝑑, the first child evaluates the 𝑑 1 radix-𝑏 FFT levels associated with the first factor and the second child evaluates the 𝑑 2 levels associated with the second factor. Thus, over the complete transform recursion, the transform dependencies of the fixed logical transform contain exactly log𝑏 𝑛 radix-𝑏 levels, up to a constant factor depending only on the constant-size leaf convention. At each radix-𝑏 level, every item participates in 𝑂 (1) fixed-radix butterfly and twiddle operations. The twiddle of a transform split is fused into the first access of the corresponding second child group, but it is still charged to that same FFT level. Hence one logical transform incurs 𝑂 (𝑛 log𝑏 𝑛) leaf, butterfly, and twiddle work. Summed over all 𝑠 transforms, 𝑊FFT (𝑛, 𝑠) = 𝑂 (𝑛𝑠 log𝑏 𝑛).
(54)
It remains to bound scheduler overhead. A batch split performs only 𝑂 (1) descriptor and loop-control work per child. Its 𝑏 children partition the parent batch instances, and 𝑏 = 𝑂 (1). Therefore, over all batch-only recursion nodes associated with any fixed transform subproblem, the total administrative work is linear in the number of batch instances represented by that subproblem. Summing over the transform recursion gives 𝑂 (𝑛𝑠) additional work. The constant-size leaf calls also contribute 𝑂 (𝑛𝑠) work. Consequently, 𝑊local (𝑛, 𝑠) = 𝑂 (𝑛𝑠 log𝑏 𝑛) + 𝑂 (𝑛𝑠) = 𝑂 𝑛𝑠 (1 + log 𝑛) . Finally, Appendix A proves that transform splits are view transformations rather than physical transposes; thus no additional Θ(𝑛𝑠) data-rearrangement pass is incurred at each recursive transform split. □
B.3
Recursion Depth and Stack Space
Lemma 13 (Recursion depth). Along every root-to-leaf path, the scheduler performs 𝑂 1 + log log𝑏 𝑛 (55) transform splits and 𝑂 log𝑏 (𝑛𝑠) (56) batch splits. Hence its total recursion depth is 𝑂 (log(𝑛𝑠)). Proof. Let 𝑑 = log𝑏 𝑛 be the number of active transform digits. A balanced transform split leaves either ⌈𝑑/2⌉ or ⌊𝑑/2⌋ transform digits in a child. Repeated balanced splitting therefore reaches a constant-size transform after 𝑂 (1 + log 𝑑) transform splits, which gives (55).
Every batch split fixes one active batch digit. Initially there are log𝑏 𝑠 batch digits. A transform split may move transform digits into an outer-batch field, but along one root-to-leaf path the total number of such transferred digits is at most log𝑏 𝑛. Thus no more than log𝑏 (𝑛𝑠) active batch digits can be fixed, proving (56). Summing the two bounds gives the depth bound. □ Corollary 14 (Stack space). The scheduler uses 𝑆 stack = 𝑂 (log(𝑛𝑠)) words of recursion-stack space. Proof. Each stack frame stores only a constant number of field descriptors, fixed-digit assignments, loop counters, and mode bits. The result follows from Lemma 13. □
B.4
Analytical Cache-Fitting Frontier
We now analyze cache complexity. The following frontier is used only in the proof; the scheduler never tests 𝑀, 𝐵, or the cache-fitting condition. Fix a sufficiently small constant 𝛼 > 0, and define 𝑀 ★ = 𝛼𝑀. A view is cache fitting if its mass satisfies 𝑛𝑠 ≤ 𝑀 ★ .
(57)
A maximal cache-fitting view is either the root, when the root is cache fitting, or a cache-fitting view whose parent is not cache fitting. Its purpose is to identify subproblems whose complete descendant execution can remain resident in the ideal cache; the scheduler continues recursively below the frontier, and the frontier merely groups that descendant work for analysis. We record the mass behavior of the two split types, used repeatedly below. A batch split partitions the batched transforms into 𝑏 data-disjoint children, so the child masses sum to the parent mass 𝑛𝑠. A transform split 𝑛 = 𝑛 1𝑛 2 produces child groups with parameters (𝑛 1, 𝑠𝑛 2 ) and (𝑛 2, 𝑠𝑛 1 ), each of mass 𝑛 1 (𝑠𝑛 2 ) = 𝑛 2 (𝑠𝑛 1 ) = 𝑛𝑠.
(58)
Thus a transform split preserves the mass of each child group.
B.5
Frontier Geometry
Lemma 15 (Frontier entry is through a batch split). A noncache-fitting view can first enter the cache-fitting region only through a batch split. Proof. A transform split preserves the mass of each child group, by (58), and therefore cannot turn a non-cache-fitting mass 𝑛𝑠 > 𝑀 ★ into a cache-fitting one. A batch split reduces the batch size of each child by a factor of 𝑏 and can therefore produce a cache-fitting child. □ Lemma 16 (Mass of a nonroot frontier view). Every nonroot maximal cache-fitting view satisfies 𝑀★ < 𝑛𝑠 ≤ 𝑀 ★ . (59) 𝑏 Proof. By Lemma 15, a nonroot maximal cache-fitting view is created by a batch split. Its parent has the same transform length and 𝑏 times as many batch instances. Since the child fits but the parent does not, 𝑛𝑠 ≤ 𝑀 ★ and 𝑏 𝑛𝑠 > 𝑀 ★, which is equivalent to (59). □
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
B.6
Contiguous Width at the Frontier
Hence, after 𝐾 transform splits,
Lemma 17 (Contiguous width at the frontier). Suppose the root view has initial contiguous width 𝑐 0 ≥ 𝐵. Then every maximal cache-fitting view has 𝑐 = Ω(𝐵). (60) Proof. If the distinguished field 𝐶 is never split on the root-tofrontier path, then its width remains 𝑐 = 𝑐 0 ≥ 𝐵. Otherwise, consider the last split of 𝐶 before the frontier view, and let 𝑐 denote the child width immediately after that split. The parent width was 𝑏𝑐. Since the scheduler splits 𝐶 only when 𝑂 = ∅, the parent batch size was also 𝑏𝑐, and the batch-split condition gives 𝑛 ≤ 𝑏𝑐 at that split. If the resulting child is already cache fitting, its parent is not, so 𝑀 ★ < 𝑛 𝑏𝑐 ≤ 𝑏 2𝑐 2 ; if the child is not yet cache fitting, then 𝑀 ★ < 𝑛𝑐 ≤ 𝑏𝑐 2 . In either case √ √ 𝑐 = Ω( 𝑀 ★) = Ω( 𝑀), (61) a bound stronger than (60); only (60) is used downstream. Subsequent outer-batch and transform splits do not change 𝑐, and the tall-cache assumption 𝑀 = Ω(𝐵 2 ) turns (61) into (60). □
B.7
Transform Length and Number of Transform Splits
Lemma 18 (Transform–batch relation). There is a constant 𝜅𝑏 = 𝑂 (1), depending only on the fixed radix 𝑏, such that every view reached after at least one transform split satisfies 𝑠 ≤ 𝜅𝑏 𝑛 3 .
(62)
Proof. Suppose a transform split is applied to a view with 𝑛 = 𝑛 1𝑛 2 and 𝑠 < 𝑛. The balanced digit-aligned split satisfies 1 𝑏 ≤ 𝑛 1 /𝑛 2 ≤ 𝑏. For the first child the transform length is 𝑛 1 and the batch size is 𝑠𝑛 2 ; since 𝑠 < 𝑛 1𝑛 2 , 𝑠𝑛 2 < 𝑛 1𝑛 22 ≤ 𝑏 2𝑛 31 . The second child is symmetric, 𝑠𝑛 1 < 𝑛 2𝑛 21 ≤ 𝑏 2𝑛 32 . Thus both immediate children satisfy (62) with 𝜅𝑏 = 𝑏 2 . A later batch split decreases 𝑠 without changing 𝑛, while a later transform split reestablishes the same inequality for its children. The claim follows by induction along the recursion path. □ Lemma 19 (Transform length at the frontier). Let a nonroot maximal cache-fitting view be reached along a path containing at least one transform split. Then 𝑛 = Ω(𝑀 1/4 ).
(63)
Proof. By Lemma 16, 𝑀 ★/𝑏 < 𝑛𝑠, and by Lemma 18, 𝑠 ≤ 𝜅𝑏 𝑛 3 . Hence 𝑀 ★/𝑏 < 𝜅𝑏 𝑛 4 , which proves (63). □ We index transform splits by their output length: let 𝑛 0 = 𝑛 be the root transform length, and let 𝑛 𝑗 denote the transform length of the active child group immediately after the 𝑗-th transform split on a root-to-frontier path. Lemma 20 (Transform-size shrinkage). For every 𝑗 ≥ 1, √︁ 𝑛 𝑗 ≤ 𝑏 𝑛 𝑗 −1 . (64)
𝐾
𝐾
𝐾
𝑛𝐾 ≤ 𝑏 1−1/2 𝑛 1/2 < 𝑏 𝑛 1/2 .
(65) = 𝑏𝑑 , the larger
Proof. If the length before the 𝑗-th split is 𝑛 𝑗 −1 child contains at most ⌈𝑑/2⌉ transform digits, so 𝑛 𝑗 ≤ 𝑏 ⌈𝑑/2⌉ ≤ 1−1/2 𝑗 𝑛 1/2 𝑗 , 𝑏 1/2𝑛 1/2 𝑗 −1 , which is (64). Iterating from 𝑛 0 = 𝑛 gives 𝑛 𝑗 ≤ 𝑏 and in particular (65). □ Lemma 21 (Number of transform splits before the frontier). If 𝐾 is the number of transform splits on a root-to-frontier path, then 2𝐾 = 𝑂 1 + log𝑀 𝑛 . (66) Proof. The claim is immediate when 𝐾 = 0, so assume 𝐾 ≥ 1. By Lemma 19 the frontier transform length satisfies 𝑛𝐾 ≥ 𝑐𝑏 𝑀 1/4 for a constant 𝑐𝑏 > 0 depending only on 𝑏 and 𝛼. Combining this with (65), 𝐾 𝑐𝑏 𝑀 1/4 ≤ 𝑛𝐾 < 𝑏 𝑛 1/2 . Taking base-2 logarithms, 1 log2 𝑛 > 14 log2 𝑀 − 𝑂 (1). 2𝐾 For 𝑀 larger than a suitable constant the right-hand side is at least 1 8 log2 𝑀, so log 𝑛 2𝐾 < 1 2 = 8 log𝑀 𝑛, 8 log2 𝑀 i.e. 2𝐾 = 𝑂 (log𝑀 𝑛). In the remaining regime 𝑛 = 𝑂 (𝑀), one has 1 + log𝑀 𝑛 = 𝑂 (1), and only 𝑂 (1) balanced transform splits can occur before the transform length becomes constant, so 2𝐾 = 𝑂 (1). The two cases give (66). □
B.8
Block Cover of a Frontier View
Lemma 22 (Block cover of one frontier view). A maximal cache-fitting view of mass 𝑛𝑠 intersects 𝑛𝑠 𝑂 1+ (67) 𝐵 cache blocks. Its complete depth-first descendant execution incurs the same asymptotic number of cache misses. Proof. By the active-run decomposition of Appendix A, the view is the disjoint union of 𝑛𝑠/𝑐 contiguous runs of length 𝑐. Therefore its block footprint is at most 𝑛𝑠 𝑛𝑠 𝑛𝑠 𝑛𝑠 𝑐 𝑛𝑠 1+ =𝑂 + =𝑂 =𝑂 1+ , 𝑂 𝑐 𝐵 𝑐 𝐵 𝐵 𝐵 where the third equality uses 𝑐 = Ω(𝐵) from Lemma 17, which gives 𝑛𝑠/𝑐 = 𝑂 (𝑛𝑠/𝐵). Every descendant of the frontier view accesses only items in the frontier view’s working set: batch splits select subsets of its runs, contiguous-field splits subdivide those runs, and transform splits reinterpret the same addresses through sequential child groups. Moreover the view contains at most 𝑀 ★ = 𝛼𝑀 items; choosing 𝛼 sufficiently small ensures that its block footprint, together with constant recursion metadata, fits in cache. Under depth-first execution these blocks can remain resident while the complete descendant computation is performed. Hence (67) also bounds the descendant cache cost. □
Shina Guo, Weiguo Gao, and Yuan Tang
B.9
Total Frontier Mass
Lemma 23 (Total frontier mass). Suppose every root-to-frontier path contains at most 𝐾 transform splits. Then the sum of the masses of all maximal cache-fitting views is at most 𝑛𝑠 2𝐾 . Proof. We prove the more general statement that a recursion subtree rooted at a view of mass 𝜇, with at most 𝑘 remaining transform splits on each path, has total frontier mass at most 𝜇 2𝑘 . If the root of the subtree is cache fitting, it is itself a frontier view and contributes mass 𝜇 ≤ 𝜇 2𝑘 . If the root performs a batch split, the child masses sum to 𝜇; applying the induction hypothesis to every child gives total frontier Í mass at most 2𝑘 𝑢 mass(𝑢) = 𝜇 2𝑘 . If the root performs a transform split, both child groups have mass 𝜇 by (58), and each has at most 𝑘 −1 remaining transform splits. Their combined frontier mass is therefore at most 2𝜇 2𝑘 −1 = 𝜇 2𝑘 . Applying the claim at the original root proves the lemma. □
B.10
Cache Complexity
Theorem 24 (Cache complexity of the batched local FFT). Let an admissible root view contain 𝑠 transforms of length 𝑛, with initial contiguous width 𝑐 ≥ 𝐵. Under the tall-cache assumption 𝑀 = Ω(𝐵 2 ), the scheduler incurs 𝑛𝑠 1 + log𝑀 𝑛 𝑄 local (𝑛, 𝑠) = 𝑂 (68) 𝐵 cache misses. Proof. Let 𝐾 be the maximum number of transform splits on a root-to-frontier path. Lemma 21 gives 2𝐾 = 𝑂 (1 + log𝑀 𝑛), so by Lemma 23 the total mass of the maximal cache-fitting frontier is 𝑂 𝑛𝑠 (1 + log𝑀 𝑛) . (69) By Lemma 22, a frontier view of mass 𝜇 costs 𝑂 (1 + 𝜇/𝐵) cache misses. The additive constant is absorbed into the mass term: by Lemma 17 every frontier view has contiguous width Ω(𝐵), hence mass 𝜇 = Ω(𝐵). Summing the cache costs over the frontier and using (69) gives (68). □
B.11 Application to the Two PACO Local Phases In Phase I, each processor performs 𝑁 2 /𝑝 transforms of length 𝑁 1 ; its initial contiguous batch field is precisely the local input slab, so 𝑁2 𝑛 = 𝑁 1, 𝑠 =𝑐 = . (70) 𝑝 The slab-width assumption 𝑁 2 /𝑝 ≥ 𝐵 permits direct application of Theorem 24, yielding 𝑁 𝑄𝑞,I = 𝑂 (71) 1 + log𝑀 𝑁 1 . 𝑝𝐵 Similarly, in Phase III each processor performs 𝑁 1 /𝑝 transforms of length 𝑁 2 , with 𝑁1 𝑛 = 𝑁 2, 𝑠 =𝑐 = . (72) 𝑝 Since 𝑁 1 /𝑝 ≥ 𝐵, 𝑁 𝑄𝑞,III = 𝑂 1 + log𝑀 𝑁 2 . (73) 𝑝𝐵
Finally, log𝑀 𝑁 1 + log𝑀 𝑁 2 = log𝑀 𝑁 . Adding (71) and (73), and absorbing the linear cache cost of the fused middle stage, gives 𝑁 𝑄𝑞,I + 𝑄𝑞,III = 𝑂 1 + log𝑀 𝑁 . (74) 𝑝𝐵
C
Data Movement: Local Primitives and the Fused Redistribution
This appendix specifies the data-movement machinery used in Phases I–III. It first defines the local primitives—macro-element digit reversal, its destination-indexed fusion, and cache-oblivious rectangular transposition—and then composes them into the explicit per-processor procedure that realizes PACO’s fused middle redistribution. The relevant permutations operate on contiguous macro-elements, not on individual scalars. The cache line size is 𝐵 entries. In both local applications every macro-element spans at least one line: by the slab-width conditions (5) the two uses below have macro-element length 𝑁 2 /𝑝 and 𝑁 1 /𝑝, respectively, and both are ≥ 𝐵. Sections 3.3 and 4.2 already established the fused middle stage at the level needed for the main result: its coordinate trace (18), its endto-end correctness (Theorem 3), its balanced migration counts, and its per-processor local cost (Theorem 4). The second half of this appendix supplies the single deferred component, namely the explicit cache-oblivious procedure that realizes the mapping, together with (i) a line-by-line correctness argument and (ii) the workspace and transpose-aggregation accounting promised in Sect. 4.2. It reuses the source ownership (1), the factor-swapped target ownership (3), e 𝑐] = 𝐴[𝜋𝑚1 (𝑎), 𝑐], the digit-reversal notathe Phase-I identity 𝐴[𝑎, tion of Appendix A.1, and the slab-width conditions (5).
C.1
Macro-Element Digit Reversal
Let the array hold 𝑛 = 𝑏 ℓ contiguous macro-elements, each of length 𝑠; write ℓ = log𝑏 𝑛 (the count of macro-elements, distinct from the global 𝑚 = log𝑏 𝑁 ). For 𝑖 ∈ [𝑛] with base-𝑏 digits 𝑑 ℓ −1 . . . 𝑑 0 , the digit-reversal permutation is 𝜋ℓ (𝑖) = 𝑑 0𝑏 ℓ −1 + · · · + 𝑑 ℓ −1 , an involution, so every orbit is a fixed point or a two-cycle. For the array 𝑋 [0 . . 𝑛 − 1] [0 . . 𝑠 − 1], the desired rearrangement is 𝑋 ′ [𝑖, 𝑗] = 𝑋 [𝜋ℓ (𝑖), 𝑗].
(75)
Since 𝜋ℓ is an involution, (75) is realized in place by exchanging the two macro-elements of each nontrivial two-cycle; the condition 𝑖 < 𝜋ℓ (𝑖) selects one representative per two-cycle. Macro-Digit-Reversal(𝑋, 𝑛, 𝑠) 1 2 3 4 5
// ℓ = log𝑏 𝑛; reverse the 𝑛 contiguous macro-elements for 𝑖 = 0 to 𝑛 − 1 if 𝑖 < 𝜋ℓ (𝑖) for 𝑗 = 0 to 𝑠 − 1 𝑋 [𝑖, 𝑗] ↔ 𝑋 [𝜋ℓ (𝑖), 𝑗]
Lemma 25 (Correctness). The procedure Macro-Digit-Reversal implements (75). Proof. If 𝑖 = 𝜋ℓ (𝑖), the macro-element is a fixed point and left unchanged, agreeing with (75). Otherwise 𝑖 lies in the twocycle {𝑖, 𝜋ℓ (𝑖)}, and exactly one member satisfies 𝑖 < 𝜋ℓ (𝑖); when processed, the two macro-elements are exchanged entry by entry, so afterward 𝑋 [𝑖, 𝑗] = 𝑋 old [𝜋ℓ (𝑖), 𝑗] for both members. Every nonfixed macro-element is exchanged exactly once. □
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
Lemma 26 (Cost). If 𝑠 ≥ 𝐵, the procedure Macro-Digit-Reversal has local work 𝑂 (𝑛𝑠) and cache complexity 𝑂 (𝑛𝑠/𝐵), using 𝑂 (1) auxiliary scalar storage. Proof. Each macro-element is scanned at most once, giving 𝑂 (𝑛𝑠) scalar reads and writes. Each contiguous segment of 𝑠 entries costs 𝑂 (𝑠/𝐵) misses, and neither member of a two-cycle is revisited after its exchange; summing over 𝑛 macro-elements gives 𝑂 (𝑛𝑠/𝐵). One temporary scalar suffices. □
C.2
Fusing Destination-Indexed Unary Operations
A macro-element rearrangement may be fused with a unary operation whose coefficient depends on the destination. Let 𝑔(𝑖, 𝑗) be the coefficient applied to the entry written to (𝑖, 𝑗); the fused transformation is 𝑋 ′ [𝑖, 𝑗] = 𝑔(𝑖, 𝑗) 𝑋 [𝜋ℓ (𝑖), 𝑗].
(76)
For a nontrivial two-cycle {𝑖, 𝑘 }, 𝑘 = 𝜋ℓ (𝑖), 𝑖 < 𝑘, the writes are 𝑋 [𝑖, 𝑗] = 𝑔(𝑖, 𝑗)𝑋 old [𝑘, 𝑗] and 𝑋 [𝑘, 𝑗] = 𝑔(𝑘, 𝑗)𝑋 old [𝑖, 𝑗]; for a fixed point, 𝑋 [𝑖, 𝑗] = 𝑔(𝑖, 𝑗)𝑋 [𝑖, 𝑗]. Fused-Macro-Digit-Reversal(𝑋, 𝑛, 𝑠, 𝑔) 1 for 𝑖 = 0 to 𝑛 − 1 2 if 𝑖 < 𝜋ℓ (𝑖) 3 for 𝑗 = 0 to 𝑠 − 1 4 𝑥 = 𝑋 [𝑖, 𝑗] 5 𝑦 = 𝑋 [𝜋ℓ (𝑖), 𝑗] 6 𝑋 [𝑖, 𝑗] = 𝑔(𝑖, 𝑗) · 𝑦 7 𝑋 [𝜋ℓ (𝑖), 𝑗] = 𝑔(𝜋ℓ (𝑖), 𝑗) · 𝑥 8 elseif 𝑖 = = 𝜋ℓ (𝑖) 9 for 𝑗 = 0 to 𝑠 − 1 10 𝑋 [𝑖, 𝑗] = 𝑔(𝑖, 𝑗) · 𝑋 [𝑖, 𝑗] Lemma 27 (Correctness). The procedure Fused-Macro-DigitReversal implements (76). Proof. For a nontrivial two-cycle the two writes place each old macro-element at its digit-reversed destination and scale by the destination-indexed coefficient; for a fixed point the update is exactly (76). As in Lemma 25, every two-cycle is processed once. □ Corollary 28 (Cost of fused rearrangement). If 𝑠 ≥ 𝐵 and evaluating 𝑔 takes 𝑂 (1) work, the procedure Fused-Macro-DigitReversal has local work 𝑂 (𝑛𝑠) and cache complexity 𝑂 (𝑛𝑠/𝐵). Proof. The fused operation adds 𝑂 (1) arithmetic per entry to the scans of Lemma 26, with neither an extra traversal nor extra storage. □
Application to twiddle multiplication. In the middle stage the destination-indexed operation is twiddle multiplication: with destination macro-element index 𝑣 and within-element (local) coordinate 𝑐 loc , 𝑔(𝑣, 𝑐) = 𝜔 𝑁−𝑣𝑐 , where 𝑐 = 𝑞 · (𝑁 2 /𝑝) + 𝑐 loc is the global column. Hence twiddles are applied when an entry is written to its destination; no separate pass is required.
C.3
Use in Phases I and III
Both local reversals are covered by Lemma 26 and Corollary 28 with a single parameterized statement. In Phase I and III each processor reverses 𝑛 macro-elements of length 𝑠 ∈ {𝑁 2 /𝑝, 𝑁 1 /𝑝}; by condition (5), 𝑠 ≥ 𝐵, so rearranging the 𝑛 macro-elements costs 𝑂 (𝑛𝑠) work and 𝑂 (𝑛𝑠/𝐵) misses. Instantiating 𝑠 = 𝑁 2 /𝑝 (Phase I) and 𝑠 = 𝑁 1 /𝑝 (Phase III), with 𝑛 = 𝑁 1 and 𝑛 = 𝑁 2 respectively, gives 𝑂 (𝑁 1 𝑁 2 /𝑝) = 𝑂 (𝑁 /𝑝) work and 𝑂 (𝑁 /(𝑝𝐵)) misses per reversal. The macro-element length lower bound 𝑠 ≥ 𝐵 is exactly what ensures spatial locality, so no cache-oblivious traversal of a scalar digit-reversal permutation is needed.
C.4
Local Rectangular Transposition
After communication, each received block is a rectangle of dimensions 𝑁1 𝑁2 𝑟 ×𝑐 = × . (77) 𝑝 𝑝 PACO transposes it with the cache-oblivious rectangular-transpose algorithm of Frigo et al. [8], whose standard bound is: transposing an 𝑟 × 𝑐 rectangle takes 𝑂 (𝑟𝑐) work and 𝑂 (1 + 𝑟𝑐/𝐵) misses. For a block of dimensions (77) this is 𝑂 𝑁𝑝1 𝑁2 2 = 𝑂 𝑝𝑁2 work, 𝑂 1 + 𝑝𝑁2 𝐵 misses, (78) and the reusable transpose workspace is (𝑁 1 /𝑝)(𝑁 2 /𝑝) = 𝑁 /𝑝 2 entries. This subsection is the single source for the per-block transpose accounting; its aggregation over the 𝑝 blocks of a processor is carried out in Appendix C.7 (and used in Theorem 4).
C.5
The Per-Processor Fused Redistribution Procedure
Processor 𝑞 enters the middle stage owning an 𝑁 1 × (𝑁 2 /𝑝) Phase-I slab 𝑋𝑞 . The procedure Para-Co-Drp𝑏 composes the three primitives above: it packs the slab by reversing the 𝑁 1 -field (Appendix C.1) while fusing the top-level twiddle (Appendix C.2), exchanges one balanced set of blocks in a single BSP ℎ-relation, transposes each received block cache-obliviously (Appendix C.4), and unpacks by reversing the 𝑁 2 -field. Each field reversal is an in-place sequence of contiguous macro-element swaps; the fused twiddle multiplication is folded into the packing scan, costing no separate pass.
Shina Guo, Weiguo Gao, and Yuan Tang
Para-Co-Drp𝑏 (𝑋𝑞 , 𝑌𝑞 , 𝑁 1, 𝑁 2, 𝑝, 𝑞)
1
2 3 4 5 6 7 8 9 10 11 12 13 14
// 𝑋𝑞 is the local input array, 𝑌𝑞 is the output // 𝑞 is the local processor id // Local pack: 𝑁 1 -field reversal with the top-level twiddle view 𝑋𝑞 as 𝑁 1 row macro-elements of length 𝑁 2 /𝑝 // 𝑔(𝑣, 𝑐 loc ) = 𝜔 𝑁−𝑣𝑐 , where 𝑐 = 𝑞 · (𝑁 2 /𝑝) + 𝑐 loc (global column) , 𝑁 = 𝑁 1𝑁 2 Fused-Macro-Digit-Reversal(𝑋𝑞 , 𝑁 1, 𝑁 2 /𝑝, 𝑔) // One balanced BSP all-to-all exchange partition the packed slab into blocks 𝑆𝑞,𝑞 ′ of shape (𝑁 1 /𝑝) × (𝑁 2 /𝑝), 𝑞 ′ ∈ [𝑝] place the self block 𝑆𝑞,𝑞 in its raw output region 𝑅𝑞,𝑞 for 𝑞 ′ = 0 to 𝑝 − 1 if 𝑞 ′ ≠ 𝑞 send 𝑆𝑞,𝑞 ′ to raw region 𝑅𝑞 ′ ,𝑞 on processor 𝑞 ′ // Local rectangular transpose of each received block allocate or reuse one scratch array 𝑇 of 𝑁 /𝑝 2 items for 𝑟 = 0 to 𝑝 − 1 CO-Transpose(𝑅𝑞,𝑟 ,𝑇 , 𝑁 1 /𝑝, 𝑁 2 /𝑝) copy 𝑇 to the output region 𝑌𝑞,𝑟 // Local unpack: 𝑁 2 -field digit reversal view the local array as 𝑁 2 row macro-elements of length 𝑁 1 /𝑝 Macro-Digit-Reversal(𝑌𝑞 , 𝑁 2, 𝑁 1 /𝑝) return the local 𝑁 2 × (𝑁 1 /𝑝) target slab
The 𝑝 raw incoming regions 𝑅𝑞,· together form the mandatory target slab; the only auxiliary storage is the single transpose scratch 𝑇 of 𝑁 /𝑝 2 items.
C.6
Line-by-Line Correctness
C.7
Cost and Workspace
Each primitive was costed in Appendices C.1–C.4. The two fused field reversals (packing 𝑁 1 macro-elements of length 𝑁 2 /𝑝; unpacking 𝑁 2 of length 𝑁 1 /𝑝) each cost 𝑂 (𝑁 /𝑝) work and 𝑂 (𝑁 /(𝑝𝐵)) misses by Corollary 27, since both macro-element widths are ≥ 𝐵 by (5); the fused twiddle adds Θ(𝑁 /𝑝) arithmetic with no extra cache pass. The 𝑝 block transposes cost 𝑂 (𝑁 /𝑝) work and, by (74) together with 𝑝 ≤ 𝑁 /(𝑝𝐵 2 ) = 𝑂 (𝑁 /(𝑝𝐵)), aggregate to 𝑂 (𝑁 /(𝑝𝐵)) misses. Summing reproduces (22)–(24). The only nonoutput workspace is one reusable 𝑁 /𝑝 2 -item transpose buffer plus 𝑂 (log 𝑁 ) recursion-stack words, giving (28); the balanced counts 𝜇 = 𝑁 (1 − 1/𝑝) and ℎ = (𝑁 /𝑝)(1 − 1/𝑝) are as in Theorem 4.
D
Lower Bounds: Deferred Proofs
The main text establishes three of PACO’s four lower bounds in full: the one-round redistribution bound (Proposition 3), the exact migration bound for the prescribed fused permutation (Proposition 4), and the maximum per-processor work bound (Proposition 5). The maximum per-processor cache lower bound (Proposition 6) was stated there with only a proof roadmap; its detailed segment argument is deferred to this appendix. We import the FFT-DAG input-boundary theorem [15], convert boundary growth into block-cache misses, obtain a finite blockaware bound for a single radix-2 FFT DAG, extend it to disjoint batches, and specialize it to PACO’s two local phases, completing Proposition 6. As a byproduct we also record a stronger constantone asymptotic miss bound. Throughout, arithmetic vertices of the fixed-radix Cooley–Tukey DAG are never recomputed.
Remark 1 (Fixed-radix scope). Let 𝑏 = 2𝑡 , 𝑡 = 𝑂 (1). If every radix-𝑏 butterfly is expanded into its 𝑡 constituent radix-2 layers, the resulting dependency graph is the standard radix-2 FFT DAG. Hence all lower bounds of this appendix—the single-transform bounds of Appendices D.3 and D.4, the batched bounds of Appendix D.6, and the Lemma 29 (The procedure implements 𝜋𝑚 ). For an item stored at Phase-I physical coordinate (𝑎, 𝑐) on its source processor 𝑞, write 𝑐 = PACO local bound of Appendix D.7—apply to the expanded fixed-radix 𝑞(𝑁 2 /𝑝)+𝑐 loc . Then Para-Co-Drp𝑏 places it at (𝑢, 𝑣) = 𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎) , DAG. If a complete radix-𝑏 DFT is instead treated as one indivisible atomic vertex, the graph structure and leading constants differ; that i.e. 𝜋𝑚 (𝑎 ++ 𝑐) = 𝜋𝑚2 (𝑐) ++ 𝜋𝑚1 (𝑎), owned by ⌊𝜋𝑚1 (𝑎)/(𝑁 1 /𝑝)⌋. atomic model is outside the scope of these claims. Proof. The row reversal moves the item to packed row 𝑣 = 𝜋𝑚1 (𝑎) with column 𝑐 unchanged; the fused scan forms 𝑍 [𝑣, 𝑐] = D.1 Imported FFT-DAG Input-Boundary 𝐴[𝑣, 𝑐] 𝜔 𝑁−𝑣𝑐 (the local-pack step of Sect. 4.2). Decomposing 𝑣 = Theorem 𝑞 ′ (𝑁 1 /𝑝) + 𝑣 loc identifies the unique destination 𝑞 ′ = ⌊𝑣/(𝑁 1 /𝑝)⌋, Let 𝐺 2 (𝑁 ) denote the radix-2 FFT DAG on 𝑁 = 2𝑚 inputs; it has which is the owner under (3), so the item is assigned to block 𝑆𝑞,𝑞 ′ 𝑚 + 1 levels of 𝑁 vertices and 𝑁 log2 𝑁 noninput computation at block-local coordinate (𝑣 loc, 𝑐 loc ) and routed to region 𝑅𝑞 ′ ,𝑞 on vertices. For a set 𝑈 of noninput vertices, its input boundary is ′ 𝑞 . Because raw regions are ordered by source processor, the pair in(𝑈 ) = {𝑣 ∉ 𝑈 : ∃ 𝑢 ∈ 𝑈 with (𝑣, 𝑢) ∈ 𝐸}. (79) (𝑞, 𝑐 loc ) recovers 𝑐 = 𝑞(𝑁 2 /𝑝) + 𝑐 loc . The block transpose sends (𝑣 loc, 𝑐 loc ) ↦→ (𝑐 loc, 𝑣 loc ), and the 𝑁 2 -field reversal maps the reconTheorem 30 (FFT input-boundary growth). For every intestructed source row 𝑐 to 𝑢 = 𝜋𝑚2 (𝑐), leaving 𝑣 = 𝜋𝑚1 (𝑎) intact. The ger 𝑘 ≥ 1, every set 𝑈 of noninput vertices of 𝐺 2 (𝑁 ) with |𝑈 | ≥ result (𝑢, 𝑣) = (𝜋𝑚2 (𝑐), 𝜋𝑚1 (𝑎)) is precisely the trace (18); its owner 𝑘 log2 (4𝑘) satisfies |in(𝑈 )| ≥ 𝑘. is ⌊𝑣/(𝑁 1 /𝑝)⌋ by (3). □ Theorem 30 is imported from the boundary-growth literature [15]; By Lemma 29 the procedure Para-Co-Drp𝑏 produces, for each we use only its graph-theoretic consequence. fixed 𝑣 = 𝑘 1 , the local layout 𝑍b[𝑢, 𝑣] = 𝑍 [𝑣, 𝜋𝑚2 (𝑢)] of (17): exactly the digit-reversed input layout consumed by B𝑚2 = Ω𝑁2 𝑃𝜋𝑚2 . D.2 From Input Boundaries to Block Cache Hence the composition (10) of Theorem 3 applies verbatim, and Misses PACO returns the canonical coefficients 𝑌 [𝑘 2, 𝑘 1 ] = 𝑦 [𝑘 1 + 𝑘 2 𝑁 1 ] Lemma 31 (Input boundary to block-cache misses). Consider under the target ownership (3). one contiguous segment of a legal schedule, and let 𝑈 be the set of
We verify that the procedure Para-Co-Drp𝑏 realizes exactly the coordinate trace (18); the algebraic correctness of the whole pipeline then follows from Theorem 3.
PACO: A Fully Cache-Oblivious Parallel FFT with One Global Redistribution
previously uncomputed FFT vertices evaluated during it. If |in(𝑈 )| ≥ 𝑘, the segment incurs at least (𝑘 − 𝑀)/𝐵 + cache misses, where ⌈𝑥⌉ + := max{0, ⌈𝑥⌉}.
D.5
Proof. At the start of the segment at most 𝑀 items reside in cache, so at most 𝑀 of the required boundary values are resident; at least 𝑘 − 𝑀 are absent when 𝑘 > 𝑀. Each miss loads one block of at most 𝐵 items and thus supplies at most 𝐵 missing values, so at least ⌈(𝑘 − 𝑀)/𝐵⌉ misses are required when 𝑘 > 𝑀, and zero otherwise. □
Lemma 34 (Boundary growth for disjoint batches). For every integer 𝑘 ≥ 1, every set 𝑈 of noninput vertices of 𝐺 2(𝑠 ) (𝑛) with |𝑈 | ≥ 4𝑘 log2 (4𝑘) satisfies |in(𝑈 )| ≥ 𝑘.
Lemma 31 permits a strong algorithmic model: the algorithm may know 𝑀 and 𝐵, choose arbitrary slow-memory layouts, pack 𝐵 useful values per loaded block, keep the most useful 𝑀 values in cache, and maintain multiple slow-memory copies. None of these eliminates the need to load a block when a required boundary value is absent.
D.3
Finite Block-Aware Lower Bound for One FFT DAG
Theorem 32 (Finite block-aware cache lower bound). For every integer 𝑘 > 𝑀, every no-recomputation topological schedule of 𝐺 2 (𝑁 ) satisfies kl m j 𝑁 log2 𝑁 𝑘 −𝑀 𝑄 miss (𝑁 , 𝑀, 𝐵) ≥ 𝑘 log (4𝑘 . (80) 𝐵 ) 2
Proof. Partition the schedule into consecutive complete segments, each computing exactly 𝑘 log2 (4𝑘) previously uncomputed vertices (ignore a final incomplete segment). Since 𝐺 2 (𝑁 ) has 𝑁 log2 𝑁 noninput vertices, there are at least ⌊𝑁 log2 𝑁 /(𝑘 log2 (4𝑘))⌋ complete segments. Each has |in(𝑈 )| ≥ 𝑘 by Theorem 30 and thus at least ⌈(𝑘 − 𝑀)/𝐵⌉ misses by Lemma 31; disjoint segments add. □
D.4
A Stronger Constant-One Asymptotic Bound
This subsection is not needed for PACO’s upper-bound proof; it records the stronger miss-only bound noted in the main text. Let 𝜙 (𝑀) satisfy 𝜙 (𝑀) → ∞,
log 𝜙 (𝑀) = 𝑜 (log 𝑀),
(81)
choose 𝑘 = ⌈𝑀𝜙 (𝑀)⌉, and assume additionally 𝑀𝜙 (𝑀) log2 4𝑀𝜙 (𝑀) = 𝑜 (𝑁 log2 𝑁 ).
(82)
Theorem 33 (Asymptotic constant-one cache lower bound). Under (81)–(82), every no-recomputation schedule of 𝐺 2 (𝑁 ) satisfies 𝑄 miss (𝑁 , 𝑀, 𝐵) ≥ (1 − 𝑜 (1)) 𝑁𝐵 log𝑀 𝑁 . Proof. Substitute 𝑘 = ⌈𝑀𝜙 (𝑀)⌉ into Theorem 32; by (82) the floor, ceiling, and final incomplete segment are lower-order, so 𝑁 log2 𝑁 𝑀 log2 𝑀 𝑄 miss ≥ (1 − 𝑜 (1)) · · 1− · . 𝐵 log2 𝑀 𝑘 log2 (4𝑘) The choice of 𝑘 gives 1 − 𝑀/𝑘 = 1 − 1/𝜙 (𝑀) + 𝑜 (1) = 1 − 𝑜 (1), while log 𝜙 (𝑀) = 𝑜 (log 𝑀) gives log2 𝑀/log2 (4𝑘) = 1−𝑜 (1); finally log2 𝑁 /log2 𝑀 = log𝑀 𝑁 . □ Thus the normalized miss-only lower-bound constant is at least 1 − 𝑜 (1); this does not assert that PACO or any known schedule attains an upper bound of the form (1 + 𝑜 (1)) 𝑁𝐵 log𝑀 𝑁 .
Disjoint Unions of FFT DAGs
Ã𝑠 Let 𝐺 2(𝑠 ) (𝑛) = 𝑖=1 𝐺 2 (𝑛) (𝑖 ) be the disjoint union of 𝑠 radix-2 FFT DAGs of length 𝑛.
Proof. Contrapositive. Let 𝑈𝑖 = 𝑈 ∩ 𝑉 (𝐺 2 (𝑛) (𝑖 ) ) and 𝑟𝑖 = Í |in(𝑈𝑖 )|, so |in(𝑈 )| = 𝑖 𝑟𝑖 . Suppose |in(𝑈 )| < 𝑘. For every nonempty 𝑈𝑖 , 𝑟𝑖 ≥ 1, and the contrapositive of Theorem 30 with parameter 𝑟𝑖 + 1 gives |𝑈𝑖 | < (𝑟𝑖 + 1) log2 (4(𝑟𝑖 + 1)). Since 𝑟𝑖 ≥ 1, 𝑟𝑖 + 1 ≤ 2𝑟𝑖 and 4(𝑟𝑖 + 1) ≤ 8𝑟𝑖 , so |𝑈𝑖 | < 2𝑟𝑖 log2 (8𝑟𝑖 ) ≤ 2𝑟𝑖 log2 (8𝑘). Summing, Í |𝑈 | < 2( 𝑖 𝑟𝑖 ) log2 (8𝑘) < 2𝑘 log2 (8𝑘). As 2 log2 (8𝑘) ≤ 4 log2 (4𝑘) for 𝑘 ≥ 1, |𝑈 | < 4𝑘 log2 (4𝑘). □ The constant 4 is chosen for a uniform statement and is not claimed tight.
D.6
Cache Lower Bound for a Batch of Local FFTs
Theorem 35 (Finite lower bound for a batch of FFTs). Let one processor evaluate 𝑠 data-disjoint length-𝑛 radix-2 FFT DAGs without recomputation. For every integer 𝑘 > 𝑀, j kl m 𝑠𝑛 log2 𝑛 𝑘 −𝑀 𝑄 batch (𝑛, 𝑠) ≥ 4𝑘 log (4𝑘 . (83) 𝐵 ) 2
Proof. The disjoint union has 𝑠𝑛 log2 𝑛 noninput vertices. Partition into complete segments of 4𝑘 log2 (4𝑘) previously uncomputed vertices; Lemma 34 gives |in(𝑈 )| ≥ 𝑘 per segment, and Lemma 31 gives ⌈(𝑘 − 𝑀)/𝐵⌉ misses each. Sum over segments. □ Corollary 36 (Asymptotic batched lower bound). There are universal constants 𝑐, 𝑐 ′ > 0 with 𝑀 𝑛𝑠 1 + log𝑀 𝑛 − 𝑐 ′ . (84) 𝑄 batch (𝑛, 𝑠) ≥ 𝑐 𝐵 𝐵 If the cache is initially empty, or if 𝑛𝑠 ≥ 𝜅𝑀 for a sufficiently large universal constant 𝜅, then 𝑄 batch (𝑛, 𝑠) = Ω 𝑛𝑠 𝐵 (1 + log𝑀 𝑛) . Proof. Take 𝑘 = 2𝑀 in Theorem 35. Then 𝑘 > 𝑀, and 𝑘 − 𝑀 = 𝑀, 4𝑘 = 8𝑀, so log2 (4𝑘) = 3 + log2 𝑀. Hence 𝑠𝑛 log2 𝑛 𝑀 . 𝑄 batch (𝑛, 𝑠) ≥ 8𝑀 (3 + log2 𝑀) 𝐵 𝑠𝑛 log2 𝑛 . Using ⌈𝑀/𝐵⌉ ≥ 𝑀/𝐵 ≥ 0 and ⌊𝐴⌋ ≥ 8𝑀 (3 + log2 𝑀) 𝐴 − 1 (the product bound holds trivially when 𝐴 < 1, since the left side is nonnegative), 𝑀 𝑀 𝑀 = 𝐴 − . 𝑄 batch ≥ (𝐴 − 1) 𝐵 𝐵 𝐵 For 𝑀 ≥ 2 we have 3 + log2 𝑀 ≤ 4 log2 𝑀, so by the change of base log𝑀 𝑛 = log2 𝑛/log2 𝑀,
Write 𝐴 :=
𝐴
log2 𝑛 𝑀 𝑠𝑛 𝑠𝑛 log𝑀 𝑛 1 𝑠𝑛 = · ≥ · = log𝑀 𝑛. 𝐵 8𝐵 3 + log2 𝑀 8𝐵 4 32 𝐵
1 This gives the stated bound with 𝑐 1 = 32 and 𝑐 2 = 1:
𝑄 batch ≥ 𝑐 1
𝑛𝑠 𝑀 log𝑀 𝑛 − 𝑐 2 . 𝐵 𝐵
Shina Guo, Weiguo Gao, and Yuan Tang
Independently, the 𝑛𝑠 inputs are distinct, so with at most 𝑀 initially resident, 𝑄 batch ≥ (𝑛𝑠 − 𝑀)+ /𝐵. Combining via max{𝑥, 𝑦} ≥ (𝑥 + 𝑦)/2 gives (84). If the cache is empty the compulsory scan contributes Ω(𝑛𝑠/𝐵); if 𝑛𝑠 ≥ 𝜅𝑀 the linear term absorbs the 𝑀/𝐵 term. □
D.7
Application to PACO Local FFT Phases
In Phase I every processor evaluates 𝑠 1 = 𝑁 2 /𝑝 data-disjoint transforms of length 𝑛 1 = 𝑁 1 ; in Phase III it evaluates 𝑠 2 = 𝑁 1 /𝑝 of length 𝑛 2 = 𝑁 2 . Both batches have size 𝑛 1𝑠 1 = 𝑛 2𝑠 2 = 𝑁 /𝑝. Theorem 37 (PACO local cache lower bound). Under the exact slab decomposition and no arithmetic recomputation, PACO’s local = Ω 𝑁 (1 + log 𝑁 ) , and hence two local FFT phases satisfy 𝑄 max 𝑀 𝑝𝐵 so does the complete PACO execution.
Proof. Fix the constant 𝜅 of Corollary 36 (playing the role of the constant 𝜖 in the Proposition 6 roadmap). Large-cache case 𝑁 /𝑝 ≥ 𝜅𝑀. Applying Corollary 36 to each phase 𝑁 𝑁 gives Ω 𝑝𝐵 (1 + log𝑀 𝑁 1 ) and Ω 𝑝𝐵 (1 + log𝑀 𝑁 2 ) ; the phases run sequentially and log𝑀 𝑁 1 + log𝑀 𝑁 2 = log𝑀 𝑁 , so their bounds sum to the claim. Small-cache case 𝑁 /𝑝 < 𝜅𝑀. This is exactly the elementary case of the Proposition 6 roadmap: 𝑝 | 𝑁 1 and 𝑝 | 𝑁 2 give 𝑁 1, 𝑁 2 ≤ 𝑁 /𝑝, hence 𝑁 ≤ (𝑁 /𝑝) 2 < 𝜅 2 𝑀 2 and 1 + log𝑀 𝑁 = 𝑂 (1); Phase I begins with an empty cache and touches Θ(𝑁 /𝑝) distinct inputs, incurring 𝑁 Ω(𝑁 /(𝑝𝐵)) = Ω 𝑝𝐵 (1 + log𝑀 𝑁 ) . The exact slab decomposition assigns identical batch sizes to every processor, so the bound holds for the maximum per-processor cost; the fused middle stage can only add misses. □