Implementation of Parallel FFTs onCluster of Intel Xeon Phi Processors
Daisuke Takahashi
Center for Computational Sciences
University of Tsukuba, Japan

2017/12/7

2017 CCS-EPCC Workshop

1Outline• Background• Objectives• Six-Step FFT Algorithm• In-Cache FFT Algorithm and Vectorization• Computation-Communication Overlap• Automatic Tuning of Parallel 1-D FFT• Performance Results• Conclusion
2017/12/7

2017 CCS-EPCC Workshop

2

Background• The fast Fourier transform (FFT) is widely used inscience and engineering.
• Parallel FFTs on distributed-memory parallel
computers require intensive all-to-allcommunication, which affects their performance.
• How to overlap the computation and the all-to-all
communication is an issue that needs to beaddressed for parallel FFTs.
• Moreover, we need to select the optimal
parameters according to the computational
environment and the problem size.
2017/12/7

2017 CCS-EPCC Workshop

3

Objectives• Several FFT libraries with automatic tuning havebeen proposed.
– FFTW, SPIRAL, and UHFFT• An Implementation of parallel 1-D FFT on Xeon Phi(Knights Corner) cluster has been presented [Parket al. 2013].
• However, to the best of our knowledge, parallel 1-D
FFT with automatic tuning on Xeon Phi (KnightsLanding) cluster has not yet been reported.
• We propose an implementation of a parallel 1-D
FFT with automatic tuning on Xeon Phi cluster.
2017/12/7

2017 CCS-EPCC Workshop

4Approach• The parallel 1-D FFT implemented is based on thesix-step FFT algorithm [Bailey 90], which requires
two multicolumn FFTs and three datatranspositions.
• Using this method, we have implemented anautomatic tuning facility for selecting the optimal
parameters of the computation-communication
overlap, the radices, and the block size.

2017/12/7

2017 CCS-EPCC Workshop

5Discrete Fourier Transform (DFT)• 1-D discrete Fourier transform (DFT) isgiven by𝑦𝑦 𝑘𝑘

𝑛𝑛−1

𝑗𝑗𝑗𝑗
= � 𝑥𝑥(𝑗𝑗)𝜔𝜔𝑛𝑛 ,𝑗𝑗=00 ≤ 𝑘𝑘 ≤ 𝑛𝑛 − 1 ,where 𝜔𝜔𝑛𝑛 = 𝑒𝑒 −2𝜋𝜋𝑖𝑖/𝑛𝑛 and 𝑖𝑖 = −1.

2017/12/7

2017 CCS-EPCC Workshop

6

2-D Formulation• If 𝑛𝑛 has factors 𝑛𝑛1 and 𝑛𝑛2 (𝑛𝑛 = 𝑛𝑛1 × 𝑛𝑛2 ), then theindices 𝑗𝑗 and 𝑘𝑘 can be expressed as

𝑗𝑗 = 𝑗𝑗1 + 𝑗𝑗2 𝑛𝑛1 , 𝑘𝑘 = 𝑘𝑘2 + 𝑘𝑘1 𝑛𝑛2 .

• Substituting the indices 𝑗𝑗 and 𝑘𝑘, we derive thefollowing equation:𝑦𝑦 𝑘𝑘2 , 𝑘𝑘1

𝑛𝑛1 −1 𝑛𝑛2 −1

𝑗𝑗2 𝑘𝑘2 𝑗𝑗1 𝑘𝑘2 𝑗𝑗1 𝑘𝑘1= � � 𝑥𝑥 𝑗𝑗1 , 𝑗𝑗2 𝜔𝜔𝑛𝑛2 𝜔𝜔𝑛𝑛1𝑛𝑛2 𝜔𝜔𝑛𝑛1 .
𝑗𝑗1 =0 𝑗𝑗2 =0• An 𝑛𝑛-point FFT can be decomposed into an 𝑛𝑛1 point FFT and an 𝑛𝑛2 -point FFT.
2017/12/7

2017 CCS-EPCC Workshop

7

Six-Step FFT Algorithm• This derivation leads to the following six-step
FFT algorithm [Bailey90]:– Step 1: Transpose– Step 2: Perform 𝑛𝑛1 individual 𝑛𝑛2 -point
multicolumn FFTs

