# Presentation: 153127

## CS 267: Applications of Parallel ComputersLecture 15:Parallel Matrix Multiplication

- Horst D. Simon
**http://www.cs.berkeley.edu/~strive/cs267**

- CS267, Yelick

## Parallel Numerical Algorithms	

- Lecture schedule:
- 10/15: Dense Matrix Products
- BLAS 1: Vector operations
- BLAS 2: Matrix-Vector operations
- BLAS 3: Matrix-Matrix operations
- Use of Performance models in algorithm design
- 10/17: Dense Matrix Solvers (Jim Demmel)
- 10/22: Sparse Matrix Products
- 10/24: Sparse Direct Solvers (Esmond Ng)
- 10/29: Iterative Solvers and Preconditioners (Esmond Ng)

- CS267, Yelick

**Notes:**
Not all problems solved on parallel machines involve simulation, but many do

## Parallel Vector Operations

- Some common vector operations for vectors x,y,z:
- Vector add: z = x + y
- Trivial to parallelize if vectors are aligned
- AXPY: z = a*x+y (where a is scalar)
- Broadcast a, followed by independent * and +
- Dot product:  s = j x[j] * y[j]
- Independent * followed by + reduction

- CS267, Yelick

## Broadcast and reduction

- Broadcast of 1 value to p processors in log p time
- Reduction of p values to 1 in log p time
- Takes advantage of associativity in +,*, min, max, etc.

**8**

**  ****1   3  1  0  4 -6 3   2**

- Add-reduction

- Broadcast

- CS267, Yelick

## Broadcast Algorithms

- Sequential or “centralized” algorithm
- P0 sends value to P-1 other processors in sequence
- O(P) algorithm
- Note: variations in UPC/Titanium model based on whether P0 writes to all others, or others read from P0
- Tree-based algorithm
- May vary branching factor
- O(log P) algorithm
- If broadcasting large data blocks, may break into pieces and pipeline

- Broadcast

- P0

- P4

- P6

- P0  P1 P2 P3 P4 P5 P6 P7

- CS267, Yelick

## Lower Bound on Parallelism

- Given only binary operations on our machine.
- In 1 time step, only 1 operation can be performed
- output depends on at most 2 inputs
- In 2 time steps, output depends on at most 4 inputs
- Adding  a time step increases possible inputs by at most 2x
- In k=log n time steps, output depends on at most n=2k inputs
- A function of n inputs requires at least log n parallel steps.

- CS267, Yelick

## Scan (or Parallel prefix), A Digression

- Definition: the Scan operation take a binary associative operator    , and an array of n elements
-                [a0, a1, a2, ... an-1]
-     and produces the array
-               [a0, (a0    a1), ... (a0     a1    ...     an-1)]
- Example: add scan of
-         [1, 2, 0, 4, 2, 1, 1, 3]  is  [1, 3, 3, 7, 9, 10, 11, 14]
- Can be implemented in O(n) time by a serial algorithm
- Generalized to “segmented scans” given a mask of positions at which total is reset.

- CS267, Yelick

## Applications of scans

- There are several applications of scans, some more obvious than others
- lexically compare string of characters
- add multi-precision numbers (represented as array of numbers)
- evaluate polynomials
- implement bucket sort and radix sort
- solve tridiagonal systems
- to dynamically allocate processors
- to search for regular expression (e.g., grep)

- CS267, Yelick

## E.g., Using Scans for Array Compression

- Given an array of n elements
-                [a0, a1, a2, ... an-1]
-     and an array of flags
-               [1,0,1,1,0,0,1,...]
-     compress the flagged elements
-               [a0, a2, a3, a6, ...]
- Compute a “prescan” i.e., a scan that doesn’t include the element at position i in the sum
- 		[0,1,1,2,3,3,4,...]
- Gives the index of the ith element in the compressed array

- CS267, Yelick

## Implementing Scans

- Naïve approach: use multiple trees to calculate partial sums
- works in O(log n) time
- but uses roughly n2/2 processors
- Tree summation 2 phases
- up sweep
- get values L and R from left and right child
- save L in local variable Mine
- compute Tmp = L + R and pass to parent
- down sweep
- get value Tmp from parent
- send Tmp to left child
- send Tmp+Mine to right child

- CS267, Yelick

## Implementing Scans: Example

- Tree summation is used to efficiently implement scans

**6**

**5**

**4**

**3**

**2**

**4**

**1**

**6**

**5**

**4**

**3**

**2**

**4**

**1**

**Up sweep:**
**  ****mine = left**
**  ****tmp = left + right**

