Linear Algebra Libraries for High- Performance Computing: Scientific - - PowerPoint PPT Presentation

▶
linear algebra libraries for high performance computing
SMART_READER_LITE
LIVE PREVIEW

Linear Algebra Libraries for High- Performance Computing: Scientific - - PowerPoint PPT Presentation

Feb 26, 2023 107 likes •877 views

Linear Algebra Libraries for High- Performance Computing: Scientific Computing with Multicore and Accelerators Presenters Prof. Jack Dongarra (8:30 10:00) University of Tennessee & Oak Ridge National Lab Dr. Jakub Kurzak (10:30

slide-1
SLIDE 1

1

Linear Algebra Libraries for High- Performance Computing: Scientific Computing with Multicore and Accelerators

Presenters

  • Prof. Jack Dongarra (8:30 – 10:00)

University of Tennessee & Oak Ridge National Lab

  • Dr. Jakub Kurzak (10:30 – 12:00)

University of Tennessee

  • Prof. James Demmel (1:30 – 3:00)

University of California Berkeley

  • Dr. Michael Heroux (3:30 – 5:00)

Sandia National Laboratory

slide-2
SLIDE 2

Overview ¡of ¡Dense ¡Numerical ¡Linear ¡ Algebra ¡Libraries ¡

  • BLAS: ¡kernel ¡for ¡dense ¡linear ¡algebra ¡
  • LAPACK: ¡sequen;al ¡dense ¡linear ¡algebra ¡
  • ScaLAPACK: ¡parallel ¡distributed ¡dense ¡linear ¡

algebra ¡

2

Linear ¡Algebra ¡PACKage ¡

slide-3
SLIDE 3

3

What do you mean by performance?

¨ What is a xflop/s? Ø xflop/s is a rate of execution, some number of floating point operations per second.

Ø Whenever this term is used it will refer to 64 bit floating point

  • perations and the operations will be either addition or

multiplication.

Ø Tflop/s refers to trillions (1012) of floating point operations per second and Ø Pflop/s refers to 1015 floating point operations per second. ¨ What is the theoretical peak performance? Ø The theoretical peak is based not on an actual performance from a benchmark run, but on a paper computation to determine the theoretical peak rate of execution of floating point operations for the machine. Ø The theoretical peak performance is determined by counting the number of floating-point additions and multiplications (in full precision) that can be completed during a period of time, usually the cycle time of the machine. Ø For example, an Intel Xeon 5570 quad core at 2.93 GHz can complete 4 floating point operations per cycle or a theoretical peak performance of 11.72 GFlop/s per core or 46.88 Gflop/s for the socket.

slide-4
SLIDE 4

4

What Is LINPACK?

¨

LINPACK is a package of mathematical software for solving problems in linear algebra, mainly dense linear systems of linear equations.

¨

LINPACK: “LINear algebra PACKage”

Ø Written in Fortran 66 ¨

The project had its origins in 1974

¨

The project had four primary contributors: myself when I was at Argonne National Lab, Jim Bunch from the University of California-San Diego, Cleve Moler who was at New Mexico at that time, and Pete Stewart from the University of Maryland.

¨

LINPACK as a software package has been largely superseded by LAPACK, which has been designed to run efficiently on shared- memory, vector supercomputers.

slide-5
SLIDE 5

5

Computing in 1974

¨ High Performance Computers:

Ø IBM 370/195, CDC 7600, Univac 1110, DEC PDP-10, Honeywell 6030

¨ Fortran 66 ¨ Trying to achieve software portability ¨ Run efficiently ¨ BLAS (Level 1)

Ø Vector operations

¨ Software released in 1979

Ø About the time of the Cray 1

slide-6
SLIDE 6

6

LINPACK Benchmark?

¨ The Linpack Benchmark is a measure of a

computer’s floating-point rate of execution.

Ø It is determined by running a computer program that solves a dense system of linear equations.

¨ Over the years the characteristics of the

benchmark has changed a bit.

Ø In fact, there are three benchmarks included in the Linpack Benchmark report.

¨ LINPACK Benchmark

Ø Dense linear system solve with LU factorization using partial pivoting Ø Operation count is: 2/3 n3 + O(n2) Ø Benchmark Measure: MFlop/s Ø Original benchmark measures the execution rate for a Fortran program on a matrix of size 100x100.

slide-7
SLIDE 7

7

Accidental Benchmarker

¨ Appendix B of the Linpack Users’ Guide Ø Designed to help users extrapolate execution time for Linpack software package ¨ First benchmark report from 1977; Ø Cray 1 to DEC PDP-10

slide-8
SLIDE 8

8

High Performance Linpack (HPL)

Benchmark Matrix Optimizations Parallel Name dimension allowed Processing Linpack 100 100 compiler –a Linpack 1000b 1000 hand, code replacement –c Linpack Parallel 1000 hand, code replacement Yes HPLinpackd arbitrary hand, code replacement Yes

a Compiler parallelization possible. b Also known as TPP (Toward Peak Performance) or Best Effort c Multiprocessor implementations allowed. d Highly-Parallel LINPACK Benchmark is also known as NxN Linpack

Benchmark or High Parallel Computing (HPC).

slide-9
SLIDE 9

A brief history of (Dense) Linear Algebra software

¨ But the BLAS-1 weren’t enough

Ø Consider AXPY ( y = α·x + y ): 2n flops on 3n read/writes Ø Computational intensity = (2n)/(3n) = 2/3 Ø Too low to run near peak speed (read/write dominates)

¨ So the BLAS-2 were developed (1984-1986)

Ø Standard library of 25 operations (mostly) on matrix/vector pairs

Ø “GEMV”: y = α·A·x + β·x, “GER”: A = A + α·x·yT,

x = T-1·x

Ø Up to 4 versions of each (S/D/C/Z), 66 routines, 18K LOC

Ø Why BLAS 2 ? They do O(n2) ops on O(n2) data Ø So computational intensity still just ~(2n2)/(n2) = 2

Ø OK for vector machines, but not for machine with caches

9

slide-10
SLIDE 10

A brief history of (Dense) Linear Algebra software

¨ The next step: BLAS-3 (1987-1988)

Ø Standard library of 9 operations (mostly) on matrix/matrix pairs

Ø “GEMM”: C = α·A·B + β·C, C = α·A·AT + β·C, B = T-1·B Ø Up to 4 versions of each (S/D/C/Z), 30 routines, 10K LOC

Ø Why BLAS 3 ? They do O(n3) ops on O(n2) data Ø So computational intensity (2n3)/(4n2) = n/2 – big at last!

Ø Good for machines with caches, other mem. hierarchy levels

¨ How much BLAS1/2/3 code so far (all at

www.netlib.org/blas)

Ø Source: 142 routines, 31K LOC, Testing: 28K LOC Ø Reference (unoptimized) implementation only

Ø Ex: 3 nested loops for GEMM

10

slide-11
SLIDE 11

Memory Hierarchy

¨ By taking advantage of the principle of locality:

Ø Present the user with as much memory as is available in the cheapest technology. Ø Provide access at the speed offered by the fastest technology.

Control Datapath Secondary Storage (Disk) Processor Registers Main Memory (DRAM) Level 2 and 3 Cache (SRAM) On-Chip Cache 1s 10,000,000s (10s ms) 100,000 s (.1s ms) Speed (ns): 10s 100s 100s Gs Size (bytes): Ks Ms Tertiary Storage (Disk/Tape) 10,000,000,000s (10s sec) 10,000,000 s (10s ms) Ts Distributed Memory Remote Cluster Memory

slide-12
SLIDE 12

12

Why Higher Level BLAS?

¨ Can only do arithmetic on data at the top of

the hierarchy

¨ Higher level BLAS lets us do this

Registers L 1 Cache L 2 Cache Local Memory Remote Memory Secondary Memory

slide-13
SLIDE 13

13

Level 1, 2 and 3 BLAS

¨ Level 1 BLAS

Vector-Vector

  • perations

¨ Level 2 BLAS

Matrix-Vector

  • perations

¨ Level 3 BLAS

Matrix-Matrix

  • perations

+ * * + *

slide-14
SLIDE 14

Level ¡1, ¡2 ¡and ¡3 ¡BLAS ¡

Before ¡(2007) ¡

3. 4 G H z EM 64T Xeon M KL8. 1 Peak: 6. 8 G f l

  • p/

s gcc

  • f
  • m i

t

  • f

r am e- poi nt er

  • f

unr

  • l

l

  • al

l

  • l
  • ops
  • O 3

1 2 3 4 5 6 7 1 19 37 55 73 91 109 127 145 163 181 199 217 235 253 271 289 307 325 343 361 379 397 415 433 451 469 487 505 O r der G f l

  • p/

s daxpy dgem v dgem m

slide-15
SLIDE 15

Level ¡1, ¡2 ¡and ¡3 ¡BLAS ¡

Now ¡(2011) ¡ ¡

0 ¡ 2 ¡ 4 ¡ 6 ¡ 8 ¡ 10 ¡ 12 ¡ 200 ¡ 400 ¡ 600 ¡ 800 ¡ 1000 ¡ 2000 ¡ 3000 ¡ 4000 ¡ 5000 ¡

GFLOPS ¡ Matrix ¡Size ¡

Level ¡3 ¡BLAS: ¡DGEMM ¡ Level ¡2 ¡BLAS: ¡DGEMV ¡ Level ¡1 ¡BLAS: ¡DAXPY ¡

AMD Opteron 8439 SE Processor (6 cores total @ 2.8Ghz) Using 1 core 11.2 Gflop/s theoretical peak

slide-16
SLIDE 16

02/25/2009 CS267 Lecture 8 16

slide-17
SLIDE 17

A brief history of (Dense) Linear Algebra software

¨ LAPACK – “Linear Algebra PACKage” - uses BLAS-3 (1989 – now)

Ø Ex: Obvious way to express Gaussian Elimination (GE) is adding multiples of one row to other rows – BLAS-1

Ø How do we reorganize GE to use BLAS-3 ? (details later)

Ø Contents of LAPACK (summary)

Ø Algorithms we can turn into (nearly) 100% BLAS 3 Ø Linear Systems: solve Ax=b for x Ø Least Squares: choose x to minimize ||Ax - b||2 Ø Algorithms that are only 50% BLAS 3 (so far) Ø “Eigenproblems”: Find λ and x where Ax = λ x Ø Singular Value Decomposition (SVD): (ATA)x=σ2x Ø Generalized problems (eg Ax = λ Bx) Ø Error bounds for everything Ø Lots of variants depending on A’s structure (banded, A=AT, etc)

Ø How much code? (Release 3.4, Nov 2011) (www.netlib.org/ lapack)

Ø Source: 1674 routines, 490K LOC, Testing: 448K LOC

17

slide-18
SLIDE 18

A brief history of (Dense) Linear Algebra software

¨ Is LAPACK parallel?

Ø Only if the BLAS are parallel (possible in shared memory)

¨ ScaLAPACK – “Scalable LAPACK” (1995 – now)

Ø For distributed memory – uses MPI Ø More complex data structures, algorithms than LAPACK

Ø Only (small) subset of LAPACK’s functionality available

Ø All at www.netlib.org/scalapack

18

slide-19
SLIDE 19

LAPACK

¨

http://www.netlib.org/lapack/

¨

LAPACK (Linear Algebra Package) provides routines for Ø solving systems of simultaneous linear equations, Ø least-squares solutions of linear systems of equations, Ø eigenvalue problems, Ø and singular value problems.

¨

LAPACK relies on BLAS

¨

The associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are related computations such as reordering of the Schur factorizations and estimating condition numbers.

¨

Dense and banded matrices are handled, but not general sparse matrices. In all areas, similar functionality is provided for real and complex matrices, in both single and double precision.

LAPACK is in FORTRAN Column Major LAPACK is SEQUENTIAL LAPACK is a REFERENCE implementation

19

slide-20
SLIDE 20

A ¡new ¡genera;on ¡of ¡algorithms? ¡

Algorithms ¡follow ¡hardware ¡evoluIon ¡along ¡Ime. ¡ LINPACK ¡(80’s) ¡ (Vector ¡opera;ons) ¡ ¡ ¡ ¡ ¡ Rely ¡on ¡ ¡ ¡ ¡-­‑ ¡Level-­‑1 ¡BLAS ¡opera;ons ¡ LAPACK ¡(90’s) ¡ (Blocking, ¡cache ¡friendly) ¡ ¡ ¡ ¡ ¡ Rely ¡on ¡ ¡ ¡ ¡-­‑ ¡Level-­‑3 ¡BLAS ¡opera;ons ¡

20

slide-21
SLIDE 21

¡subrouIne ¡dgesv( ¡n, ¡nrhs, ¡A, ¡lda, ¡ipiv, ¡b, ¡ldb, ¡info ¡) ¡ input: ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

  • utput: ¡

¡ ¡ ¡ ¡ ¡ ¡ ¡

Example ¡with ¡GESV ¡

Solve ¡a ¡system ¡of ¡linear ¡equa;ons ¡using ¡a ¡LU ¡factoriza;on ¡

A ¡ B ¡

n n n nrhs i p ¡ i ¡ v n n n

L U X ¡

n nrhs i p ¡ i ¡ v n info info Solu;on ¡of ¡Ax=b ¡

21

slide-22
SLIDE 22

Func;onali;es ¡in ¡LAPACK ¡

Type ¡of ¡Problem ¡ Acronyms ¡ Linear ¡system ¡of ¡equa;ons ¡ ¡ SV ¡ Linear ¡least ¡squares ¡problems ¡ ¡ LLS ¡ Linear ¡equality-­‑constrained ¡least ¡squares ¡problem ¡ ¡ LSE ¡ General ¡linear ¡model ¡problem ¡ ¡ GLM ¡ Symmetric ¡eigenproblems ¡ ¡ SEP ¡ ¡ Nonsymmetric ¡eigenproblems ¡ ¡ NEP ¡ Singular ¡value ¡decomposi;on ¡ ¡ SVD ¡ Generalized ¡symmetric ¡definite ¡eigenproblems ¡ ¡ GSEP ¡ Generalized ¡nonsymmetric ¡eigenproblems ¡ ¡ GNEP ¡ Generalized ¡(or ¡quo;ent) ¡singular ¡value ¡decomposi;on ¡ ¡ GSVD ¡(QSVD) ¡ ¡

22

slide-23
SLIDE 23

LAPACK ¡Socware ¡

  • First ¡release ¡in ¡February ¡1992 ¡
  • Version ¡3.4.0 ¡released ¡in ¡November ¡2011 ¡
  • LICENSE: ¡Mod-­‑BSD, ¡freely-­‑available ¡socware ¡package ¡-­‑ ¡Thus, ¡it ¡can ¡be ¡included ¡in ¡

commercial ¡socware ¡packages ¡(and ¡has ¡been). ¡We ¡only ¡ask ¡that ¡proper ¡credit ¡be ¡ given ¡to ¡the ¡authors. ¡

  • Open ¡SVN ¡repository ¡ ¡
  • MulI-­‑OS ¡

– *nix, ¡Mac ¡OS/X, ¡Windows ¡

  • MulI-­‑build ¡support ¡(cmake) ¡

– ¡make, ¡xcode, ¡nmake, ¡VS ¡studio, ¡Eclipse, ¡etc.. ¡

  • LAPACKE: ¡Standard ¡C ¡language ¡APIs ¡for ¡LAPACK ¡(In ¡collabora;on ¡with ¡INTEL) ¡

– ¡2 ¡layers ¡of ¡interface ¡

  • High-­‑Level ¡Interface ¡: ¡Workspace ¡alloca;on ¡and ¡NAN ¡Check ¡
  • Low-­‑Level ¡Interface ¡
  • Prebuilt ¡Libraries ¡for ¡Windows ¡
  • Extensive ¡test ¡suite ¡
  • Forum ¡and ¡User ¡support: ¡hip://icl.cs.utk.edu/lapack-­‑forum/ ¡

23

slide-24
SLIDE 24

Latest ¡Algorithms ¡

Since ¡release ¡3.0 ¡of ¡LAPACK ¡ – Hessenberg ¡QR ¡algorithm ¡with ¡the ¡small ¡bulge ¡mul;-­‑shic ¡QR ¡algorithm ¡together ¡with ¡ aggressive ¡early ¡defla;on. ¡[2003 ¡SIAM ¡SIAG ¡LA ¡Prize ¡winning ¡algorithm ¡of ¡Braman, ¡Byers ¡and ¡ Mathias] ¡ – Improvements ¡of ¡the ¡Hessenberg ¡reduc;on ¡subrou;nes. ¡[G. ¡Quintana-­‑Orn ¡and ¡van ¡de ¡Geijn] ¡ – New ¡MRRR ¡eigenvalue ¡algorithms ¡[2006 ¡SIAM ¡SIAG ¡LA ¡Prize ¡winning ¡algorithm ¡of ¡Dhillon ¡and ¡ Parlei] ¡ – New ¡par;al ¡column ¡norm ¡upda;ng ¡strategy ¡for ¡QR ¡factoriza;on ¡with ¡column ¡pivo;ng. ¡[Drmač ¡ and ¡Bujanovic] ¡ – Mixed ¡Precision ¡Itera;ve ¡Refinement ¡for ¡exploi;ng ¡fast ¡single ¡precision ¡hardware ¡for ¡GE, ¡PO ¡ [Langou’s] ¡ – Variants ¡of ¡various ¡factoriza;on ¡(LU, ¡QR, ¡Chol) ¡[Du] ¡ – RFP ¡(Rectangular ¡Full ¡Packed) ¡ ¡format ¡[Gustavson, ¡Langou] ¡ – XBLAS ¡and ¡Extra ¡precise ¡itera;ve ¡refinement ¡for ¡GESV ¡[Demmel ¡et ¡al.]. ¡ – New ¡fast ¡and ¡accurate ¡Jacobi ¡SVD ¡[2009 ¡SIAM ¡SIAG ¡LA ¡Prize, ¡Drmač ¡and ¡Veselić] ¡ – Pivoted ¡Cholesky ¡[Lucas] ¡ – Beier ¡mul;shic ¡Hessenberg ¡QR ¡algorithm ¡with ¡early ¡aggressive ¡defla;on ¡[Byers] ¡ – Complete ¡CS ¡decomposi;on ¡[Suion] ¡ – Level-­‑3 ¡BLAS ¡symmetric ¡indefinite ¡solve ¡and ¡symmetric ¡indefinite ¡inversion ¡[Langou’s] ¡ – Since ¡LAPACK ¡3.3, ¡all ¡rou;nes ¡in ¡are ¡now ¡thread-­‑safe ¡

¡

24

3.1 ¡ 3.2 ¡ 3.3 ¡

slide-25
SLIDE 25

LAPACK ¡3.4.0 ¡

  • xGEQRT: ¡QR ¡factoriza;on ¡(improved ¡interface). ¡

¡Contribu)on ¡by ¡Rodney ¡James, ¡UC ¡Denver. ¡ ¡xGEQRT ¡is ¡analogous ¡to ¡xGEQRF ¡with ¡a ¡modified ¡interface ¡which ¡enables ¡beier ¡ ¡performance ¡when ¡the ¡blocked ¡reflectors ¡need ¡to ¡be ¡reused. ¡The ¡companion ¡subrou;nes ¡ ¡xGEMQRT ¡apply ¡the ¡reflectors. ¡

  • xGEQRT3: ¡recursive ¡QR ¡factoriza;on. ¡

