CS 267 Applications of Parallel Processors Lecture 13: Parallel Matrix Multiply

1. CS 267 Applications of Parallel ProcessorsLecture 13: Parallel Matrix Multiply Kathy Yelick http://www.cs.berkeley.edu/~dmartin/cs267

2. Outline - Recap - Sources of large dense linear systems - BLAS for linear algebra - Parallel Matrix multiply

3. Model overview • Work-depth • PRAM • Latency/Bandwidth model • a is the 1-time cost per message (latency) • b is the per byte cost of communication • Use this today • LogP model • correction: gap should be greater than overhead • more on this with parallel sorting • Topology-specific models

4. Dense Linear Algebra in Electromagentics

5. Computational Electromagnetics • developed during 1980s, driven by defense applications • determine the RCS (radar cross section) of airplane • reduce signature of plane (stealth technology) • other applications are antenna design, medical equipment • two fundamental numerical approaches: • MOM methods of moments ( frequency domain), and • finite differences (time domain)

6. Computational Electromagnetics - discretize surface into triangular facets using standard modeling tools - amplitude of currents on surface are unknowns - integral equation is discretized into a set of linear equations image: NW Univ. Comp. Electromagnetics Laboratory http://nueml.ece.nwu.edu/

7. Computational Electromagnetics (MOM) After discretization the integral equation has the formZ J = V where Z is the impedance matrix, J is the unknown vector of amplitudes, and V is the excitation vector. (see Cwik, Patterson, and Scott, Electromagnetic Scattering on the Intel Touchstone Delta, IEEE Supercomputing ‘92, pp 538 - 542)

8. Computational Electromagnetics (MOM) The main steps in the solution process are A) computing the matrix elements B) factoring the dense matrix C) solving for one or more excitations (RHS) D) computing the fields scattered from the object

9. Analysis of MOM for Parallel Implementation Task Work Parallelism Parallel Speed Fill O(n**2) embarrassing low Factor O(n**3) moderately diff. very high Solve O(n**2) moderately diff. high Field Calc. O(n) embarrassing high For most scientific applications the biggest gain in performance is from parallelism within each task.

10. Results for Parallel Implementation on Delta Task Time (hours) Fill 9.20 Factor 8.25 Solve 2.17 Field Calc. 0.12 The problem solved was for a matrix of size 48,672. (The world record in 1991.)

11. Current Records for Solving Dense Systems Year System Size Machine 1950's O(100) 1991 55,296 CM-2 1992 75,264 Intel 1993 75,264 Intel 1994 76,800 CM-5 1995 128,600 Paragon XP 1996 215,000 ASCI Red (Tflop) source: Alan Edelman http://www-math.mit.edu/~edelman/records.html

12. Sources for large dense linear systems - Not many basic factorizations outside CEM - Large dense eigen problems used in chemistry - Alternatives often debated Choice for algorithms in existing codes are not the result of careful planning and design. - Reflect the start-of-the-art at the time, - May be purely coincidental.

13. Solving Large Dense Linear Systems Gaussian elimination to solve Ax=b where A is a dense matrix • Add multiples of each row to subsequent rows in order to create zeros below the diagonal • End up with an upper triangular matrix U. • Solve a linear system with U by substitution, starting with the last variable. Solving these systems uses basic vector and matrix operations called BLAS. see Demmel http://HTTP.CS.Berkeley.EDU/~demmel/cs267/lecture12/lecture12.html

14. Parallel Matrix Multiply

15. Parallel Matrix Multiply • Computing C=C+A*B • Using basic algorithm: 2*n3 Flops • Variables are: • Data layout • Topology of machine • Scheduling communication • Use of performance models for algorithm design

16. p0 p1 p2 p3 p4 p5 p6 p7 1D Layout • Assume matrices are nxn and n is divisible by p • A(i) refers to the n by n/p block column that processor i owns (similiarly for B(i) and C(i)) • B(i,j) is the n/p by n/p sublock of B(i) • in rows j*n/p through (j+1)*n/p

17. Matrix Multiply: 1D Layout on Bus • Algorithm uses the formula C(i) = C(i) + A*B(i) = C(i) + S A(j)*B(j,i) • First consider a bus-connected machine without broadcast: only one pair of processors can communicate at a time • Second consider a bus-connected machine with broadcast: may send from one to many in single step j <p j =0

18. MatMul on 1D Bus without Broadcast Naïve algorithm: C(myproc) = C(myproc) + A(myproc)*b(myproc,myproc) for i = 0 to p-1 for j = 0 to p-1 except i if (myproc == i) send A(i) to processor j // message passing if (myproc == j) receive A(i) from processor i C(myproc) = C(myproc) + A(i)*B(i,myproc) barrier Cost of inner loop: computation: 2*n*(n/p)2 = 2*n3/p2 communication: a*p2 + b*p*n2 // approximately

