Energy Efficient, High-Performance Solvers through Small Dense - - PowerPoint PPT Presentation
Energy Efficient, High-Performance Solvers through Small Dense Matrix Computations on GPUs Azzam Haidar and Stan Tomov with J. Dongarra, P. Luszczek, I. Yamazaki, M. Gates, H. Anzt, and T. Dong University of Tennessee, Knoxville GPU
with J. Dongarra, P. Luszczek, I. Yamazaki, M. Gates, H. Anzt, and T. Dong
GPU Technology Conference March 17 – 20, 2015 San Jose, CA
– Hybrid GPU+CPU dense linear algebra – Sparse linear algebra – Batched linear algebra – Power considerations and GPU-only algorithms
Comp mput ute int ntens nsive (dense LU) vs. Memory-bound
mputa utati tion n (SpMV)
– Tiled linear algebra algorithms – Multifrontal methods – Preconditioners (using DLA) in sparse iterative solvers, many applications, …
LU, QR, or Cholesky
Sparse / Dense Matrix System
TRSMs, QRs, or LUs
TRSMs, TRMMs
Updates (Schur complement) GEMMs, SYRKs, TRMMs DAG-based factorization Batched LA
And many other BLAS/LAPACK, e.g., for application specific solvers, preconditioners, and matrices
– Matrix Algebra for GPU and Multicore Architecture – LAPACK/ScaLAPACK
– http://icl.cs.utk.edu/magma/
– UTK, UC Berkeley, UC Denver, INRIA (France), KAUST (Saudi Arabia) – Community effort, similarly to LAPACK/ScaLAPACK
MAGMA Hybrid Algorithms
(heterogeneity friendly) Rely on
(for nested parallelism)
HIBRID ALGORITHMS
MAGMA uses hybrid algorithms where the computation is split into tasks of varying granularity and their execution scheduled over the hardware components. Scheduling can be static or dynamic. In either case, small non-parallelizable tasks,
PERFORMANCE & ENERGY EFFICIENCY
– Representing linear algebra algorithms as collections
– Properly scheduling tasks' execution over multicore and GPU hardware components
– One- and two-sided factorizations and solvers – Iterative linear and eigensolvers
– 1) High level; 2) Leveraging prior developments; 3) Exceeding in performance homogeneous solutions
Hybrid CPU+GPU algorithms
(small tasks for multicores and large tasks for GPUs)
for( j=0, j<n; j+=nb) { jb = min(nb, n-j); zherk( MagmaUpper jb, j, one, dA(0 if (j+jb < n) zgemm( MagmaCo dA(0,j), ldd zpotrf( MagmaUpper if (info != 0) *info += j; If (j+jb) < n) { ztrsm( MagmaLeft, jb, n } }
Left-looking hybrid Cholesky
LAPACK
MAGMA
From sequential to parallel hybrid
representation of MAGMA algorithms almost identical to LAPACK
Note:
1 for( j=0, j<n; j+=nb) { 2 jb = min(nb, n-j); 3 magma_zherk( MagmaUpper, MagmaConjTrans, jb, j, one, dA(0,j), ldda, one, dA(j,j), ldda, queue); 4 magma_zgetmatrix_async( jb, jb, dA(j,j), ldda, work, jb, queue, &event); 5 if (j+jb < n) 6 magma_zgemm( MagmaConjTrans, MagmaNoTrans, jb, n-j-jb, j, one, dA(0,j), ldda, dA(0,j+jb), ldda, one, dA(j, j+jb), ldda, que 7 magma_event_sync( event ); 8 zpotrf( MagmaUpperStr, &jb, work, &jb, info); 9 if (info != 0) 10 *info += j; 11 magma_zsetmatrix_async(jb, jb, work, jb, dA(j, j), ldda, queue, &event); 12 If (j+jb) < n) { 13 magma_event_sync( event ); 14 magma_ztrsm( MagmaLeft, MagmaUpper, MagmaConjTrans, MagmaNo jb, n-j-jb, one, dA(j,j), ldda, dA(j,j+jb), ldda, queue); } }
1 for( j=0, j<n; j+=nb) { 2 jb = min(nb, n-j); 3 magma_zherk( MagmaUpper, MagmaConjTrans, jb, j, one, dA(0,j), ldda, one, dA(j,j), ldda, queue); 4 magma_zgetmatrix_async( jb, jb, dA(j,j), ldda, work, jb, queue, &event); 5 if (j+jb < n) 6 magma_zgemm( MagmaConjTrans, MagmaNoTrans, jb, n-j-jb, j, one, dA(0,j), ldda, dA(0,j+jb), ldda, one, dA(j, j+jb), ldda, que 7 magma_event_sync( event ); 8 zpotrf( MagmaUpperStr, &jb, work, &jb, info); 9 if (info != 0) 10 *info += j; 11 magma_zsetmatrix_async(jb, jb, work, jb, dA(j, j), ldda, queue, &event); 12 If (j+jb) < n) { 13 magma_event_sync( event ); 14 magma_ztrsm( MagmaLeft, MagmaUpper, MagmaConjTrans, MagmaNo jb, n-j-jb, one, dA(j,j), ldda, dA(j,j+jb), ldda, queue); } }
for( j=0, j<n; j+=nb) { jb = min(nb, n-j); zherk( MagmaUpper jb, j, one, dA(0 if (j+jb < n) zgemm( MagmaCo dA(0,j), ldd zpotrf( MagmaUpper if (info != 0) *info += j; If (j+jb) < n) { ztrsm( MagmaLeft, jb, n } }
Left-looking hybrid Cholesky
LAPACK
MAGMA
From sequential to parallel hybrid
representation of MAGMA algorithms almost identical to LAPACK
Note: MAGMA runtime environment
computation to the GPU
1 for( j=0, j<n; j+=nb) { 2 jb = min(nb, n-j); 3 magma_zherk( MagmaUpper, MagmaConjTrans, jb, j, one, dA(0,j), ldda, one, dA(j,j), ldda, queue); 4 magma_zgetmatrix_async( jb, jb, dA(j,j), ldda, work, jb, queue, &event); 5 if (j+jb < n) 6 magma_zgemm( MagmaConjTrans, MagmaNoTrans, jb, n-j-jb, j, one, dA(0,j), ldda, dA(0,j+jb), ldda, one, dA(j, j+jb), ldda, que 7 magma_event_sync( event ); 8 zpotrf( MagmaUpperStr, &jb, work, &jb, info); 9 if (info != 0) 10 *info += j; 11 magma_zsetmatrix_async(jb, jb, work, jb, dA(j, j), ldda, queue, &event); 12 If (j+jb) < n) { 13 magma_event_sync( event ); 14 magma_ztrsm( MagmaLeft, MagmaUpper, MagmaConjTrans, MagmaNo jb, n-j-jb, one, dA(j,j), ldda, dA(j,j+jb), ldda, queue); } }
for( j=0, j<n; j+=nb) { jb = min(nb, n-j); zherk( MagmaUpper jb, j, one, dA(0 if (j+jb < n) zgemm( MagmaCo dA(0,j), ldd zpotrf( MagmaUpper if (info != 0) *info += j; If (j+jb) < n) { ztrsm( MagmaLeft, jb, n } }
Left-looking hybrid Cholesky
LAPACK
MAGMA
14’ 14’ 14’ 14’ 14 13 11 7 6 4 3’ 3’ 3’ 3’ 3
CPU GPU
7 13 8
CUDA Queues
3 4 6 7 11 13 14 4 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 6’ 11 6’ 6’ 6’ 6’
PCI time
Tasks
6 4 3 7 14 13 11
Offloaded to the GPU Offloaded to the GPU
Computed
From sequential to parallel hybrid
representation of MAGMA algorithms almost identical to LAPACK
Note:
CPU task #8 and CPU-GPU communications are overlapped with GPU computations
MAGMA runtime environment
100 200 300 400 500 600 700 SP Solve DP Solve (MP Iter.Ref.) DP Solve
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores)
Solving general dense linear systems using mixed precision iterative refinement
50 100 150 200 250 300 350 DGETRF
Solving large problems that do not fit in the GPU memory
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores) Matrices of size that do not fit in a specified GPU memory
50 100 150 200 250 300 350 DGETRF
Solving large problems that do not fit in the GPU memory
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores) Out-of-GPU-memory Algorithms can now solve large problems
1) Copy A2 to the GPU 2) Update A2 using A1 (a panel of A1 at a time) 3) Factor the updated A2 using existing hybrid code 4) Copy factored A2 to the CPU Trivially extended to multiGPUs: A2 is “larger” with 1-D block cyclic distribution, again reusing existing algorithms Factored sub-matric A1 on CPU To be factored sub-matrix A2 on GPU . . .
Untouched part of the matrix
GPU GPU 1 GPU 2 GPU
nb
100 200 300 400 500 600 700 800 900 1000 1 GPU CPU (MKL)
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores)
100 200 300 400 500 600 700 800 900 1000 2 GPUs 1 GPU CPU (MKL)
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores)
100 200 300 400 500 600 700 800 900 1000 3 GPUs 2 GPUs 1 GPU CPU (MKL)
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores)
100 200 300 400 500 600 700 800 900 1000 Kepler (K20X) 3 GPUs 2 GPUs 1 GPU CPU (MKL)
Keeneland GPU M2090 (14 MP @1.3 GHz, peak 583 GFlop/s) CPU Intel Xeon X5660@2.80GHz (2 x 6 cores)
200 400 600 800 1000 1200 1400 1600 1800 2000 5000 10000 15000 20000 25000 30000 Gflop/s Matrix Size N X N 1 K80 .5 K80
Host
Ivy Bridge (2 x 10 cores @3.0 GHz) DP Peak 480 GFlop/s DGEMM Peak ~ 430 GFlop/s
GPUs
K80 ( 2 x 1664 CUDA cores @823.5 MHz) DGEMM Peak ~ 1000 GFlop/s
95% of dgemm peak 89% of dgemm peak
500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 10000 20000 30000 40000 50000 60000 Gflop/s Matrix Size N X N 4 K80 2 K80 1 K80 .5 K80
Host
Ivy Bridge (2 x 10 cores @3.0 GHz) DP Peak 480 GFlop/s DGEMM Peak ~ 430 GFlop/s
GPUs
K80 ( 2 x 1664 CUDA cores @823.5 MHz) DGEMM Peak ~ 1000 GFlop/s
1,530 Gflop/s per K80
– Quantum mechanics (Schrödinger equation) – Quantum chemistry – Principal component analysis (in data mining) – Vibration analysis (of mechanical structures) – Image processing, compression, face recognition – Eigenvalues of graph, e.g., in Google’s page rank . . .
first stage second stage
Characteristics
transformation extra cost.
Keeneland system, using one node 3 NVIDIA GPUs (M2090@ 1.1 GHz, 5.4 GB) 2 x 6 Intel Cores (X5660 @ 2.8 GHz, 23 GB)
flops formula: 4n3/3*time
Higher is better
Acceleration w/ 3 GPUs:
15 X vs. 12 Intel cores
Towards extreme-scale all-electron quantum simulations for materials design. In preparation 2015
A novel hybrid CPU-GPU generalized eigensolver for electronic structure calculations based on fine grained memory aware tasks, International Journal of High Performance Computing Applications , vol. 28, no 2, pp. 196—209, May 2014.
Leading edge multi-GPU algorithms for generalized eigenproblems for electronic structure calculations. International Supercomputing Conference IEEE-ISC 2013.
Hybrid-GPU is 1.5 times faster than CPU Hybrid-GPU is 2 times more energy efficient
Cray XC30:
n = 16000,
H = Q1T A Q1
T = Q2T H Q2
and back-transform to eigenvectors X of A
A H T Z X
Accelerating computation of eigenvectors in the nonsymmetric eigenvalue problem. VecPar 2014 the 11th International Meeting High Performance Computing for Computational Science.
Toward a scalable multi-GPU eigensolver via compute-intensive kernels and efficient communication. International Conference on Supercomputing, ICS 2013, Oregon, USA
matrix of size 40K using six GPUs K20c.
implementation is able to schedule and balance the tasks on the GPU devices (six).
Trace of Dynamic MAGMA DPOTRF using QUARK on six GPUs
Unified Development for Mixed Multi-GPU and Multi-Coprocessor Environments using a Lightweight Runtime Environment. IEEE IPDPS 2014, 28th IEEE International Parallel and Distributed Processing Symposium.
Model-Driven One-Sided Factorizations on Multicore Accelerated Systems. Supercomputing frontiers and innovations journal 2014.
– ScaLAPACK 2D block-cyclic data distribution – Lightweight “local” (node) scheduling with QUARK + MPI communications – Match in performance previous results using “tile” algorithms
200 400 600 800 1000 1200
Kernel peakx(#cores+#GPUs) Cholesky factorization
20 40 60 80 100 120 1 2 4 8 16 32 64 100
Nu Number of Nodes DGEMM UB
Single node performance ( N= 34 K) Weak scaling ( N= 340 K)
75% 75%
Systems”, ACM International Conference on Supercomputing (ICS 2012), San Servolo Island, Venice, Italy, June 2012.
SpMM on block of 32 vectors in double precision
(1) Replacing GMRES’ SpMV Matrix Powers Kernel (MPK): vk+1 = A vk for k = j, ... , j+s (2) BLAS-2 BLAS-3 based orthogonalization (next …)
MPK to gener erate e 100 vec ector
vario ious s s Over erall all perfor
mance imp mprovemen ment on up to 3 GPUs
Multiple GPUs”, Proc. of IPDPS 2014, Phoenix, Arizona, May 19—23, 2014. SpMV Orth
2 X
Acceleration
– CholQR obtains BLAS-3 performance, but error is bounded by εκ(V)2 – Remove the “square” by selectively using double-double (doubled) precision
CA-GMRES on GPUs”, VECPAR 2014 (Best paper award), Eugene, Oregon, USA, June 30—July 3, 2014.
Motivation: power and performance advantages
GPU: NVIDIA K40c CPU: Intel Sandy Bridge 15 MP x 192 @ 0.88 GHz 2 x 8 cores @ 2.6 GHz DP Peak: 1690 GFlop/s DP Peak: 332.8 Gflop/s
=> to minimi mize energy y consum sumption, ption, use the entirely ely GPU U imp mpleme ement ntat ation ion
MAGMA QR in double precision on a K40 GPU: GPU only vs. Hybrid
Motivation: power and performance advantages
QR factorization in single precision on the Jetson TK1
Mot
ivat ation:
We present here a feasibility design study, the idea is to target the new high end technologies. Obser servati tions
urrent rent situati uation
for GPU
have often been used in combination with CPUs, where the CPU handles the small and difficult tasks to be parallelized
We present here a feasibility design study, the idea is to target the new high- end technologies. Our ur goals: s:
algebra subroutines tuned for the modern processor architecture and system designs and that outperform multicore CPUs in both performance and energy efficiency.
model of the new emerging accelerators and coprocessors.
This will provide the developers of applications, compilers, and runtime systems with the option of expressing computation as a loop, or a single call to a routine from the new batch operation standard.
Al Algori rith thmic mic basics: cs:
er Ax=b follow the Lapack style algorithmic design blocking algorithm
Hardw dware re charac acteri rist stics cs and limit: t:
is 32 threads)
P a n e l Pi Trailing matrix update Factored part of A
Classical sical strat ategies egies design sign
handle it Difficulties culties
becomes expensive Proposition position
Classical sical strat ategies egies design sign
Difficulties culties
phases to be efficient Proposition position
Classical sical strat ategies egies design sign
shared memory – load it with data and reuse that data in computations as much as possible. Difficulties culties
performance for classical GPU kernels but is not that appealing for batched algorithm for different reasons.
Difficulties culties
performance of memory bound operations, since only one thread- block will be mapped to that SMX at a time (low occupancy)
threads used in the thread block will be limited, resulting in low
Proposition: position: cus ustom
sign per operat rations ions type e
swap kernel 60% gemm kernel 15%
swap kernel 60% gemm kernel 15%
swap kernel 60% gemm kernel 15%
swap kernel 60% gemm kernel 15%
swap kernel 60% gemm kernel 15%
swap kernel 60% gemm kernel 15% Bottlenecks:
matrix is stored in transpose form. However if matrix is stored in transpose form the swap is fast BUT the other components become very slow. Proposition:
dlaswp, so one optimization is to use shared memory to load a chunk of the nb rows
swap kernel 60% gemm kernel 15% gemm kernel 30% swap kernel 10%
panel: classical getf2 38%
P a n e L Trailing matrix update Factored part of A
32
panel: classical getf2 38% Bottlenecks:
efficient since dealing with gemm’s of small sizes
Proposition:
nested blocking technique that block also the panel.
P a n e L Trailing matrix update Factored part of A
32
panel: classical getf2 38%
panel: classical blocked getf2 8%
batched dgemm
batched dgemm
batched dgemm
batched dgemm
Bottlenecks:
are well suited for small matrix sizes (128) but stagnate for larger sizes (>128) Proposition:
for large matrix size (>128) and thus we propose to use both streamed and batched according to the size of the trailing matrix
batched dgemm streamed dgemm
Comparison with CPU:
after matrix, where each factorization is using the multi-thread version of the MKL Library. Expected to have low in performance because each matrix is small – it does not exhibit parallelism and so the multithreaded code is not able to feed with work all 16 SB threads used.
matrices are small (< 512) and at least 16 of them fit in the L3 cache level. One of the best technique is to use each thread to factorize independently a matrix. This way 16 factorizations are conducted independently in parallel.
High gher is bet etter
High gher is bet etter
High gher is bet etter
dgeqrf of 1000 batched matrices of size 1024x1024
CPU does es not t include de DRAM AM power er
– Simulates evolution of the nuclear kinetics where for any single time step on a single zone there is need to solve a small problem with LU – LU factorizations, including the backward substitution, account for 75% of the time – Zones can be solved independently – Number of zones can grow with domain size and dimension to tens of thousands – Can set for example to have 3900 zones per GPU – Problem sizes are small, in the range 14 to 150
Titan supercomputer at ORNL – hybrid nodes with CPU : 16 AMD Opteron Cores GPU : Nvidia K20x
Tomov, and J. Dongarra. LU Factorization of Small Matrices: Accelerating Batched DGETRF on the GPU. HPCC 2014, The 16th IEEE International Conferences
References
Batched Matrix Computations on Hardware Accelerators Based on GPUs. International Journal of High Performance Computing Applications 2014.
Optimization for Performance and Energy for Batched Matrix Computations on GPUs 20th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, Workshop on General Purpose Processing Using
LU Factorization of Small Matrices: Accelerating Batched DGETRF on the GPU. HPCC 2014, The 16th IEEE International Conferences on high performance computing and communications.
A Fast Batched Cholesky Factorization on a GPU ICPP 2014, The 43rd International Conference on Parallel Processing 2014.
– Variable sizes – Mixed-precision techniques – Sparse direct multifrontal solvers & preconditioners – Applications
– autotuning
Multicore Architectures] team http://icl.cs.utk.edu/magma/
Scalable Multicore Architectures] team http://icl.cs.utk.edu/plasma
– University of Tennessee, Knoxville – University of California, Berkeley – University of Colorado, Denver – INRIA, France – KAUST, Saudi Arabia