¡Contribu)on ¡by ¡Rodney ¡James, ¡UC ¡Denver. ¡ ¡The ¡recursive ¡QR ¡factoriza;on ¡enable ¡cache-­‑oblivious ¡and ¡enable ¡high ¡performance ¡on ¡tall ¡ ¡and ¡skinny ¡matrices. ¡

  • xTPQRT: ¡Communica;on-­‑Avoiding ¡QR ¡sequen;al ¡kernels. ¡

¡Contribu)on ¡by ¡Rodney ¡James, ¡UC ¡Denver. ¡ ¡These ¡subrou;nes ¡are ¡useful ¡for ¡upda;ng ¡a ¡QR ¡factoriza;on ¡and ¡are ¡used ¡in ¡sequen;al ¡and ¡ ¡parallel ¡Communica;on ¡Avoiding ¡QR. ¡These ¡subrou;nes ¡support ¡the ¡general ¡case ¡Triangle ¡ ¡on ¡top ¡of ¡Pentagone ¡which ¡includes ¡as ¡special ¡cases ¡the ¡so-­‑called ¡triangle ¡on ¡top ¡of ¡triangle ¡ ¡and ¡triangle ¡on ¡top ¡of ¡square. ¡This ¡is ¡the ¡right-­‑looking ¡version ¡of ¡the ¡subrou;nes ¡and ¡the ¡ ¡rou;ne ¡is ¡blocked.The ¡T ¡matrices ¡and ¡the ¡block ¡size ¡are ¡part ¡of ¡the ¡interface. ¡The ¡ ¡companion ¡subrou;nes ¡xTPMQRT ¡apply ¡the ¡reflectors. ¡

  • xSYEVK: ¡LDLT ¡with ¡rook ¡pivo;ng ¡and ¡fast ¡Bunch-­‑Parlei ¡pivo;ng. ¡

