Report

Solving Challenging Numerical Linear Algebra Algorithms Using GPU Accelerators Hatem Ltaief KAUST Supercomputing Laboratory Stanimire Tomov University of Tennessee, Knoxville ISC’12 Tutorial, Hamburg June 17, 2012 Outline MAGMA: LAPACK for GPUs Methodology overview • Use both GPUs and multicore CPUs MAGMA: from single to multiGPU support • One-sided factorizations and linear solvers • Two-sided factorizations and eigensolvers Dynamic scheduling approaches to DLA MAGMA algorithms with dynamic scheduling Conclusions MAGMA: LAPACK for GPUs MAGMA • Matrix algebra for GPU and multicore architecture • To provide LAPACK/ScaLAPACK on hybrid architectures • http://icl.cs.utk.edu/magma/ MAGMA BLAS • A subset of BLAS for GPUs, highly optimized for NVIDIA GPGPUs • Fast GEMM for Fermi (CUBLAS3.2) [IJHPCA’10] MAGMA developers & collaborators • UTK, UC Berkeley, UC Denver, INRIA (France), KAUST (Saudi Arabia) • Community effort, similarly to LAPACK/ScaLAPACK A New Generation of DLA Software MAGMA Hybrid Algorithms (heterogeneity friendly) Rely on - hybrid scheduler - hybrid kernels MAGMA Software Stack MAGMA 1.1 50+ hybrid LAPACK algorithms have been developed (total of 200+ routines) Every algorithm is in 4 precisions (s/c/d/z) There are 3 mixed precision algorithms (zc & ds) These are hybrid algorithms, expressed in terms of BLAS MAGMA BLAS A subset of GPU BLAS, optimized for Tesla and Fermi GPUs MAGMA Methodology A methodology to use all available resources: MAGMA uses HYBRIDIZATION methodology based on – Representing linear algebra algorithms as collections of TASKS and DATA DEPENDENCIES among them – Properly SCHEDULING tasks' execution over Hybrid CPU+GPU algorithms (small tasks for multicores and large tasks for GPUs) multicore and GPU hardware components Successfully applied to fundamental linear algebra algorithms – One and two-sided factorizations and solvers – Iterative linear and eigen-solvers Productivity – 1) High-level; 2) Leveraging prior developments; 3) Exceeding in performance homogeneous solutions Hybrid Algorithms One-sided factorizations (LU, QR, Cholesky) Hybridization – Panels (Level 2 BLAS) are factored on CPU using LAPACK – Trailing matrix updates (Level 3 BLAS) are done on the GPU using “look-ahead” A Hybrid Algorithm Example Left-looking hybrid Cholesky factorization in MAGMA 1.0 The difference with LAPACK – the 3 additional lines in red Line 10 (done on CPU) is overlapped with work on the GPU LU Factorization (Single GPU) From single to multiple GPUs support Data distribution • 1-D block-cyclic distribution nb Algorithm • GPU holding current panel is sending it to CPU • All updates are done in parallel on the GPUs • Look-ahead is done with GPU holding the next panel GPU GPU GPU GPU 0 1 2 0 ... LU factorization (multiGPUs) LU factorization (multiGPUs) Matrix out of GPU memory Out of GPU Memory Algorithms Perform left-looking factorizations on sub-matrices that fit in the GPU memory (using existing algorithms) The rest of the matrix stays on the CPU Left-looking versions minimize writing on the CPU 1) Copy A2 to the GPU Untouched 2) Update A2 using A1 (a panel of A1 at a time) part of the 3) Factor the updated A2 using To be Factored factored matrix existing sub-matric sub... A1 on CPU matrix hybrid code A2 on 4) Copy factored A2 to the CPU GPU Trivially extended to multiGPUs: A2 is “larger” with 1-D block cyclic distribution, again reusing existing algorithms Hybrid Algorithms Two-sided factorizations (to bidiagonal, tridiagonal, and upper Hessenberg forms) for eigen- and singular-value problems Hybridization – Trailing matrix updates (Level 3 BLAS) are done on the GPU (similar to the one-sided factorizations) – Panels (Level 2 BLAS) are hybrid – operations with memory footprint restricted to the panel are done on CPU – The time consuming matrix-vector products involving the entire trailing matrix are done on the GPU MultiGPUs Two-Sided Factorizations Need HP multiGPU Level 2 BLAS (e.g., 50% of flops in the tridiagonal reduction) Performance of DSYMV on multi M2090s T. Dong, J. Dongarra, S. Tomov, I. Yamazaki, T. Schulthess, and R. Solca, Symmetric dense matrix-vector multiplication on multiple GPUs and its application to symmetric dense and sparse eigenvalue problems, ICL Technical report, 03/2012. Hybrid Two-Sided Factorizations MAGMA Tridiagonalization in DP 50 % of the flops are in SYMV Memory bound, i.e. does not scale well on multicore CPUs Use the GPU’s high memory bandwidth and optimized SYMV 8 x speedup over 12 Intel cores (X5660 @2.8 GHz) Keeneland system, using one node 3 NVIDIA GPUs ([email protected] 1.1 GHz, 5.4 GB) 2 x 6 Intel Cores (X5660 @ 2.8 GHz, 23 GB) T. Dong, J. Dongarra, S. Tomov, I. Yamazaki, T. Schulthess, and R. Solca, Symmetric dense matrix-vector multiplication on multiple GPUs and its application to symmetric dense and sparse eigenvalue problems, ICL Technical report, 03/2012. Further GPU Kernel Optimizations Fermi 2070 A. Abdelfattah, J. Dongarra, D. Keyes and H. Ltaief, Optimizing Memory-Bound Numerical Kernels on GPU Hardware Accelerators, VECPAR, Japan, 2012. Further GPU Kernel Optimizations Fermi 2070 Further GPU Kernel Optimizations Fermi 2070 From Static to Dynamic Scheduling… Static may stall in situations where work is available Hand tuned optimizations Hardware heterogeneity Kernel heterogeneity Separation of concerns Dynamic Runtime System Punch Lines Productivity! Productivity! Productivity! Oh… Did I say Productivity? Block Algorithms Panel-Update Sequence Transformations are blocked/accumulated within the Panel (Level 2 BLAS) Transformations applied at once on the trailing submatrix (Level 3 BLAS) Parallelism hidden inside the BLAS Fork-join Model Block QR Factorization Fork-Join Paradigm Leveraging Block Algorithms… Column-major data layout Tile data layout Lessons Learnt from PLASMA PLASMA: Parallel Linear Algebra for Scalable Multi-core Architectures: http://icl.cs.utk.edu/plasma/ Tile Algorithms on homogeneous x86 cores Parallelism is brought to the fore May require the redesign of linear algebra algorithms Tile data layout translation Remove unnecessary synchronization points between PanelUpdate sequences DAG execution where nodes represent tasks and edges define dependencies between them Dynamic runtime system environment QUARK Example: Tile QR Factorization First panel factorization and corresponding updates DAG for a 4x4 tiles matrix Let’s go crazy! DAG of 20x20 tile QR Factoriza Dynamic Scheduling Conceptually similar to out-of-order processor scheduling because it has: • Dynamic runtime DAG scheduler • Out-of-order execution flow of fine-grained tasks • Task scheduling as soon as dependencies are satisfied • Producer-Consumer Data Flow Programming Model: five decades old concept • Think "how things connect" rather than "how things happen” • Assembly line • Inherently parallel Matrices Over Runtime Systems at Exascale (MORSE) Mission statement: "Design dense and sparse linear algebra methods that achieve the fastest possible time to an accurate solution on large-scale Hybrid systems” Runtime challenges due to the ever growing hardware complexity Algorithmic challenges to exploit the hardware capabilities at most Integrated into MAGMA software stack MAGMA-MORSE: x86 + MultiGPUs Lessons Learned from PLASMA! CUDA-based hybrid systems New high performance numerical kernels StarPU Runtime System (Augonnet et. Al, INRIA, Bordeaux) Both: x86 and GPUs => Hybrid Computations Similar to LAPACK in functionality Achieving High Level of Productivity From Sequential Nested-Loop Code to Parallel Execution for (k = 0; k < min(MT, NT); k++){ zgeqrt(A[k;k], ...); for (n = k+1; n < NT; n++) zunmqr(A[k;k], A[k;n], ...); for (m = k+1; m < MT; m++){ ztsqrt(A[k;k],,A[m;k], ...); for (n = k+1; n < NT; n++) ztsmqr(A[m;k], A[k;n], A[m;n], ...); } } Achieving High Level of Productivity From Sequential Nested-Loop Code to Parallel Execution for (k = 0; k < min(MT, NT); k++){ starpu_Insert_Task(&cl_zgeqrt, k , k, ...); for (n = k+1; n < NT; n++) starpu_Insert_Task(&cl_zunmqr, k, n, ...); for (m = k+1; m < MT; m++){ starpu_Insert_Task(&cl_ztsqrt, m, k, ...); for (n = k+1; n < NT; n++) starpu_Insert_Task(&cl_ztsmqr, m, n, k, ...); } } Hybrid Architecture Targeted ⇒ PCI Interconnect 16X 64Gb/s, very thin pipe! ⇒ Fermi C2050 448 cuda cores 515 Gflop/s Performance Charts Cholesky QR LU Symmetric Matrix Inversion Cholesky Performance 8 Intel x86 cores + 3 Tesla GPUs E. Agullo, C. Augonnet, J. Dongarra, H. Ltaief, S. Thibault, R. Namyst, S. Thibault, S. Tomov, Software for GPUs, GPU Computing GEMs, vol.2, 2011. QR Performance 8 AMD x86 cores + 4 Tesla GPUs E. Agullo, C. Augonnet, J. Dongarra, M. Faverge, H. Ltaief, S. Thibault, S. Tomov, IEEE International Parallel and Distributed Processing Symposium, 2011. QR Performance +~200 Gflop/s but 12 cores = ~150 Gflop/s 8 AMD x86 cores + 4 Tesla GPUs E. Agullo, C. Augonnet, J. Dongarra, M. Faverge, H. Ltaief, S. Thibault, S. Tomov, IEEE International Parallel and Distributed Processing Symposium, 2011. Performance Breakdown Kernel CPU GPU Speedup sgeqrt 9 Gflops 60 Gflops ~6 stsqrt 12 Gflops 67 Gflops ~6 sormqr 8.5 Gflops 227 Gflops ~ 27 stsmqr 10 Gflops 285 Gflops Task distribution observed on StarPU: ~ 27 • sgeqrt : 20% of tasks on GPUs • stsmqr : 92.5% of tasks on GPUs Taking advantage of heterogeneity ! • Only do what you are good for • Don’t do what you are not good for LU Performance 8 Intel x86 cores + 3 Tesla GPUs E. Agullo, C. Augonnet, J. Dongarra, M. Faverge, J. Langou, H. Ltaief and S. Tomov, ACS/IEEE International Conference on Computer Systems and Applications (best paper award), 2011 Symmetric Matrix Inversion A−1, Seriously??? YES! Critical component of the variance-covariance matrix computation in statistics Three steps: • Cholesky factorization (DPOTRF) • Inverting the Cholesky factor (DTRTRI) • Calculating the product of the inverted Cholesky factor with its transpose (DLAUUM) Built on previous work from E. Agullo, H. Bouwmeester, J. Dongarra, J. Kurzak, J. Langou, and L. Rosenberg Scheduling Algorithms as DAGs 44 A−1 Performance 8 Intel x86 cores + 2 Fermi GPUs H. Ibeid, D. Kaushik, D. Keyes and H. Ltaief, HIPC'11, India Summary and Future Directions Two methodologies for solving challenging DLA Static Scheduling (performance) Dynamic Scheduling (productivity) LAPACK compliant API Source codes freely available in MAGMA What’s next? • Extended numerical functionality • Distributed Memory Systems Colloborators / Support MAGMA [Matrix Algebra on GPU and Multicore Architectures] team http://icl.cs.utk.edu/magma/ PLASMA [Parallel Linear Algebra for Scalable Multicore Architectures] team http://icl.cs.utk.edu/plasma Collaborating partners University of Tennessee, Knoxville University of California, Berkeley University of Colorado, Denver INRIA, France KAUST, Saudi Arabia Questions?