**Down sweep:**
**  ****tmp = parent (root is 0)**
**  ****right = tmp + mine**

**4**

**6**

**9**

**5**

**4**

**0**

**6**

**0**

**0**

**3**

**3         4   6        6      10      11  12      15**

**3         1   2        0      4         1    1        3**

**+X = 3         1   2        0      4         1    1        3**

**4**

**4**

**6**

**6**

**10**

**11**

**6**

**11**

**12**

- CS267, Yelick

## Parallel Matrix-Vector Product

- y = y + A*x, where A is a dense  matrix
- Layout:
- 1D by rows
- Algorithm:
- Foreach processor i
-    Broadcast X(i)
-    Compute A(i)*x
- A(i) refers to the n by n/p block row that processor i owns
- Algorithm uses the formula
- Y(i) = Y(i) + A(i)*X = Y(i) + j A(i)*X(j)

- x

- y

- Po
- P1
- P2
- P3

- Po   P1    P2     P3

- CS267, Yelick

## Matrix-Vector Product

- A column layout of the matrix eliminates the broadcast
- But adds a reduction to update the destination
- A blocked layout uses a broadcast and reduction, both on a subset of processors

- P0     P1    P2    P3

- P0      P1     P2    P3

- P4      P5     P6    P7

- P8     P9     P10   P11

- P12    P13   P14   P15

- CS267, Yelick

## 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
- Time (message) = latency + #words * time-per-word
-                             =  + n*

- CS267, Yelick

## Latency Bandwidth Model

- Network of fixed number P of processors
- fully connected
- each with local memory
- Latency ()
- accounts for varying performance with number of messages
- Inverse bandwidth ()
- accounts for performance varying with volume of data
- Efficiency (in any model):
- serial time / (p *  parallel time)
- perfect (linear) speedup  efficiency = 1

- CS267, Yelick

## Matrix Multiply with 1D Column Layout

- Assume matrices are n x n 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
- Algorithm uses the formula
- C(i) = C(i) + A*B(i) = C(i) + j A(j)*B(j,i)

- CS267, Yelick

## Matrix Multiply: 1D Layout on Bus or Ring

- Algorithm uses the formula
- C(i) = C(i) + A*B(i) = C(i) + j A(j)*B(j,i)
- First consider a bus-connected machine without broadcast:  only one pair of processors can communicate at a time (ethernet)
- Second consider a machine with processors on a ring: all processors may communicate with nearest neighbors simultaneously

- CS267, Yelick

## MatMul: 1D layout on 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
-             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:  + *n2  /p

- CS267, Yelick

## Naïve MatMul (continued)

**Cost of inner loop**:
-        computation: 2*n*(n/p)2 = 2*n3/p2
-        communication:  + *n2 /p        ... approximately
- Only 1 pair of processors (i and j) are active on any iteration,
-   and 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* + p*n2*
-   this is worse than the serial time and grows with p

- CS267, Yelick

## Matmul for 1D layout on a Processor Ring

- Pairs of processors can communicate simultaneously

**Copy A(myproc) into Tmp**
**C(myproc) = C(myproc) + Tmp*B(myproc , myproc)**
**for j = 1 to p-1**
**      ****Send Tmp to processor myproc+1 mod p**
**      ****Receive Tmp from processor myproc-1 mod p**
**      ****C(myproc) = C(myproc) + Tmp*B( myproc-j mod p , myproc)**

-  Same idea as for gravity in simple sharks and fish algorithm
- May want double buffering in practice for overlap
- Ignoring deadlock details in code
-  Time  of inner loop = 2*( + *n2/p) + 2*n*(n/p)2

- CS267, Yelick

## Matmul for 1D layout on a Processor Ring

- Time  of inner loop = 2*(a + b*n2/p) + 2*n*(n/p)2
- Total Time  = 2*n* (n/p)2  +  (p-1) * Time of inner loop
-                      ~ 2*n3/p  + 2*p*a + 2*b*n2
- Optimal for 1D layout on Ring or Bus, even with with Broadcast:
-  Perfect speedup for arithmetic
-  A(myproc) must move to each other processor, costs at least
-                (p-1)*cost of sending n*(n/p) words
- Parallel Efficiency = 2*n3 / (p * Total Time)
-                                 = 1/(1 +  * p2/(2*n3) +  * p/(2*n) )
-                                 = 1/ (1 + O(p/n))
-  Grows to 1 as n/p increases (or  and  shrink)

- CS267, Yelick

## MatMul with 2D Layout

