Advanced Linear Algebra: Foundations to Frontiers

Subsection 12.1.1 Simple Implementation of matrix-matrix multiplication

The current coronavirus crisis hit UT-Austin on March 14, 2020, a day we spent quickly making videos for Week 11. We have not been back to the office, to create videos for Week 12, since then. We will likely add such videos as time goes on. For now, we hope that the notes suffice.

Let \(A \text{,}\) \(B \text{,}\) and \(C \) be \(m \times k \text{,}\) \(k \times n \text{,}\) and \(m \times n \) matrices, respectively. We can expose their individual entries as

and

The computation \(C := A B + C \text{,}\) which adds the result of the matrix-matrix multiplication \(A B \) to a matrix \(C \text{,}\) is defined element-wise as

for all \(0 \leq i \lt m \) and \(0 \leq j \lt n \text{.}\) We add to \(C \) because this will make it easier to play with the orderings of the loops when implementing matrix-matrix multiplication. The following pseudo-code computes \(C := A B + C \text{:}\)

The outer two loops visit each element of \(C \) and the inner loop updates \(\gamma_{i,j} \) with (12.1.1). We use C programming language macro definitions in order to explicitly index into the matrices, which are passed as one-dimentional arrays in which the matrices are stored in column-major order.

#define alpha( i,j ) A[ (j)*ldA + i ]   // map alpha( i,j ) to array A 
#define beta( i,j )  B[ (j)*ldB + i ]   // map beta( i,j )  to array B
#define gamma( i,j ) C[ (j)*ldC + i ]   // map gamma( i,j ) to array C

void MyGemm( int m, int n, int k, double *A, int ldA,
	     double *B, int ldB, double *C, int ldC )
{
  for ( int i=0; i<m; i++ )
    for ( int j=0; j<n; j++ )
      for ( int p=0; p<k; p++ )
        gamma( i,j ) += alpha( i,p ) * beta( p,j );
}

make IJP

On Robert's laptop, Homework 12.1.1.1 yields the graph

The performance of a matrix-matrix multiplication implementation is measured in billions of floating point operations per second (GFLOPS). We know that it takes \(2 m n k \) flops to compute \(C := A B + C \) when \(C \) is \(m \times n \text{,}\) \(A\) is \(m \times k \text{,}\) and \(B \) is \(k \times n \text{.}\) If we measure the time it takes to complete the computation, \(T( m, n, k ) \text{,}\) then the rate at which we compute is given by

For our implementation, this yields

make IPJ
make JIP
   ...

make All_Orderings

Homework 12.1.1.4.

In directory Assignments/Week12/C/ execute

make JPI

and view the results with the Live Script in Assignments/Week12/C/data/Plot_Opener.mlx. (This may take a little while, since the Makefile now specifies that the largest problem to be executed is \(m=n=k=1500\text{.}\))

Next, change that Live Script to also show the performance of the reference implementation provided by the BLAS-like Library Instantion Software (BLIS): Change

% Optionally show the reference implementation performance data 
if ( 0 )

to

% Optionally show the reference implementation performance data 
if ( 1 )

and rerun the Live Script. This adds a plot to the graph for the reference implementation.

What do you observe? Now are you happy with the improvements you made by reordering the loops>

Solution

On Robert's laptop:

Left: Plotting only simple implementations. Right: Adding the performance of the reference implementation provided by BLIS.

Note: the performance in the graph on the left may not exactly match that in the graph earlier in this unit. My laptop does not always attain the same performance. When a processor gets hot, it "clocks down." This means the attainable performance goes down. A laptop is not easy to cool, so one would expect more fluctuation than when using, for example, a desktop or a server.

Here is a video from our course "LAFF-On Programming for High Performance", which explains what you observe. (It refers to "Week 1" or that course. It is part of the launch for that course.)