CPU Levels in the memory hierarchy Level 1 Level 2... Increasing distance from the CPU in access time Level n Size of the memory at each level 1: 2.2

Size: px

Start display at page:

Download "CPU Levels in the memory hierarchy Level 1 Level 2... Increasing distance from the CPU in access time Level n Size of the memory at each level 1: 2.2"

まさとしうづき
5 years ago
Views:

1 FFT 1 Fourier fast Fourier transform FFT FFT FFT 1 FFT FFT 2 Fourier 2.1 Fourier FFT Fourier discrete Fourier transform DFT DFT n 1 y k = j=0 x j ω jk n, 0 k n 1 (1) x j y k ω n = e 2πi/n i = 1 (1) n DFT O(n 2 ) FFT O(n log n) n DFT FFT [4, 16] FFT Cooley-Tukey [6] Stockham [5, 13] スーパーコンピューティングニュース

2 CPU Levels in the memory hierarchy Level 1 Level 2... Increasing distance from the CPU in access time Level n Size of the memory at each level 1: 2.2 FFT FFT [16, 4] FFT FFT radix p p 1 Y (k) = X(j)Ω j ωp jk (2) j=0 Ω twiddle factor [4] 1 p FFT X(j) Ω j p DFT[10] (2) [12, 15] memory hierarchy 1 locality スーパーコンピューティングニュース

3 CPU L1 Cache L2 Cache Main Memory 2: RISC RISC L1 Cache L2 Cache 3 スーパーコンピューティングニュース

4 C SUBROUTINE ZAXPY(N,A,X,Y) IMPLICIT REAL*8 (A-H,O-Z) COMPLEX*16 A,X(*),Y(*) DO I=1,N Y(I)=Y(I)+A*X(I) END DO RETURN END 3: ZAXPY ZAXPY FFT FFT ZAXPY A X plus Y 4 Intel Xeon 3.06 GHz FSB 533 MHz 512 KB L2 cache PC2100 DDR- SDRAM Intel C Compiler 8.0 Intel Pentium4 SIMD Single Instruction Multiple Data SSE2 [8] x87 Intel MKL Math Kernel Library Version [9] BLAS Basic Linear Algebra Subprograms ZAXPY 3 1 iteration 4 load 4 store 2 4 L2 N 8192 SSE2 with SSE2 3 GFLOPS Xeon 3.06 GHz 6.12 GFLOPS L2 x87 4 Six-Step FFT six-step FFT [3, 16] six-step FFT FFT スーパーコンピューティングニュース

5 スーパーコンピューティングニュース

6 スーパーコンピューティングニュース

7 1 COMPLEX*16 X(N1,N2),Y(N2,N1),U(N2,N1) 2 DO I=1,N1 3 DO J=1,N2 4 Y(J,I)=X(I,J) 5 END DO 6 END DO 7 DO I=1,N1 8 CALL IN CACHE FFT(Y(1,I),N2) 9 END DO 10 DO I=1,N1 11 DO J=1,N2 12 Y(J,I)=Y(J,I)*U(J,I) 13 END DO 14 END DO 15 DO J=1,N2 16 DO I=1,N1 17 X(I,J)=Y(J,I) 18 END DO 19 END DO 20 DO J=1,N2 21 CALL IN CACHE FFT(X(1,J),N1) 22 END DO 23 DO I=1,N1 24 DO J=1,N2 25 Y(J,I)=X(I,J) 26 END DO 27 END DO 5: six-step FFT 6 six-step FFT 7 6 NB NP WORK 7 X WORK Y 1 16 WORK WORK X WORK multicolumn FFT six-step FFT two-pass [3, 16] six-step FFT n FFT O(n log n) O(n) Step 2 Step 4 column FFT L1 n column FFT L1 [1, 2] column FFT L1 column FFT two-pass three-pass FFT six-step FFT スーパーコンピューティングニュース