- Consider processors in 2D grid (physical or logical)
- Processors can communicate with 4 nearest neighbors
- Broadcast along rows and columns
- Assume p is square s x s grid

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

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

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

- =

- *

- CS267, Yelick

## Cannon’s Algorithm

- ... C(i,j) = C(i,j) +   A(i,k)*B(k,j)
- ...  assume s = sqrt(p) is an integer
-    forall  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)
-    forall  i=0 to s-1               ...  “skew” B
-          up-circular-shift B column i of B by i
-           ... so that B(i,j) overwritten by B((i+j)mod s), j)
-    for k=0 to s-1        ... sequential
-           forall i=0 to s-1 and j=0 to s-1    ... all processors in parallel
-                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

**k**

- CS267, Yelick

**C(1,2) = ****A(1,0)**** * ****B(0,2)**** + ****A(1,1)**** * ****B(1,2)**** + ****A(1,2)**** * ****B(2,2)**

- CS267, Yelick

## Initial Step to Skew Matrices in Cannon

- Initial blocked input
- After skewing before initial block multiplies

- A(0,1)

- A(0,2)

- A(1,0)

- A(2,0)

- A(1,1)

- A(1,2)

- A(2,1)

- A(2,2)

- A(0,0)

- B(0,1)

- B(0,2)

- B(1,0)

- B(2,0)

- B(1,1)

- B(1,2)

- B(2,1)

- B(2,2)

- B(0,0)

- A(0,1)

- A(0,2)

- A(1,0)

- A(2,0)

- A(1,1)

- A(1,2)

- A(2,1)

- A(2,2)

- A(0,0)

- B(0,1)

- B(0,2)

- B(1,0)

- B(2,0)

- B(1,1)

- B(1,2)

- B(2,1)

- B(2,2)

- B(0,0)

- CS267, Yelick

## Skewing Steps in Cannon

- First step
- Second
- Third

- A(0,1)

- A(0,2)

- A(1,0)

- A(2,0)

- A(1,2)

- A(2,1)

- B(0,1)

- B(0,2)

- B(1,0)

- B(2,0)

- B(1,1)

- B(1,2)

- B(2,1)

- B(2,2)

- B(0,0)

- A(0,1)

- A(0,2)

- A(1,0)

- A(2,0)

- A(1,1)

- A(1,2)

- A(2,1)

- A(2,2)

- A(0,0)

- B(0,1)

- B(0,2)

- B(1,0)

- B(2,0)

- B(1,1)

- B(1,2)

- B(2,1)

- B(2,2)

- B(0,0)

- A(1,1)

- A(2,2)

- A(0,0)

- CS267, Yelick

## Cost of Cannon’s Algorithm

-   forall  i=0 to s-1              ...  recall s = sqrt(p)
-          left-circular-shift row i of A by i    ... cost = s*( + *n2/p)
-    forall  i=0 to s-1
-          up-circular-shift B column i of B by i ... cost = s*( + *n2/p)
-    for k=0 to s-1
-           forall  i=0 to s-1 and j=0 to s-1
-                C(i,j) = C(i,j) + A(i,j)*B(i,j)   ... cost = 2*(n/s)3 = 2*n3/p3/2
-                left-circular-shift each row of A by 1   ... cost = + *n2/p
-                up-circular-shift each row of B by 1     ... cost =  + *n2/p

** ****Total Time = 2*n****3****/p +  4********* ****s************* + 4**********n****2****/s  **
** ****Parallel Efficiency = 2*n****3**** / (p * Total Time)**
**                                  ****= 1/( 1 + ******** * 2*(s/n)****3**** + ******** * 2*(s/n) )**
**                                  ****= 1/(1 + O(sqrt(p)/n)) **
** ****Grows to 1 as n/s = n/sqrt(p) = sqrt(data per processor) grows**
** ****Better than 1D layout, which had Efficiency = 1/(1 + O(p/n))**

- CS267, Yelick

## Drawbacks to Cannon

- Hard to generalize for
- p not a perfect square
- A and B not square
- Dimensions of A, B not perfectly divisible by s=sqrt(p)
- A and B not “aligned” in the way they are stored on processors
- block-cyclic layouts
- Memory hog (extra copies of local matrices)

- CS267, Yelick

## SUMMA Algorithm

- SUMMA = Scalable Universal Matrix Multiply
- Slightly less efficient, but simpler and easier to generalize
- Presentation from van de Geijn and Watts
- www.netlib.org/lapack/lawns/lawn96.ps
- Similar ideas appeared many times
- Used in practice in PBLAS = Parallel BLAS
- www.netlib.org/lapack/lawns/lawn100.ps