𝑗𝑗1 𝑘𝑘2– Step 3: Perform twiddle factor (𝜔𝜔𝑛𝑛1 𝑛𝑛2 ) multiplication– Step 4: Transpose– Step 5: Perform 𝑛𝑛2 individual 𝑛𝑛1 -point
multicolumn FFTs
– Step 6: Transpose
2017/12/7

2017 CCS-EPCC Workshop

8Parallel 1-D FFT Algorithm Based on
Six-Step
FFT𝑁𝑁
2

GlobalTranspose

𝑁𝑁1𝑃𝑃0𝑃𝑃1𝑃𝑃2 𝑃𝑃3𝑁𝑁1

𝑁𝑁2 𝑃𝑃0 𝑃𝑃1 𝑃𝑃2 𝑃𝑃3
2017/12/7

𝑁𝑁1

𝑁𝑁2

GlobalTranspose

𝑁𝑁12017 CCS-EPCC Workshop

Performtwiddle factor
𝐽𝐽 𝐾𝐾(𝜔𝜔𝑁𝑁11𝑁𝑁22 )multiplication

𝑁𝑁2

GlobalTranspose

9In-Cache FFT Algorithm andVectorization• For in-cache FFT, we used radix-2, 3, 4, 5, 8, 9,
and 16 FFT algorithms based on the mixed-radix
FFT algorithms [Temperton 83].
• Automatic vectorization was used to access theIntel AVX-512 instructions on the Knights Landingprocessor.
• Although higher radix FFTs require more floatingpoint registers to hold intermediate results, theKnights Landing processor has 32 ZMM 512-bit
registers.
2017/12/7

2017 CCS-EPCC Workshop

10Optimization of Parallel 1-D FFT onKnights Landing Processor
COMPLEX*16 X(N1,N2),Y(N2,N1)!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(I,J,JJ)
DO II=1,N1,NB
DO JJ=1,N2,NB
DO I=II,MIN(II+NB-1,N1)
DO J=JJ,MIN(JJ+NB-1,N2)
Y(J,I)=X(I,J)
END DOEND DOTo expand the outermost loop,
END DOthe double-nested loop can be
END DO
collapsed into a single-nested loop.
!$OMP PARALLEL DO
DO I=1,N1
CALL IN_CACHE_FFT(Y(1,I),N2)
END DO…
2017/12/7

2017 CCS-EPCC Workshop

11Computation-Communication
Overlap [Idomura et al. 2014]!$OMP PARALLEL!$OMP MASTER← MPI communication is performedon the master thread!$OMP END MASTER ← No barrier synchronization
MPI communication!$OMP DO SCHEDULE(DYNAMIC)
DO I=1,N← Computation is performed
by a thread other than theEND DO!$OMP DO ← Implicit barrier master thread
DO I=1,N synchronization← Computation is performed
Computation using theafter completion of theresult of communication
END DOMPI communication
Computation!$OMP END PARALLEL
2017/12/7

2017 CCS-EPCC Workshop

12Pipelined ComputationCommunication Overlap
Withoutoverlap

Overlap(NDIV=2)

Computation
Comp.

Communication

Comm.
Comp. Comm.

Overlap(NDIV=4)

2017/12/7

2017 CCS-EPCC Workshop

13Automatic Tuning of Parallel 1-D FFT• The automatic tuning process consists ofthree steps:– Selection of the number of divisions NDIV for thecomputation-communication overlap– Selection of the radices (𝑁𝑁1 and 𝑁𝑁2 )– Selection of the block size NB

2017/12/7

2017 CCS-EPCC Workshop

14

Selection of Number of Divisions for
Computation-Communication Overlap• When the number of divisions for computationcommunication overlap is increased, the overlap ratio
also increases.
• On the other hand, the performance of all-to-all
communication decreases due to reducing themessage size.
• Thus, a tradeoff exists between the overlap ratio andthe performance of all-to-all communication.
• The default overlapping parameter of the original FFTE
6.2alpha is NDIV=4.
• In our implementation, the overlapping parameter
NDIV is varied between 1, 2, 4, 8 and 16.
2017/12/7

2017 CCS-EPCC Workshop

15