¡Contribu)on ¡by ¡Craig ¡Lucas. ¡ ¡These ¡subrou;nes ¡enables ¡beier ¡stability ¡than ¡the ¡Bunch-­‑Kaufman ¡pivo;ng ¡scheme ¡ ¡(xSYEV) ¡currently ¡used ¡in ¡LAPACK. ¡The ¡computa;onal ¡;me ¡is ¡slightly ¡higher. ¡

25

slide-26
SLIDE 26

Resources ¡

Reference ¡Code: ¡ – Reference ¡code: ¡(current ¡version ¡3.3.1) ¡ hip://www.netlib.org/lapack/lapack.tgz ¡ – LAPACK ¡build ¡for ¡windows ¡(current ¡version ¡3.3.1) ¡ ¡ ¡hip://icl.cs.utk.edu/lapack-­‑for-­‑windows/lapack ¡ ¡ – LAPACKE: ¡Standard ¡C ¡language ¡APIs ¡for ¡LAPACK ¡(in ¡collabora;on ¡with ¡INTEL): ¡ hip://www.netlib.org/lapack/#_standard_c_language_apis_for_lapack ¡ – Remi’s ¡wrappers ¡(wrapper ¡for ¡Matlab ¡users): ¡ hip://icl.cs.utk.edu/~delmas/lapwrapmw.htm ¡ ¡ ¡ Vendor ¡Libraries: ¡

