# Presentation: 153130

## CS 267: Applications of Parallel ComputersLecture 3:Matrix Multiplication(cont.)

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

- CS267, Yelick

## Recap of Last Lecture

- Scaled speedup: operate at the memory boundary
- Memory systems on modern processors are complicated.
- The performance of a simple program can depend on the details of the micro-architecture.

- CS267, Yelick

## Outline

- Matrix Multiply Cache Optimizations
- Bag of Tricks
- Parallel Architectures and Programming Models

- CS267, Yelick

## Note on Matrix Storage

- A matrix is a 2-D array of elements, but memory addresses are “1-D”
- Conventions for matrix layout
- by column, or “column major” (Fortran default)
- by row, or “row major” (C default)

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 9

- 10

- 11

- 12

- 13

- 14

- 15

- 16

- 17

- 18

- 19

- 0

- 4

- 8

- 12

- 16

- 1

- 5

- 9

- 13

- 17

- 2

- 6

- 10

- 14

- 18

- 3

- 7

- 11

- 15

- 19

- Column major

- Row major

- CS267, Yelick

## Assume just 2 levels in the hierarchy, fast and slow
All data initially in slow memory
m = number of memory elements (words) moved between fast and slow memory 
tm = time per slow memory operation
f = number of arithmetic operations
tf = time per arithmetic operation << tm
q = f / m  average number of flops per slow element access
Minimum possible time = f* tf when all data in fast memory
Actual time 

Larger q means Time closer to minimum f * tf 

- All data initially in slow memory
- m = number of memory elements (words) moved between fast and slow memory
- tm = time per slow memory operation
- f = number of arithmetic operations
- tf = time per arithmetic operation << tm
- q = f / m  average number of flops per slow element access
- Minimum possible time = f* tf when all data in fast memory
- Actual time
- Larger q means Time closer to minimum f * tf

** ****Using a Simple Model of Memory to Optimize**

- f * tf + m * tm = f * tf * (1 + tm/tf  * 1/q)

- CS267, Yelick

## Warm up: Matrix-vector multiplication

- {implements y = y + A*x}
- for i = 1:n
- 		for j = 1:n
- 			y(i) = y(i) + A(i,j)*x(j)

- =

- +

- *

- y(i)

- y(i)

- A(i,:)

- x(:)

- CS267, Yelick

## Warm up: Matrix-vector multiplication

- {read x(1:n) into fast memory}
- {read y(1:n) into fast memory}
- for i = 1:n
- 	    {read row i of A into fast memory}
- 	    for j = 1:n
- 		   y(i) = y(i) + A(i,j)*x(j)
- {write y(1:n) back to slow memory}

-  m = number of slow memory refs = 3n + n2
-  f   = number of arithmetic operations = 2n2
-  q  = f / m ~= 2
-  Matrix-vector multiplication limited by slow memory speed

- CS267, Yelick

## “Naïve” Matrix Multiply

- {implements C = C + A*B}
- for i = 1 to n
-        for j = 1 to n
- 		for k = 1 to n
-            		C(i,j) = C(i,j) + A(i,k) * B(k,j)

**=**

**+**

*****

- C(i,j)

- C(i,j)

- A(i,:)

- B(:,j)

- CS267, Yelick

## “Naïve” Matrix Multiply

- {implements C = C + A*B}
- for i = 1 to n
-   {read row i of A into fast memory}
-    for j = 1 to n
-        {read C(i,j) into fast memory}
-        {read column j of B into fast memory}
-        for k = 1 to n
-            C(i,j) = C(i,j) + A(i,k) * B(k,j)
-        {write C(i,j) back to slow memory}

**=**

**+**

*****

- C(i,j)

- A(i,:)

- B(:,j)

- C(i,j)

- CS267, Yelick

## “Naïve” Matrix Multiply