Selection of Radices• As long as condition 𝑁𝑁 = 𝑁𝑁1 𝑁𝑁2 is satisfied, we
can select 𝑁𝑁1 and 𝑁𝑁2 arbitrarily for 𝑁𝑁1 ,𝑁𝑁2 ≥ 𝑃𝑃.
• We need to select the best combination andorder of 𝑁𝑁1 and 𝑁𝑁2 for computing a parallel 1-DFFT.
• If 𝑁𝑁 and 𝑃𝑃 are a power of two, 𝑁𝑁1 is varied
between 𝑃𝑃, 2𝑃𝑃, … , 𝑁𝑁, then 𝑁𝑁2 = 𝑁𝑁/ 𝑁𝑁1 .
• In this case, the size of search space islog 2 ( 𝑁𝑁/𝑃𝑃).
2017/12/7

2017 CCS-EPCC Workshop

16

Selection of Block Size
• The default blocking parameter of the original
FFTE 6.2alpha is NB=32.
• Although the optimal block size may dependon the problem size, the block size NB canalso be varied.
• In our implementation, the block size NB isvaried between 4, 8, 16, 32 and 64.

2017/12/7

2017 CCS-EPCC Workshop

17Performance Results• To evaluate the parallel 1-D FFT with automatic tuning (AT),we compared its performance with that of the FFTW 3.3.6pl1, the FFTE 6.2alpha (http://www.ffte.jp/) and the FFTE
6.2alpha with AT.
• The performance was measured on the Oakforest-PACS at
Joint Center for Advanced HPC (JCAHPC).
– 8208 nodes, Peak 25.008 PFlops– CPU: Intel Xeon Phi 7250 (68 cores, Knights Landing 1.4 GHz)– Interconnect: Intel Omni-Path Architecture– Compiler: Intel Fortran compiler 17.0.1.132– Compiler option: “mpiifort -O3 -xMIC-AVX512 -qopenmp”– MPI library: Intel MPI 2017.1.132– flat/quadrant, MCDRAM only, KMP_AFFINITY=compact– Each MPI process has 64 cores and 64 threads.

2017/12/7

2017 CCS-EPCC Workshop

18

Results of automatic tuning of parallel 1-D FFTs
(Oakforest-PACS, 128 nodes)
FFTE 6.2alpha

FFTE 6.2alpha with AT

N

N1

N2

NB NDIV GFlops

N1

N2

NB NDIV GFlops

16M

4K

4K

32

457.8

4K

4K

32

1109.4

64M

8K

8K

32

4116.98K

8K

162154.8

256M

16K

16K

32

473.3

8K

32K

16

1513.8

1G32K32K

324541.7

32K

32K

64

4554.9

4G

64K

64K

324217.0

64K

64K

32

16516.5

2017/12/7

2017 CCS-EPCC Workshop

19

700
600
500
400
300
200
10003230

282624FFTE
6.2alpha (nooverlap)
FFTE
6.2alpha(NDIV=4)FFTE
6.2alpha withAT
FFTW 3.3.6pl1

2220

GFlopsPerformance of parallel 1-D FFTs
（Oakforest-PACS，128 nodes）

Length of transform (log_2 N)
2017/12/7

2017 CCS-EPCC Workshop

20

2500
MPI_Alltoall

2000
1500
1000
500

4M

64
K
512K

8K

1K

1280

16

Bandwidth (MB/sec)

Performance of all-to-all communication（Oakforest-PACS，128 nodes）

Message size (bytes)
2017/12/7

2017 CCS-EPCC Workshop

21

Breakdown of execution time in FFTE 6.2alpha (nooverlap, Oakforest-PACS, N=2^26×number of nodes)

3Communication

Time (sec)

2.5Computation

2
1.5
1
0.5
0
12

4

8

1632 64 128

Number of nodes
2017/12/7

2017 CCS-EPCC Workshop

22

Conclusion• We proposed an implementation of parallel 1-D FFTwith automatic tuning on Xeon Phi cluster.
• We used a computation-communication overlap
method that introduces a communication thread withOpenMP.
• An automatic tuning facility for selecting the optimal
parameters of the computation-communication overlap,
the radices, and the block size was implemented.
• The performance results demonstrate that theproposed implementation of a parallel 1-D FFT with
automatic tuning is efficient for improving theperformance on Xeon Phi cluster.
2017/12/7

2017 CCS-EPCC Workshop

23