- CS267, Yelick

## SUMMA 

*****

** ****=**

**I**

**J**

**A(I,k)**

**k**

**k**

**B(k,J)**

-   I, J represent all rows, columns owned by a processor
-  k is a single row or column
-  or a block of b rows or columns
- C(I,J) = C(I,J) + k A(I,k)*B(k,J)
-  Assume a pr by pc processor grid (pr = pc = 4 above)
- Need not be square

- C(I,J)

- CS267, Yelick

## SUMMA 

- For k=0 to n-1     ... or n/b-1 where b is the block size
-                             ...  = # cols in A(I,k) and # rows in B(k,J)
-      for all I = 1 to pr   ... in parallel
-            owner of A(I,k) broadcasts it to whole processor row
-      for all J = 1 to pc  ... in parallel
-             owner of B(k,J) broadcasts it to whole processor column
-      Receive A(I,k) into Acol
-      Receive B(k,J) into Brow
-      C( myproc , myproc ) = C( myproc , myproc) + Acol * Brow

*****

** ****=**

**I**

**J**

**A(I,k)**

**k**

**k**

**B(k,J)**

- C(I,J)

- CS267, Yelick

## SUMMA performance 

**For k=0 to n/b-1**
**     ****for all I = 1 to s   ****...  s = sqrt(p)**
**           ****owner of A(I,k) broadcasts it to whole processor row**
**               ****... ****time = log s *( ******** + ******** * b*n/s), using a tree**
**     ****for all J = 1 to  s**
**            ****owner of B(k,J) broadcasts it to whole processor column**
**              ****... ****time = log s *( ******** + ******** * b*n/s), using a tree**
**     ****Receive A(I,k) into Acol**
**     ****Receive B(k,J) into Brow**
**     ****C( myproc , myproc ) = C( myproc , myproc) + Acol * Brow**
**              ****... ****time = 2*(n/s)****2*****b**

- Total time = 2*n3/p   +    * log p * n/b   +   * log p * n2 /s

- To simplify analysis only, assume s = sqrt(p)

- CS267, Yelick

## SUMMA performance 

-  Total time = 2*n3/p   +   a * log p * n/b   +   b * log p * n2 /s
-  Parallel Efficiency =
-        1/(1 +  * log p * p / (2**n2)  +  * log p * s/(2*n) )
-  ~Same  term as Cannon, except for log p factor
-           log p grows slowly so  this is ok
-  Latency (a) term can be larger, depending on b
-           When b=1, get   * log p * n
-           As b grows to n/s, term shrinks to
-                      * log p * s  (log p times Cannon)
-  Temporary storage grows like 2*b*n/s
-  Can change b to tradeoff latency cost with memory

- CS267, Yelick

## ScaLAPACK Parallel Library

**ScaLAPACK Parallel Library**

- CS267, Yelick

## PDGEMM =  PBLAS routine

**PDGEMM =  PBLAS routine**
**      ****for matrix multiply**
**Observations:**
**    ****For fixed N, as P increases**
**       ****Mflops increases, but**
**       ****less than 100% efficiency**
**    ****For fixed P, as N increases,**
**       ****Mflops (efficiency) rises**

**DGEMM = BLAS routine**
**      ****for matrix multiply**
**Maximum speed for PDGEMM**
**     ****= # Procs * speed of DGEMM**
**Observations (same as above):**
**     ****Efficiency always at least 48%**
**     ****For fixed N, as P increases,**
**         ****efficiency drops **
**     ****For fixed P, as N increases,**
**         ****efficiency increases**

- CS267, Yelick

## Recursive Layouts

- For both cache hierarchies and parallelism, recursive layouts may be useful
- Z-Morton, U-Morton, and X-Morton Layout
- Also Hilbert layout and others
- What about the user’s view?
- Fortunately, many problems can be solved on a permutation
- Never need to actually change the user’s layout

- CS267, Yelick

## Summary of Parallel Matrix Multiplication

- 1D Layout
- Bus without broadcast - slower than serial
- Nearest neighbor communication on a ring (or bus with broadcast): Efficiency = 1/(1 + O(p/n))
- 2D Layout
- Cannon
- Efficiency = 1/(1+O(p1/2 /n))
- Hard to generalize for general p, n, block cyclic, alignment
- SUMMA
- Efficiency = 1/(1 + O(log p * p / (b*n2) + log p * p1/2 /n))
- Very General
- b small => less memory, lower efficiency
- b large => more memory, high efficiency
- Recursive layouts
- Efficiency = 1/(1 + O(p1/3 /n) )  ??