¡more ¡or ¡less ¡same ¡as ¡the ¡BLAS: ¡MKL, ¡ACML, ¡VECLIB, ¡ESSL, ¡etc… ¡(WARNING: ¡some ¡implementa;ons ¡are ¡just ¡a ¡subset ¡of ¡LAPACK) ¡ ¡

DocumentaIon: ¡ – LAPACK ¡Users’ ¡guide: ¡ hip://www.netlib.org/lapack/lug/ ¡ – LAPACK ¡Working ¡notes ¡(in ¡par;cular ¡LAWN ¡41) ¡ hip://www.netlib.org/lapack/lawns/downloads/ ¡ – LAPACK ¡release ¡notes ¡ hip://www.netlib.org/lapack/lapack-­‑3.1.0.changes ¡ – LAPACK ¡NAG ¡example ¡and ¡auxiliary ¡rou;nes ¡ hip://www.nag.com/lapack-­‑ex/lapack-­‑ex.html ¡ ¡ – CRC ¡Handbook ¡of ¡Linear ¡Algebra, ¡Leslie ¡Hogben ¡ed, ¡Packages ¡of ¡SubrouInes ¡for ¡Linear ¡Algebra, ¡Bai, ¡Demmel, ¡Dongarra, ¡Langou, ¡and ¡Wang, ¡ Sec;on ¡75: ¡pages ¡75-­‑1,75-­‑24, ¡CRC ¡Press, ¡2006. ¡ hip://www.netlib.org/netlib/utk/people/JackDongarra/PAPERS/CRC-­‑LAPACK-­‑2005.pdf ¡ Support: ¡ – LAPACK ¡forum: ¡(more ¡than ¡1000 ¡topics) ¡ hip://icl.cs.utk.edu/lapack-­‑forum/ ¡ – LAPACK ¡mailing-­‑list: ¡ lapack@cs.utk.edu ¡ – LAPACK ¡mailing-­‑list ¡archive: ¡ hip://icl.cs.utk.edu/lapack-­‑forum/archives/ ¡