- Number of slow memory references on unblocked matrix multiply
- 	m = n3 read each column of B  n  times
-          + n2 read each row of A once
-          + 2n2 read and write each element of C once
-         = n3 + 3n2
- So q = f / m = 2n3 / (n3 + 3n2)
-         ~= 2 for large n, no improvement over matrix-vector multiply

**=**

**+**

*****

- C(i,j)

- C(i,j)

- A(i,:)

- B(:,j)

- CS267, Yelick

## Blocked (Tiled) Matrix Multiply

- Consider A,B,C to be N by N matrices of b by b subblocks where b=n / N is called the block size
- 	   for i = 1 to N
-  	      for j = 1 to N
-        	{read block C(i,j) into fast memory}
-        	for k = 1 to N
-            	       {read block A(i,k) into fast memory}
-            	       {read block B(k,j) into fast memory}
-           	        C(i,j) = C(i,j) + A(i,k) * B(k,j) {do a matrix multiply on blocks}
-       	 {write block C(i,j) back to slow memory}

**=**

**+**

*****

- C(i,j)

- C(i,j)

- A(i,k)

- B(k,j)

- CS267, Yelick

## Blocked (Tiled) Matrix Multiply

- Recall:
-    m is amount memory traffic between slow and fast memory
-    matrix has nxn elements, and NxN blocks each of size bxb
-    f is number of floating point operations, 2n3 for this problem
-    q = f / m is our measure of algorithm efficiency in the memory system
- So:
- 	m =  N*n2    read each block of B  N3 times (N3 * n/N * n/N)
-          + N*n2    read each block of A  N3 times
-          + 2n2       read and write each block of C once
-         =  (2N + 2) * n2
- So q = f / m = 2n3 / ((2N + 2) * n2)
-         ~= n / N = b  for large n
- So we can improve performance by increasing the blocksize b
- Can be much faster than matrix-vector multiply (q=2)

- CS267, Yelick

## Limits to Optimizing Matrix Multiply

- The blocked algorithm has ratio q ~= b
-    The large the block size, the more efficient our algorithm will be
-    Limit:   All three blocks from A,B,C must fit in fast memory (cache), so we cannot make these blocks arbitrarily large:
-                   3b2 <= M, so q ~= b <= sqrt(M/3)
- There is a lower bound result that says we cannot do any better than this (using only algebraic associativity)
- Theorem (Hong & Kung, 1981): Any reorganization of this algorithm
- (that uses only algebraic associativity) is limited to q = O(sqrt(M))

- CS267, Yelick

## Basic Linear Algebra Subroutines

- Industry standard interface (evolving)
- Vendors, others supply optimized implementations
- History
- BLAS1 (1970s):
- vector operations: dot product, saxpy (y=*x+y), etc
- m=2*n, f=2*n, q ~1 or less
- BLAS2 (mid 1980s)
- matrix-vector operations: matrix vector multiply, etc
- m=n^2, f=2*n^2, q~2, less overhead
- somewhat faster than BLAS1
- BLAS3 (late 1980s)
- matrix-matrix operations: matrix matrix multiply, etc
- m >= 4n^2, f=O(n^3), so q can possibly be as large as n, so BLAS3 is potentially much faster than BLAS2
- Good algorithms used BLAS3 when possible (LAPACK)
- See www.netlib.org/blas, www.netlib.org/lapack

- CS267, Yelick

## BLAS speeds on an IBM RS6000/590

- BLAS 3

- BLAS 2

- BLAS 1

- BLAS 3 (n-by-n matrix matrix multiply) vs BLAS 2 (n-by-n matrix vector multiply) vs BLAS 1 (saxpy of  n vectors)

- Peak speed = 266 Mflops

- Peak

- CS267, Yelick

## Locality in Other Algorithms

- The performance of any algorithm is limited by q
- In matrix multiply, we increase q by changing computation order
- increased temporal locality
- For other algorithms and data structures, even hand-transformations are still an open problem
- sparse matrices (reordering, blocking)
- trees (B-Trees are for the disk level of the hierarchy)
- linked lists (some work done here)

