FFTE: A High-Performance FFT Library
andParallel Implementation of ClassicalGram-Schmidt Orthogonalization
Daisuke Takahashi
Center for Computational Sciences/Graduate School of Systems and Informaiton Engineering
University of Tsukuba

2007/10/30

External Review1Outline• FFTE: A High-Performance FFT Library– Features, Design and Approach– Performance Results• Parallel Implementation of Classical Gram-Schmidt
Orthogonalization Using Matrix Multiplication(with T. Yokozawa, T. Boku and M. Sato)– Column-wise Blocking CGS Algorithm– Recursive Blocking CGS Algorithm– Performance Results

2007/10/30

External Review

2FFTE: A High-Performance
FFT Library• FFTE is a Fortran subroutine library for computing
the Fast Fourier Transform (FFT) in one or moredimensions.
• It includes complex, mixed-radix and paralleltransforms.
• FFTE is typically faster than other publicallyavailable FFT implementations, and is evencompetitive with vendor-tuned libraries.

2007/10/30

External Review

3

FFTE Library - Features• High speed– Supports Intel’s SSE2/SSE3 instructions.

• Parallel transforms– Shared / Distributed memory parallel computers(OpenMP, MPI and OpenMP + MPI)

• High portability– Fortran77 + OpenMP + MPI– Intel’s intrinsics for SSE2/SSE3 instructions.

• HPC Challenge Benchmark– FFTE’s 1-D parallel FFT routine has been incorporated
into the HPC Challenge (HPCC) benchmark.
2007/10/30

External Review

4FFTE Library - Design• Performance– One goal for large FFTs is to minimize the number ofcache misses.

• Ease of use: routine interfaces– Similar to sequential SGI SCSL or Intel MKL routines• Portability– Communication: MPI– Computation: Fortran77 + OpenMP

2007/10/30

External Review

5

FFTE Library - Approach• Many FFT routines work well when data sets fit intoa cache.
• When a problem size exceeds the cache size,
however, the performance of these FFT routines
decreases dramatically.
• Some previously presented three-dimensional FFT
algorithms require– Three multicolumn FFTs.
– Three data transpositions.
→ The chief bottlenecks in cache-based processors.

• We combine the multicolumn FFTs andtranspositions to reduce the number of cache
misses.
2007/10/30

External Review

6Parallel Block 3-D FFT Algorithm
Partialn2 n3nBTranspose

n1 P0 P1 P2 P3

All-to-Allnn23
Partial comm.
TransposenBn2n2 n1
n3 P0 P1 P2 P3
2007/10/30

n3n1 P0 P1 P2 P3nn23Partial
Transpose

n1External Review

7Performance Results• To evaluate the implemented parallel 3-D FFTs, wecompared– The implemented parallel 3-D FFT,
named FFTE (ver 3.2, supports SSE2, using MPI)
– FFTW (ver. 2.1.5, not support SSE2, using MPI)

• Target parallel machine:– A 16-node dual PC SMP cluster(dual Xeon 2.8GHz, 2GB DDR-SDRAM / node, Linux2.4.20smp).
– Interconnected through a Myrinet-2000 switch.
– MPI-SCore [http://www.pccluster.org] was used.
– We used an intra-node MPI library for the PC SMPcluster.
2007/10/30

External Review

8Performance of parallel 3-D FFTs
（dual Xeon PC SMP cluster, N1xN2xN3 = 2^24xP)

5000
MFLOPS

4000
3000
2000
1000
01
2007/10/30

248
16
Number of Nodes
External Review

FFTE (SSE2,1CPU)FFTE (SSE2,2CPUs)
FFTE (x87,1CPU)FFTE (x87,2CPUs)
FFTW(1CPU)FFTW(2CPUs)
9

1098765432
10Comp.
Comm.

16x1
16x2

8x2

4x2
8x1

2x2
4x12x1.

1x1
1x2Time (sec)Breakdown of parallel 3-D FFTs
(dual Xeon PC SMP Cluster, N1xN2xN3=2^24xP)

Number of Processors
2007/10/30

External Review

10Parallel Implementation of ClassicalGram-Schmidt Orthogonalization• Gram-Schmidt orthogonalization process is oneof the fundamental algorithms in linear algebra.
– Density-funtional theory (DFT) code includes GramSchmidt orthogonalization of a large set of wave
functions.

• Two basic computational variants of the GramSchmidt process exist:– Classical Gram-Schmidt algorithm (CGS)– Modified Gram-Schmidt algorithm (MGS)• Often selected for practical application• More stable than the CGS algorithm
2007/10/30

External Review

11Background• Modified Gram-Schmidt (MGS) algorithm– Cannot be expressed by Level-2 BLAS.
– Parallel implementation requires additionalcommunication.

• Classical Gram-Schmidt (CGS) algorithm– Can be expressed by Level-2 BLAS.
– CGS with DGKS correction [Daniel et al. ‘76]
is one of the most efficient way to perform theorthogonalization process.
2007/10/30

