Cs 267 applications of parallel processors lecture 13 parallel matrix multiply
This presentation is the property of its rightful owner.
Sponsored Links
1 / 33

CS 267 Applications of Parallel Processors Lecture 13: Parallel Matrix Multiply PowerPoint PPT Presentation


  • 97 Views
  • Uploaded on
  • Presentation posted in: General

CS 267 Applications of Parallel Processors Lecture 13: Parallel Matrix Multiply. Kathy Yelick http://www.cs.berkeley.edu/~dmartin/cs267. Outline. - Recap - Sources of large dense linear systems - BLAS for linear algebra - Parallel Matrix multiply. Model overview. Work-depth PRAM

Download Presentation

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

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

CS 267 Applications of Parallel ProcessorsLecture 13: Parallel Matrix Multiply

Kathy Yelick

http://www.cs.berkeley.edu/~dmartin/cs267


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

Outline

- Recap

- Sources of large dense linear systems

- BLAS for linear algebra

- Parallel Matrix multiply


Model overview

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


Dense linear algebra in electromagentics

Dense Linear Algebra in Electromagentics


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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)


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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/


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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)


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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.


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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.)


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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.


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


Parallel matrix multiply

Parallel Matrix Multiply


Parallel matrix multiply1

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


1d layout

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


Matrix multiply 1d layout on bus

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


Matmul on 1d bus without broadcast

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


Na ve matmul continued

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


Better matmul on a bus

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


Performance of better algorithm

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


Performance of better algorithm1

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


Matmul 1d layout and broadcast

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)


Matmul with 2d layout

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)


Cannon s algorithm

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


Communication in cannon

Communication in Cannon


Blas review

BLAS Review


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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)


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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.


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

Performance of BLAS

Level 3

Level 2

Level 1


Cs 267 applications of parallel processors lecture 13 parallel matrix multiply

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


  • Login