- CS267, Yelick

## Tiling Alone Might Not Be Enough

- Naïve and a “naïvely tiled” code

- CS267, Yelick

## Optimizing in Practice

- Tiling for registers
- loop unrolling, use of named “register” variables
- Tiling for multiple levels of cache
- Exploiting fine-grained parallelism in processor
- superscalar; pipelining
- Complicated compiler interactions
- Hard to do by hand (but you’ll try)
- Automatic optimization an active research area
- BeBOP: www.cs.berkeley.edu/~richie/bebop
- PHiPAC: www.icsi.berkeley.edu/~bilmes/phipac
-           in particular tr-98-035.ps.gz
- ATLAS: www.netlib.org/atlas

- CS267, Yelick

## PHiPAC: Portable High Performance ANSI C

- Speed of n-by-n matrix multiply on Sun Ultra-1/170, peak = 330 MFlops

- CS267, Yelick

## ATLAS (DGEMM n = 500)

- ATLAS is faster than all other portable BLAS implementations and it is comparable with machine-specific libraries provided by the vendor.

- Source: Jack Dongarra

- CS267, Yelick

## Recap

- Simple performance models can aid in understanding
- Two ratios are key to efficiency
- algorithmic: q = f/m = # floating point opns / # slow memory opns
- tm/tf = time for slow memory operation / time for floating point operation
- A common technique for improving cache performance (lowering q) is called blocking
- Applied to matrix multiplication.

- CS267, Yelick

## Search Over Block Sizes

- Performance models are useful for high level algorithms
- Helps in developing a blocked algorithm
- Models have not proven very useful for block size selection
- too complicated to be useful
- See work by Sid Chatterjee for detailed model
- too simple to be accurate
- Multiple multidimensional arrays, virtual memory, etc.
- Some systems use search
- Atlas – being incorporated into Matlab
- BeBOP – http://www.cs.berkeley.edu/~richie/bebop

- CS267, Yelick

## What the Search Space Looks Like

- A 2-D slice of a 3-D register-tile search space. The dark blue region was pruned.
- (Platform: Sun Ultra-IIi, 333 MHz, 667 Mflop/s peak, Sun cc v5.0 compiler)

- Number of columns in register block

- Number of rows in register block

- CS267, Yelick

## Strassen’s Matrix Multiply

- The traditional algorithm (with or without tiling) has O(n^3) flops
- Strassen discovered an algorithm with asymptotically lower flops
- O(n^2.81)
- Consider a 2x2 matrix multiply
- normally 8 multiplies, Strassen does it with 7 multiplies (but many more adds)

- Let M = [m11 m12] = [a11 a12] * [b11 b12]
-              [m21 m22] = [a21 a22]   [b21 b22]
- Let p1 = (a12 - a22) * (b21 + b22)                               p5 = a11 * (b12 - b22)
-       p2 = (a11 + a22) * (b11 + b22)                              p6 = a22 * (b21 - b11)
-       p3 = (a11 - a21) * (b11 + b12)                               p7 = (a21 + a22) * b11
-       p4 = (a11 + a12) * b22
- Then  m11 = p1 + p2 - p4 + p6
-           m12 = p4 + p5
-           m21 = p6 + p7
-           m22 = p2 - p3 + p5 - p7

- Extends to nxn by divide&conquer

- CS267, Yelick

## Strassen (continued)

- Asymptotically faster
- Several times faster for large n in practice
- Cross-over depends on machine
- Available in several libraries
- Caveats
- Needs more memory than standard algorithm
- Can be less accurate because of roundoff error
- Current world’s record is O(n 2.376.. )
- Why does Hong/Kung theorem not apply?

- CS267, Yelick

**Notes:**
World record due to Winograd and Copperfield
  Very complicated algorithm
  Needs unrealistically enormous matrices

## Outline

- Parallelism in Modern Processors
- Memory Hierarchies
- Matrix Multiply Cache Optimizations
- Bag of Tricks