¡

26

slide-27
SLIDE 27

Organizing Linear Algebra – in books

www.netlib.org/lapack www.netlib.org/scalapack www.cs.utk.edu/~dongarra/etemplates www.netlib.org/templates

slide-28
SLIDE 28

ParallelizaIon ¡of ¡LU ¡and ¡QR. ¡ Parallelize ¡the ¡update: ¡

  • ¡Easy ¡and ¡done ¡in ¡any ¡reasonable ¡socware. ¡
  • ¡This ¡is ¡the ¡2/3n3 ¡term ¡in ¡the ¡FLOPs ¡count. ¡
  • ¡Can ¡be ¡done ¡efficiently ¡with ¡LAPACK+mul;threaded ¡BLAS ¡
  • dgemm
  • lu( )

dgetf2 dtrsm (+ dswp) dgemm \ L U A(1) A(2) L U

slide-29
SLIDE 29

Overview ¡of ¡Dense ¡Numerical ¡Linear ¡ Algebra ¡Libraries ¡

  • BLAS: ¡kernel ¡for ¡dense ¡linear ¡algebra ¡
  • LAPACK: ¡sequen;al ¡dense ¡linear ¡algebra ¡
  • ScaLAPACK: ¡parallel ¡distributed ¡dense ¡linear ¡

algebra ¡ ¡ Scalable ¡Linear ¡Algebra ¡PACKage ¡

slide-30
SLIDE 30

30

ScaLAPACK

¨ Library of software dealing with dense &

banded routines

¨ Distributed Memory - Message Passing ¨ MIMD Computers and Networks of

Workstations

¨ Clusters of SMPs

slide-31
SLIDE 31

ScaLAPACK ¡

  • hip://www.netlib.org/scalapack/ ¡
  • ScaLAPACK ¡(Scalable ¡Linear ¡Algebra ¡Package) ¡provides ¡

rou;nes ¡for ¡ ¡

– solving ¡systems ¡of ¡simultaneous ¡linear ¡equa;ons, ¡ ¡ – least-­‑squares ¡solu;ons ¡of ¡linear ¡systems ¡of ¡equa;ons, ¡ ¡ – eigenvalue ¡problems, ¡ ¡ – and ¡singular ¡value ¡problems. ¡

  • Relies ¡on ¡LAPACK ¡/ ¡BLAS ¡and ¡BLACS ¡/ ¡MPI ¡
  • Includes ¡PBLAS ¡( ¡Parallel ¡BLAS ¡) ¡ ¡

ScaLAPACK ¡is ¡in ¡ FORTRAN ¡and ¡C ¡ ScaLAPACK ¡is ¡for ¡ PARALLEL ¡ DISTRIBUTED ¡ ScaLAPACK ¡is ¡a ¡ REFERENCE ¡ implementa;on ¡

slide-32
SLIDE 32

32

Programming Style

¨ SPMD Fortran 77 with object based design ¨ Built on various modules

Ø PBLAS Interprocessor communication Ø BLACS

Ø PVM, MPI, IBM SP, CRI T3, Intel, TMC Ø Provides right level of notation.

Ø BLAS

¨ LAPACK software expertise/quality

Ø Software approach Ø Numerical methods

slide-33
SLIDE 33

33

Overall Structure of Software

¨ Object based - Array descriptor

Ø Contains information required to establish mapping between a global array entry and its corresponding process and memory location. Ø Provides a flexible framework to easily specify additional data distributions or matrix types. Ø Currently dense, banded, & out-of-core

¨ Using the concept of context

slide-34
SLIDE 34

34

PBLAS

¨ Similar to the BLAS in functionality and

naming.

¨ Built on the BLAS and BLACS ¨ Provide global view of matrix

CALL DGEXXX ( M, N, A( IA, JA ), LDA,... ) CALL PDGEXXX( M, N, A, IA, JA, DESCA,... )

slide-35
SLIDE 35

35

ScaLAPACK Structure

ScaLAPACK BLAS LAPACK BLACS PVM/MPI/... PBLAS Global Local

slide-36
SLIDE 36

36

Choosing a Data Distribution

¨ Main issues are:

Ø Load balancing Ø Use of the Level 3 BLAS

slide-37
SLIDE 37

37

Possible Data Layouts

¨ 1D block and cyclic column distributions ¨ 1D block-cycle column and 2D block-cyclic

distribution

¨ 2D block-cyclic used in ScaLAPACK for dense

matrices

slide-38
SLIDE 38

[LAPACK] ¡subrouIne ¡dgesv( ¡n, ¡nrhs, ¡a(ia,ja), ¡lda, ¡ipiv, ¡b(ib,jb), ¡ldb, ¡info ¡) ¡ ¡ ¡ input: ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

  • utput: ¡

¡ ¡ ¡ ¡ ¡ ¡ ¡

From ¡LAPACK ¡to ¡ScaLAPACK ¡

A ¡ B ¡

n ¡ n ¡ n ¡ nrhs ¡ i ¡ p ¡ i ¡ v n ¡ n ¡ n ¡

L ¡ U ¡ X ¡

n ¡ nrhs ¡ i ¡ p ¡ i ¡ v n ¡ info ¡ info ¡ LAPACK ¡Data ¡layout ¡ LAPACK ¡Data ¡layout ¡

slide-39
SLIDE 39

[LAPACK] ¡subrouIne ¡dgesv( ¡n, ¡nrhs, ¡a(ia,ja), ¡lda, ¡ipiv, ¡b(ib,jb), ¡ldb, ¡info ¡) ¡ ¡ ¡ input: ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

  • utput: ¡

¡ ¡ ¡ ¡ ¡ ¡ ¡

From ¡LAPACK ¡to ¡ScaLAPACK ¡

n ¡ n ¡ info ¡

A11 ¡ A21 ¡ A31 ¡ A12 ¡ A22 ¡ A32 ¡ A13 ¡ A23 ¡ A33 ¡

n ¡

B11 ¡ B21 ¡ B31 ¡

nrhs ¡

ip1 ¡ ip2 ¡ ip3 ¡

n ¡ info ¡