- CS267, Yelick

## Extra Slides

- CS267, Yelick

## Review: BLAS 3 (Blocked) GEPP

**for   ib = 1 to n-1 step b     ****... Process matrix b columns at a time**
**     ****end = ib + b-1                ... Point to end of block of b columns **
**     ****apply BLAS2 version of GEPP to  get A(ib:n , ib:end) = P’ * L’ * U’**
**     ****... ****let LL denote the strict lower triangular part of A(ib:end , ib:end) + I**
**     ****A(ib:end , end+1:n)**** = LL****-1**** * ****A(ib:end , end+1:n)****         ****... update next b rows of U**
**     ****A(end+1:n , end+1:n )**** = ****A(end+1:n , end+1:n )**
**                  ****- ****A(end+1:n , ib:end)**** * ****A(ib:end , end+1:n)****    **
**                                       ****... ****apply delayed updates with single matrix-multiply**
**                                       ... ****with inner dimension b**

**BLAS 3**

- CS267, Yelick

## Review: Row and Column Block Cyclic Layout

**Review: Row and Column Block Cyclic Layout**

**processors and matrix blocks**
**are distributed in a 2d array**
**pcol-fold parallelism**
**in any column, and calls to the **
**BLAS2 and BLAS3 on matrices of **
**size brow-by-bcol**
**serial bottleneck is eased**
**need not be symmetric in rows and**
**columns**

- CS267, Yelick

## Distributed GE with a 2D Block Cyclic Layout

**Distributed GE with a 2D Block Cyclic Layout**

**block size b in the algorithm and the block sizes brow **
**and bcol in the layout satisfy b=brow=bcol. **
**shaded regions indicate busy processors or **
**communication performed.**
**unnecessary to have a barrier between each **
**step of the algorithm, e.g.. step 9, 10, and 11 can be **
**pipelined**

- CS267, Yelick

## Distributed GE with a 2D Block Cyclic Layout

**Distributed GE with a 2D Block Cyclic Layout**

- CS267, Yelick

## Matrix multiply of

**Matrix multiply of**
** ****green = green - blue * pink**

- CS267, Yelick

##           expect it to be (2*n3)/(2/3*n3) = 3

**PDGESV = ScaLAPACK **
**       ****parallel LU routine**
**Since it can run no faster than its**
**    ****inner loop (PDGEMM), we measure:**
**Efficiency = **
**     ****Speed(PDGESV)/Speed(PDGEMM)**
**Observations:**
**     ****Efficiency well above 50% for large**
**          ****enough problems**
**     ****For fixed N, as P increases,**
**         ****efficiency decreases **
**          ****(just as for PDGEMM)**
**     ****For fixed P, as N increases**
**         ****efficiency increases**
**          ****(just as for PDGEMM)**
**     ****From bottom table, cost of solving**
**          ****Ax=b about half of matrix multiply**
**          ****for large enough matrices.**
**          ****From the flop counts we would**
**          ****expect it to be (2*n****3****)/(2/3*n****3****) = 3**
**          ****times faster, but communication**
**          ****makes it a little slower.**

- CS267, Yelick

## Slide 45

- CS267, Yelick

## Scales well, 

**Scales well, **
**   ****nearly full machine speed**

- CS267, Yelick

## Old Algorithm,

**Old version,**
**pre 1998 Gordon Bell Prize**
**Still have ideas to accelerate**
**Project Available!**

**Old Algorithm,**
**    ****plan to abandon**

- CS267, Yelick

## Have good ideas to speedup

**Have good ideas to speedup**
**Project available!**

**Hardest of all to parallelize**
**Have alternative,  and**
**      ****would like to compare**
**Project available!**

- CS267, Yelick

## Out-of-core means

**Out-of-core means**
**  ****matrix lives on disk;**
**  ****too big for main mem**
**Much harder to hide**
**  ****latency of disk**
**QR much easier than LU**
**  ****because no pivoting**
**  ****needed for QR**
**Moral: use QR to solve Ax=b**
**Projects available**
**  ****(perhaps very hard...)**

- CS267, Yelick

## A small software project ...

**A small software project ...**

- CS267, Yelick

## Work-Depth Model of Parallelism

- The work depth model:
- The simplest model is used
- For algorithm design, independent of a machine
- The work, W, is the total number of operations
- The depth, D, is the longest chain of dependencies
- The parallelism, P, is defined as W/D
- Specific examples include:
- circuit model, each input defines a graph with ops at nodes
- vector model, each step is an operation on a vector of elements
- language model, where set of operations defined by language

- CS267, Yelick