8 1 COMPLEX*16 X(N1,N2),Y(N2,N1),U(N1,N2) 2 COMPLEX*16 WORK(N2+NP,NB) 3 DO II=1,N1,NB 4 DO JJ=1,N2,NB 5 DO I=II,II+NB-1 6 DO J=JJ,JJ+NB-1 7 WORK(J,I-II+1)=X(I,J) 8 END DO 9 END DO 10 END DO 11 DO I=1,NB 12 CALL IN CACHE FFT(WORK(1,I),N2) 13 END DO 14 DO J=1,N2 15 DO I=II,II+NB-1 16 X(I,J)=WORK(J,I-II+1)*U(I,J) 17 END DO 18 END DO 19 END DO 20 DO JJ=1,N2,NB 21 DO J=JJ,JJ+NB-1 22 CALL IN CACHE FFT(X(1,J),N1) 23 END DO 24 DO I=1,N1 25 DO J=JJ,JJ+NB-1 26 Y(J,I)=X(I,J) 27 END DO 28 END DO 29 END DO 6: six-step FFT out-of-place Stockham [5, 13] Step 2 4 multicolumn FFT O( n) FFT Step O( n) WORK 6 In-Cache FFT multicolumn FFT column FFT in-cache FFT Stockham [5, 13] Stockham Cooley-Tukey [6] Cooley-Tukey 2 [6] 2 Stockham n = 2lm l m 2 l n/2 スーパーコンピューティングニュース

9 1. Partial transpose NB 2. NB individual N2-point FFTs NB NB N N Array X N Array WORK N padding NP 3. Partial transpose N2 NB N1 Array WORK N Array X NB 4. NB individual N1-point FFTs N1 NB N2 Array X 7: six-step FFT 2 m 1 2 X Y X Y Y X ω p = e 2πi/p c 0 = X(k + jm) c 1 = X(k + jm + lm) Y (k + 2jm) = c 0 + c 1 Y (k + 2jm + m) = ω j 2l (c 0 c 1 ) 0 j < l 0 k < m 2 FFT 2 FFT 4 8 FFT 2 FFT [14] n = 2 p (p 2) FFT n = 4 q 8 r (0 q 2, r 0) 4 8 FFT n 4 2 FFT six-step FFT multicolumn FFT [3, 16] 5 six-step FFT DO OpenMP[11]!$OMP DO スーパーコンピューティングニュース

10 1: Intel Xeon GHz FFTE 4.0 SSE3 n 1 CPU, 1 core 1 CPU, 2 cores 2 CPUs, 4 cores Time MFLOPS Time MFLOPS Time MFLOPS six-step FFT 3 20 DO WORK MPI FFT 7 Six-Step FFT six-step FFT FFT FFTE version FFT FFTW version [7] n = 2 m m FFT 10 FFT Intel Xeon GHz 4 GB DDR2- SDRAM 32 KB L1 instruction cache 32 KB L1 data cache 4 MB L2 Cache Linux fc6 Intel Fortran version 9.1 Intel C version 9.1 -O3 -xp -openmp 1 FFTE version 4.0 FFTW version six-step FFT FFTE 2CPUs 4cores n FFTW 8 FFT スーパーコンピューティングニュース

11 スーパーコンピューティングニュース

12 [8] Intel Corporation. IA-32 Intel Architecture Software Developer s Manual Volume 2: Instruction Set Reference, [9] Intel Corporation. Intel Math Kernel Library Reference Manual, [10] H. J. Nussbaumer. Fast Fourier Transform and Convolution Algorithms. Springer-Verlag, New York, second corrected and updated edition, [11] OpenMP. Simple, portable, scalable smp programming. [12] R. C. Singleton. An algorithm for computing the mixed radix fast Fourier transform. IEEE Trans. Audio Electroacoust., 17:93 103, [13] P. N. Swarztrauber. FFT algorithms for vector computers. Parallel Computing, 1:45 63, [14] D. Takahashi. A parallel 1-D FFT algorithm for the Hitachi SR8000. Parallel Computing, 29(6): , [15] C. Temperton. Self-sorting mixed-radix fast Fourier transforms. J. Comput. Phys., 52:1 23, [16] C. Van Loan. Computational Frameworks for the Fast Fourier Transform. SIAM Press, Philadelphia, PA, スーパーコンピューティングニュース

untitled

untitled A = QΛQ T A n n Λ Q A = XΛX 1 A n n Λ X GPGPU A 3 T Q T AQ = T (Q: ) T u i = λ i u i T {λ i } {u i } QR MR 3 v i = Q u i A {v i } A n = 9000 Quad Core Xeon 2 LAPACK (4/3) n 3 O(n 2 ) O(n 3 ) A {v i }