External Review

12Classical Gram-Schmidt Algorithm• The CGS orthogonalization can be performed byusing Level-2 BLAS.
• Suitable for parallelization.
• The CGS algorithm is not stable.
do j = 1, n
qj = aj
do i = 1, j − 1
q j = q j − (q i , a j )q i
end doa j Raw data vector
q j Orthogonalized vectorq j Euclid norm

(q i , a j ) Inner productqj = qj / qj
end do2007/10/30

External Review

13Naive Implementation• CGS for m by n matrix can be computed by– ( m − 1) matrix-vector multiplications (GEMV: Level-2 BLAS)– m normalizations (NRM2 and SCAL: Level-1 BLAS)

Let the matrix A is denoted as A = (a1a 2 L a n ) .
q1 = a1, q1 = q1 / q1
q2 = a2 − (q1, a2 )q1,GEMV

q2 = q2 / q2

q3 = a3 −(q1, a3 )q1 − (q2 , a3 )q2 ,NRM2, SCAL

q3 = q3 / q3

q4 = a4 − (q1, a4 )q1 − (q2, a4 )q2 − (q3, a4 )q3,

q4 = q4 / q4

q5 = a5 − (q1, a5 )q1 − (q2 , a5 )q2 − (q3, a5 )q3 − (q4 , a5 )q4,q5 = q5 / q5

q6 = a6 − (q1, a6 )q1 − (q2 , a6 )q2 − (q3, a6 )q3 − (q4 , a6 )q4 − (q5, a6 )q5,
2007/10/30

External Review

q6 = q6 / q614

CGS using Matrix Multiplication• The CGS orthogonalization of a matrix can bechanged into a matrix multiplication(GEMM: Level-3 BLAS).
q1 = a1,q1 = q1 / q1

q2 = a2 − (q1, a2 )q1,GEMM

q2 = q2 / q2

q3 = a3 −(q1, a3 )q1 − (q2 , a3 )q2 ,GEMV
NRM2, SCAL

q3 = q3 / q3

q4 = a4 − (q1, a4 )q1 − (q2, a4 )q2 − (q3, a4 )q3,q4 = q4 / q4

q5 = a5 − (q1, a5 )q1 − (q2 , a5 )q2 − (q3, a5 )q3 − (q4 , a5 )q4,q5 = q5 / q5

q6 = a6 − (q1, a6 )q1 − (q2 , a6 )q2 − (q3, a6 )q3 − (q4 , a6 )q4 − (q5, a6 )q5,
2007/10/30

External Review

q6 = q6 / q615

Column-wise Blocking CGS
Using Matrix-Vector
Multiplication （GEMV)Using Matrix-Matrix
Multiplication（GEMM)

begin CBCGS( A, Q, N, M , s, h)do j = s, h, M
qj =qj qj

do i = j + 1, j + M
Tw = Qs,s+i a s +i
q s +i = Qs −i,s w
q s +i = q s +i q s +i
end doW = QTj + M,N −( j + M ) A j,j + M

M

M : Blocking SizeN

Q j + M,N −(j + M) = WQ j + M,N −(j + M)
end doend
2007/10/30

M
External Review

N
16

Recursive Blocking CGS
Using Matrix-Vector
Multiplication （GEMV)Using Matrix-Matrix
Multiplication （GEMM)

begin RBCGS( A, Q, N , M , s, h)if (h ≤ M ) thenqs = qs qs
do i = 1, hN

Tw = Qs,s+i a s +iqs+i = Qs−i,sw

Nq s +i = q s +i q s +iN

end doelse

RBCGS( A,Q, N , M , s, h 2)
W = QsT−(h/ 2 ),s +(h/ 2 ) As −(h/ 2 ),s +(h/ 2 )Qs −(h/ 2 ),s +(h/ 2 ) = WQs −(h/ 2 ),s +(h/ 2 )NRBCGS( A,Q,N , M , s + h 2 , h 2)
end ifend
2007/10/30

External Review

8842N2N
17Evaluation Environment• A 32-node Xeon PC-cluster– Xeon 3GHz, 1GB DDR2 Memory• 1000Base-T Gigabit Ethernet• OpenMPI 1.2.3• GotoBLAS r1.19• Intel C Compiler 10.0
• Linux 2.6.20-1.2933.fc6smp• All programs were run in 64-bit mode

2007/10/30

External Review

18Performance on
Xeon 3GHz (1CPU)

7

Naive

Column-wise CGS

Recursive CGS

6G FL O P S

543
2
100
2007/10/30

2000

4000
Matrix Size
External Review

6000

8000
19

Performance on 32 node 3GHzXeonPCCluster
140
Naive
Column-wise CGS
Recursive CGS

120

GFLO PS100
8060
4020002007/10/30

10000

20000
Matrix Size
External Review

30000

40000
20