L21 ¡ L31 ¡ U12 ¡ L32 ¡ U13 ¡ U23 ¡

n ¡

X11 ¡ X21 ¡ X31 ¡ ip1 ¡ ip2 ¡ ip3 ¡

nrhs ¡ L11 ¡ U11 ¡ L22 ¡ U22 ¡ L33 ¡ U33 ¡ n ¡ ScaLAPACK ¡Data ¡layout ¡ ScaLAPACK ¡Data ¡layout ¡

slide-40
SLIDE 40

[LAPACK] ¡subrouIne ¡dgesv( ¡n, ¡nrhs, ¡a(ia,ja), ¡lda, ¡ipiv, ¡b(ib,jb), ¡ldb, ¡info ¡) ¡ ¡ [ScaLAPACK] ¡subrouIne ¡pdgesv( ¡n, ¡nrhs, ¡a, ¡ia, ¡ja, ¡desca, ¡ipiv, ¡b, ¡ib, ¡jb, ¡descb, ¡info ¡) ¡ input: ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡

  • utput: ¡

¡ ¡ ¡ ¡ ¡ ¡ ¡

From ¡LAPACK ¡to ¡ScaLAPACK ¡

n ¡ n ¡ info ¡

A11 ¡ A21 ¡ A31 ¡ A12 ¡ A22 ¡ A32 ¡ A13 ¡ A23 ¡ A33 ¡

n ¡

B11 ¡ B21 ¡ B31 ¡

nrhs ¡

ip1 ¡ ip2 ¡ ip3 ¡

n ¡ info ¡

L21 ¡ L31 ¡ U12 ¡ L32 ¡ U13 ¡ U23 ¡

n ¡

X11 ¡ X21 ¡ X31 ¡ ip1 ¡ ip2 ¡ ip3 ¡

nrhs ¡ L11 ¡ U11 ¡ L22 ¡ U22 ¡ L33 ¡ U33 ¡ n ¡ ScaLAPACK ¡Data ¡layout ¡ ScaLAPACK ¡Data ¡layout ¡

slide-41
SLIDE 41

41

Distribution and Storage

¨ Matrix is block-partitioned & maps blocks ¨ Distributed 2-D block-cyclic scheme

5x5 matrix partitioned in 2x2 blocks 2x2 process grid point of view

¨ Routines available to distribute/redistribute

data.

A

11

A

12

A

13

A

14

A

15

A

21

A

22

A

23

A

24

A

25

A

31

A

32

A

33

A

34

A

35

A

41

A

42

A

43

A

44

A

45

A

51

A

52

A

53

A

54

A

55

A

11

A

12

A

15

A

13

A

14

A

21

A

22

A

25

A

23

A

24

A

51

A

52

A

55

A

53

A

54

A

31

A

32

A

35

A

33

A

34

A

41

A

42

A

45

A

43

A

44

slide-42
SLIDE 42

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view MB NB M N Matrix is MxN Process grid is PxQ, P=2, Q=3 Blocks are MBxNB P Q

slide-43
SLIDE 43

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-44
SLIDE 44

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view 2 4 1 3 5

slide-45
SLIDE 45

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view 2 4 1 3 5 2 4 1 3 5

slide-46
SLIDE 46

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 2 4 2 4 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-47
SLIDE 47

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 2 4 2 4 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-48
SLIDE 48

2D ¡Block ¡Cyclic ¡Layout ¡

2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 2 4 2 4 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 4 5 4 5 4 5 4 5 4 4 5 4 5 4 5 4 5 4 2 3 2 3 2 3 2 3 2 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-49
SLIDE 49

2D ¡Block ¡Cyclic ¡Layout ¡

3 5 1 3 5 1 3 5 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 3 5 2 4 1 3 5 2 4 1 3 5 2 4 2 4 2 4 3 5 2 4 3 5 2 4 3 5 2 4 3 5 2 4 5 4 5 4 5 4 5 4 5 4 5 4 5 4 5 4 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-50
SLIDE 50

2D ¡Block ¡Cyclic ¡Layout ¡

4 5 2 4 1 3 5 2 4 1 3 5 4 5 2 4 1 3 5 2 4 1 3 5 4 5 2 4 1 3 5 2 4 1 3 5 4 2 4 2 4 4 5 4 5 4 5 4 4 5 4 5 4 5 4 4 5 4 5 4 5 4 2 3 2 3 2 3 2 2 3 2 3 2 3 2 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-51
SLIDE 51

2D ¡Block ¡Cyclic ¡Layout ¡

1 3 5 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 1 3 5 2 4 2 4 5 4 5 4 5 4 5 4 5 4 5 4 3 2 3 2 3 2 3 2 3 2 3 2 1 1 1 1 1 1 Matrix point of view Processor point of view

slide-52
SLIDE 52

2D ¡Block ¡Cyclic ¡Layout ¡

Matrix point of view Processor point of view

slide-53
SLIDE 53

53

Parallelism in ScaLAPACK

¨ Level 3 BLAS block

  • perations

Ø All the reduction routines

¨ Pipelining

Ø QR Algorithm, Triangular Solvers, classic factorizations

¨ Redundant

computations

Ø Condition estimators

¨ Static work

assignment

Ø Bisection

¨ Task parallelism

Ø Sign function eigenvalue computations

¨ Divide and Conquer

Ø Tridiagonal and band solvers, symmetric eigenvalue problem and Sign function

¨ Cyclic reduction

Ø Reduced system in the band solver

¨ Data parallelism

Ø Sign function

slide-54
SLIDE 54

Func;onali;es ¡in ¡LAPACK ¡

Type ¡of ¡Problem ¡ Acronyms ¡ Linear ¡system ¡of ¡equa;ons ¡ ¡ SV ¡ Linear ¡least ¡squares ¡problems ¡ ¡ LLS ¡ Linear ¡equality-­‑constrained ¡least ¡squares ¡problem ¡ ¡ LSE ¡ General ¡linear ¡model ¡problem ¡ ¡ GLM ¡ Symmetric ¡eigenproblems ¡ ¡ SEP ¡ ¡ Nonsymmetric ¡eigenproblems ¡ ¡ NEP ¡ Singular ¡value ¡decomposi;on ¡ ¡ SVD ¡ Generalized ¡symmetric ¡definite ¡eigenproblems ¡ ¡ GSEP ¡ Generalized ¡nonsymmetric ¡eigenproblems ¡ ¡ GNEP ¡ Generalized ¡(or ¡quo;ent) ¡singular ¡value ¡decomposi;on ¡ ¡ GSVD ¡(QSVD) ¡ ¡