- CS267, Yelick

## Removing False Dependencies

- Using local variables, reorder operations to remove false dependencies

**a[i] = b[i] + c;**
**a[i+1] = b[i+1] * d;**

**float f1 = b[i];**
**float f2 = b[i+1];**
**a[i] = f1 + c;**
**a[i+1] = f2 * d;**

- false read-after-write hazard
- between a[i] and b[i+1]

- With some compilers, you can declare a and b unaliased.
- Done via “restrict pointers,” compiler flag, or pragma)

- CS267, Yelick

## Exploit Multiple Registers

- Reduce demands on memory bandwidth by pre-loading into local variables

**while( ... ) {**
**   *****res++ = filter[0]*signal[0]**
**            ****+ filter[1]*signal[1]**
**            ****+ filter[2]*signal[2];**
**   ****signal++;**
**}**

**float f0 = filter[0];**
**float f1 = filter[1];**
**float f2 = filter[2];**
**while( ... ) {**
**   *****res++ = f0*signal[0]**
**            ****+ f1*signal[1]**
**            ****+ f2*signal[2];**
**   ****signal++;**
**}**

- also: register float f0 = ...;

- Example is a convolution

- CS267, Yelick

## Minimize Pointer Updates

- Replace pointer updates for strided memory addressing with constant array offsets

**f0 = *r8; r8 += 4;**
**f1 = *r8; r8 += 4;**
**f2 = *r8; r8 += 4;**

**f0 = r8[0];**
**f1 = r8[4];**
**f2 = r8[8];**
**r8 += 12;**

- Pointer vs. array expression costs may differ.
- Some compilers do a better job at analyzing one than the other

- CS267, Yelick

## Loop Unrolling

- Expose instruction-level parallelism

**float f0 = filter[0], f1 = filter[1], f2 = filter[2];**
**float s0 = signal[0], s1 = signal[1], s2 = signal[2];**
***res++ = f0*s0 + f1*s1 + f2*s2;**
**do {**
**   ****signal += 3;**
**   ****s0 = signal[0];**
**   ****res[0] = f0*s1 + f1*s2 + f2*s0;**
**   ****s1 = signal[1];**
**   ****res[1] = f0*s2 + f1*s0 + f2*s1;**
**   ****s2 = signal[2];**
**   ****res[2] = f0*s0 + f1*s1 + f2*s2;**
**   ****res += 3;**
**} while( ... );**

- CS267, Yelick

## Expose Independent Operations

- Hide instruction latency
- Use local variables to expose independent operations that can execute in parallel or in a pipelined fashion
- Balance the instruction mix (what functional units are available?)

**f1 = f5 * f9;**
**f2 = f6 + f10;**
**f3 = f7 * f11;**
**f4 = f8 + f12;**

- CS267, Yelick

## Copy optimization

- Copy input operands or blocks
- Reduce cache conflicts
- Constant array offsets for fixed size blocks
- Expose page-level locality

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 9

- 10

- 11

- 12

- 13

- 14

- 15

- Original matrix
- (numbers are addresses)

- 0

- 1

- 4

- 5

- 2

- 3

- 6

- 7

- 8

- 9

- 12

- 14

- 10

- 11

- 13

- 15

- Reorganized into
- 2x2 blocks

- CS267, Yelick

**Notes:**
This is one of the optimizations that Atlas performance that PhiPAC does not.

## Summary

- Performance programming on uniprocessors requires
- understanding of fine-grained parallelism in processor
- produce good instruction mix
- understanding of memory system
- levels, costs, sizes
- improve locality
- Blocking (tiling) is a basic approach
- Techniques apply generally, but the details (e.g., block size) are architecture dependent
- Similar techniques are possible on other data structures and algorithms
- Now it’s your turn: Homework 1 (due 9/19/02)
- Work in teams of 2 or 3 (assigned this time)

- CS267, Yelick