19. Naïve MatMul (continued) Cost of inner loop: computation: 2*n*(n/p)2 = 2*n3/p2 communication: a*p2 + b*p*n2 // approximately Only 1 pair of processors (i and j) are active on any iteration, an of those, only i is doing computation => the algorithm is almost entirely serial Running time: (p*(p-1) + 1)*computation + p*(p-1)*communication ~= 2*n3 + p2*a + p*n2*b this is worse than the serial time and grows with p

20. Better MatMul on a Bus Remove the barrier and send A(i) multiple times C(myproc) = C(myproc) + A(myproc)*b(myproc,myproc) for i = 0 to myproc-1 receive A(i) from processor i C(myproc) = C(myproc) + A(i)*B(i,myproc) for i = 0 to p-1 except myproc send A(myproc) to processor i for i = myproc+1 to p-1 receive A(i) from processor i C(myproc) = C(myproc) + A(i)*B(i,myproc) Program is indeterminate: sends/receives may happen in different orders

21. Performance of “Better” Algorithm • Intuitively, if a computation step is sufficient large compared to communication, efficiency will be good • communication: a + n*(n/p)*b • computation: 2*n3/p2 • Assume computation > communication i-j represents communication from i to j iC represents computation step by processor i

22. Performance of “Better” Algorithm • Timeline of execution • If computation<= (p-2)*communication • No bubbles in the pipeline • Time = p*(p-1)*communication + 2*computation • Time <= (p2 + p-4)*communication • If communication = computation/(p-2), close to lower bound of 2*n3/p • If communication is faster, only small impact on performance 0C 0-1 0-2 0-3 0- p-1 1C 1C 1- 0 2C 2C 2-0 p-1C

23. MatMul: 1D Layout and Broadcast • Modify the previous algorithm so that each of sends of A(i) is a broadcast (assumed 1 comm. step) • Time is now 2*n3/p + p*a +n2*b • Again, we require n>>p for good efficiency • p times less communication time • broadcast helps performance (as expected)

24. MatMul with 2D Layout • Consider processors in 2D grid (physical or logical) p(0,0) p(0,1) p(0,2) p(1,0) p(1,1) p(1,2) p(2,0) p(2,1) p(2,2)

25. Cannon’s Algorithm k<s k=0 • C(i,j) = C(i,j) + S A(i,k)*B(k,j) • Algorithm (s = sqrt(p)): • for i=0 to s-1 // skew A • left-circular-shift row i of A by i • so that A(i,j) overwritten by A(i,(j+i)mod s) • for i=0 to s-1 // skew B • up-circular shift B • so that B(i,j) overwritten by B((i+j)mod s), j) • for k=0 to s-1 • for i=0 to s-1 and j=0 to s-1 • C(i,j) = C(i,j) + A(i,j)*B(i,j) • left-circular-shift each row of A by 1 • up-circular-shift each row of B by 1

26. Communication in Cannon

27. BLAS Review

28. Fast linear algebra kernels: BLAS • Simple linear algebra kernels such as matrix-matrix multiply • More complicated algorithms can be built from these basic kernels. • The interfaces of these kernels have been standardized as the Basic Linear Algebra Subroutines (BLAS). • Early agreement on standard interface (~1980) • Led to portable libraries for vector and shared memory parallel machines. • On distributed memory, there is a less-standard interface called the PBLAS

29. Level 1 BLAS • Operate on vectors or pairs of vectors • perform O(n) operations; • return either a vector or a scalar. • saxpy • y(i) = a * x(i) + y(i), for i=1 to n. • s stands for single precision, daxpy is for double precision, caxpy for complex, and zaxpy for double complex, • sscal y = a * x, for scalar a and vectors x,y • sdot computes s = Sni=1 x(i)*y(i)

30. Level 2 BLAS • Operate on a matrix and a vector; • return a matrix or a vector; • O(n^2) operations • sgemv: matrix-vector multiply • y = y + A*x • where A is m-by-n, x is n-by-1 and y is m-by-1. • sger: rank-one update • A = A + y*x', i.e., A(i,j) = A(i,j)+y(i)*x(j) • where A is m-by-n, y is m-by-1, x is n-by-1, x' is x transpose • strsv: triangular solve • solves y=T*x for x, where T is triangular

31. Level 3 BLAS • Operate on pairs or triples of matrices • returning a matrix; • complexity is O(n**3). • sgemm: Matrix-matrix multiplication • C = C +A*B, • where C is m-by-n, A is m-by-k, and B is k-by-n • sgtrsm: multiple triangular solve • solves Y = T*X for X, • where T is a triangular matrix, and X is a rectangular matrix.

32. Performance of BLAS Level 3 Level 2 Level 1

33. Performance of BLAS • BLAS are specially optimized by the vendor • Sun BLAS uses features in the Ultrasparc • Big payoff for algorithms that can be expressed in terms of the BLAS3 instead of BLAS2 or BLAS1. • The top speed of the BLAS3 • Algorithms like Gaussian elimination organized so that they use BLAS3