slide-55
SLIDE 55

Func;onnali;es ¡in ¡ScaLAPACK ¡

Type ¡of ¡Problem ¡ Acronyms ¡ Linear ¡system ¡of ¡equa;ons ¡ ¡ SV ¡ Linear ¡least ¡squares ¡problems ¡ ¡ LLS ¡ Linear ¡equality-­‑constrained ¡least ¡squares ¡problem ¡ ¡ LSE ¡ General ¡linear ¡model ¡problem ¡ ¡ GLM ¡ Symmetric ¡eigenproblems ¡ ¡ SEP ¡ ¡ Nonsymmetric ¡eigenproblems ¡ ¡ NEP ¡ Singular ¡value ¡decomposi;on ¡ ¡ SVD ¡ Generalized ¡symmetric ¡definite ¡eigenproblems ¡ ¡ GSEP ¡ Generalized ¡nonsymmetric ¡eigenproblems ¡ ¡ GNEP ¡ Generalized ¡(or ¡quo;ent) ¡singular ¡value ¡decomposi;on ¡ ¡ GSVD ¡(QSVD) ¡ ¡

slide-56
SLIDE 56

56

Major Changes to Software

  • Must rethink the design of our software

Ø Another disruptive technology

Ø Similar to what happened with cluster computing and message passing

Ø Rethink and rewrite the applications, algorithms, and software

  • Numerical libraries for example will change

Ø For example, both LAPACK and ScaLAPACK will undergo major changes to accommodate this

slide-57
SLIDE 57

A New Generation of Software:

Parallel Linear Algebra Software for Multicore Architectures (PLASMA)

Software/Algorithms follow hardware evolution in time LINPACK (70’s) (Vector operations) Rely on

  • Level-1 BLAS
  • perations

LAPACK (80’s) (Blocking, cache friendly) Rely on

  • Level-3 BLAS
  • perations

ScaLAPACK (90’s) (Distributed Memory) Rely on

  • PBLAS Mess Passing

PLASMA (00’s) New Algorithms (many-core friendly) Rely on

  • a DAG/scheduler
  • block data layout
  • some extra kernels

Those new algorithms

  • have a very low granularity, they scale very well (multicore, petascale

computing, … )

  • removes a lots of dependencies among the tasks, (multicore, distributed

computing)

slide-58
SLIDE 58

A New Generation of Software:

Parallel Linear Algebra Software for Multicore Architectures (PLASMA)

Software/Algorithms follow hardware evolution in time LINPACK (70’s) (Vector operations) Rely on

  • Level-1 BLAS
  • perations

LAPACK (80’s) (Blocking, cache friendly) Rely on

  • Level-3 BLAS
  • perations

ScaLAPACK (90’s) (Distributed Memory) Rely on

  • PBLAS Mess Passing

PLASMA (00’s) New Algorithms (many-core friendly) Rely on

  • a DAG/scheduler
  • block data layout
  • some extra kernels

Those new algorithms

  • have a very low granularity, they scale very well (multicore, petascale

computing, … )

  • removes a lots of dependencies among the tasks, (multicore, distributed

computing)

slide-59
SLIDE 59

A New Generation of Software:

Parallel Linear Algebra Software for Multicore Architectures (PLASMA)

Software/Algorithms follow hardware evolution in time LINPACK (70’s) (Vector operations) Rely on

  • Level-1 BLAS
  • perations

LAPACK (80’s) (Blocking, cache friendly) Rely on

  • Level-3 BLAS
  • perations

ScaLAPACK (90’s) (Distributed Memory) Rely on

  • PBLAS Mess Passing

PLASMA (00’s) New Algorithms (many-core friendly) Rely on

  • a DAG/scheduler
  • block data layout
  • some extra kernels

Those new algorithms

  • have a very low granularity, they scale very well (multicore, petascale

computing, … )

  • removes a lots of dependencies among the tasks, (multicore, distributed

computing)

slide-60
SLIDE 60

A New Generation of Software:

Parallel Linear Algebra Software for Multicore Architectures (PLASMA)

Software/Algorithms follow hardware evolution in time LINPACK (70’s) (Vector operations) Rely on

  • Level-1 BLAS
  • perations

LAPACK (80’s) (Blocking, cache friendly) Rely on

  • Level-3 BLAS
  • perations

ScaLAPACK (90’s) (Distributed Memory) Rely on

  • PBLAS Mess Passing

PLASMA (00’s) New Algorithms (many-core friendly) Rely on

  • a DAG/scheduler
  • block data layout
  • some extra kernels

Those new algorithms

  • have a very low granularity, they scale very well (multicore, petascale computing, … )
  • removes a lots of dependencies among the tasks, (multicore, distributed computing)
  • avoid latency (distributed computing, out-of-core)
  • rely on fast kernels

Those new algorithms need new kernels and rely on efficient scheduling algorithms.

slide-61
SLIDE 61

Moore’s Law is Alive and Well

61

1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06 1.E+07 1970 1975 1980 1985 1990 1995 2000 2005 2010

Transistors (in Thousands)

Data from Kunle Olukotun, Lance Hammond, Herb Sutter, Burton Smith, Chris Batten, and Krste Asanoviç

slide-62
SLIDE 62

But Clock Frequency Scaling Replaced by Scaling Cores / Chip

62

1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06 1.E+07 1970 1975 1980 1985 1990 1995 2000 2005 2010

Transistors (in Thousands) Frequency (MHz) Cores

Data from Kunle Olukotun, Lance Hammond, Herb Sutter, Burton Smith, Chris Batten, and Krste Asanoviç

15 Years of exponential growth ~2x year has ended

slide-63
SLIDE 63

Performance Has Also Slowed, Along with Power

63

1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06 1.E+07 1970 1975 1980 1985 1990 1995 2000 2005 2010

Transistors (in Thousands) Frequency (MHz) Power (W) Cores

Data from Kunle Olukotun, Lance Hammond, Herb Sutter, Burton Smith, Chris Batten, and Krste Asanoviç

Power is the root cause of all this

A hardware issue just became a software problem

slide-64
SLIDE 64

64

Power Cost of Frequency

  • Power ∝ Voltage2 x Frequency (V2F)
  • Frequency ∝ Voltage
  • Power ∝Frequency3
slide-65
SLIDE 65

65

Power Cost of Frequency

  • Power ∝ Voltage2 x Frequency (V2F)
  • Frequency ∝ Voltage
  • Power ∝Frequency3
slide-66
SLIDE 66

Moore’s Law Reinterpreted

¨ Number of cores per chip

doubles every 2 year, while clock speed decreases (not increases).

Ø Need to deal with systems with millions of concurrent threads

Ø Future generation will have billions of threads!

Ø Need to be able to easily replace inter-chip parallelism with intro-chip parallelism

¨ Number of threads of

execution doubles every 2 year

500000 1000000 1500000 2000000 2500000 3000000 3500000

Cores in the Top20 Systems

slide-67
SLIDE 67

Example of typical parallel machine

Chip/Socket Core Core Core Core

slide-68
SLIDE 68

Example of typical parallel machine

Node/Board Chip/Socket Chip/Socket Chip/Socket Core Core Core Core … Core GPU GPU GPU

slide-69
SLIDE 69

Example of typical parallel machine

Cabinet Node/Board Node/Board Node/Board Chip/Socket Chip/Socket Chip/Socket Core Core Core Core … Core Shared memory programming between processes on a board and a combination of shared memory and distributed memory programming between nodes and cabinets …

slide-70
SLIDE 70

Example of typical parallel machine

Switch Cabinet Cabinet Cabinet Node/Board Node/Board Node/Board Chip/Socket Chip/Socket Chip/Socket Core Core Core Core … … Core Combination of shared memory and distributed memory programming …

slide-71
SLIDE 71

November 2011: The TOP10

Rank Site Computer Country Cores Rmax [Pflops] % of Peak Power [MW] MFlops /Watt 1 RIKEN Advanced Inst for Comp Sci K computer Fujitsu SPARC64 VIIIfx + custom Japan 705,024 10.5 93 12.7 826 2

  • Nat. SuperComputer

Center in Tianjin Tianhe-1A, NUDT Intel + Nvidia GPU + custom China 186,368 2.57 55 4.04 636 3 DOE / OS Oak Ridge Nat Lab Jaguar, Cray AMD + custom USA 224,162 1.76 75 7.0 251 4

  • Nat. Supercomputer

Center in Shenzhen Nebulea, Dawning Intel + Nvidia GPU + IB China 120,640 1.27 43 2.58 493 5 GSIC Center, Tokyo Institute of Technology Tusbame 2.0, HP Intel + Nvidia GPU + IB Japan 73,278 1.19 52 1.40 850 6 DOE / NNSA LANL & SNL Cielo, Cray AMD + custom USA 142,272 1.11 81 3.98 279 7 NASA Ames Research Center/NAS Plelades SGI Altix ICE 8200EX/8400EX + IB USA 111,104 1.09 83 4.10 265 8 DOE / OS Lawrence Berkeley Nat Lab Hopper, Cray AMD + custom USA 153,408 1.054 82 2.91 362 9 Commissariat a l'Energie Atomique (CEA) Tera-10, Bull Intel + IB France 138,368 1.050 84 4.59 229 10 DOE / NNSA Los Alamos Nat Lab Roadrunner, IBM AMD + Cell GPU + IB USA 122,400 1.04 76 2.35 446

slide-72
SLIDE 72

Critical Issues for Peta and Exascale Algorithms

¨ Synchronization-reducing algorithms

Ø Break Fork-Join model

¨ Communication-reducing algorithms

Ø Use methods which have lower bound on communication

¨ Mixed precision methods

Ø 2x speed of ops and 2x speed for data movement

¨ Autotuning

Ø Today’s machines are too complicated, build “smarts” into software have experiment to optimize.

¨ Fault resilient algorithms

Ø Implement algorithms that can recover from failures/bit flips

¨ Reproducibility of results

Ø Today we can’t guarantee this. We understand the issues, but some of our “colleagues” have a hard time with this.

slide-73
SLIDE 73

73

Major Changes to Software

  • Must rethink the design of our software

Ø Another disruptive technology

Ø Similar to what happened with cluster computing and message passing

Ø Rethink and rewrite the applications, algorithms, and software

  • Numerical libraries for example will change

Ø For example, both LAPACK and ScaLAPACK will undergo major changes to accommodate this

slide-74
SLIDE 74

LAPACK LU - Intel64 - 16 cores

20 40 60 80 100 120 140 2000 4000 6000 8000 10000 12000 14000

Gflop/s Matrix size

DGETRF - Intel64 Xeon quad-socket quad-core (16 cores) - th. peak 153.6 Gflop/s

DGEMM LAPACK

slide-75
SLIDE 75

ParallelizaIon ¡of ¡QR ¡FactorizaIon ¡ Parallelize ¡the ¡update: ¡

  • ¡Easy ¡and ¡done ¡in ¡any ¡reasonable ¡socware. ¡
  • ¡This ¡is ¡the ¡2/3n3 ¡term ¡in ¡the ¡FLOPs ¡count. ¡
  • ¡Can ¡be ¡done ¡efficiently ¡with ¡LAPACK+mul;threaded ¡BLAS ¡
  • dgemm

75

  • qr( )

dgeqf2 + dlarft dlarfb V R A(1) A(2) V R

Update of the remaining submatrix

Panel factorization

Fork - Join parallelism Bulk Sync Processing

slide-76
SLIDE 76

Ø Objectives

Ø High utilization of each core Ø Scaling to large number of cores Ø Shared or distributed memory

Ø Methodology

Ø Dynamic DAG scheduling Ø Split phases task generation and execution Ø Explicit parallelism/Implicit communication Ø Fine granularity / block data layout

Ø Arbitrary DAG with dynamic scheduling

76

Cholesky 4 x 4 Fork-join parallelism

PLASMA: Parallel Linear Algebra s/w for Multicore Architectures

DAG scheduled parallelism Time