DGEMM: Matrix Multiplication

1Learning Outcomes¶

• Explain the DGEMM benchmark: row-major order matrix multiplication.

• Explain why C DGEMM runs faster than Python DGEMM.

• Understand that in practice, library implementations like NumPy DGEMM are plenty fast because they use the optimizations described in this chapter.

How might we begin evaluating performance? Remember from our discussion of the iron law of processor performance that in order to evaluate performance, we must determine a program benchmark. Over the next few lectures, we will evaluate different optimizations of a core benchmark to many engineering, data, and image processing tasks today: matrix multiplication.

2Matrix Multiplication¶

Matrix multiplication is defined as $C = AB$ for matrices $A$, $B$, and $C$. If $\mathbb{R}$ is the set of all real numbers, the dimensions of these three matrices are: $A \in \mathbb{R}^{n \times d}, B \in \mathbb{R}^{d \times m}, C \in \mathbb{R}^{n \times m}$.

To compute each element of the resulting matrix $C$, we take the dot product of a row of $A$ and a column of $B$. The dot product of two vectors is the sum of the element-wise product of the two vectors. Figure 2 shows how we can compute the zero-th row, zero-th column element of $C$, $C_{00}$, by multiplying element-wise the (zero-indexed) zero-th row of $A$ and zero-th column of $B$, then summing everything together.

Similarly, to compute $C_{01}$, we can multiply element-wise the zero-th row of $A$ and first column of $B$, then sum everything together, as in Figure 3.

// basic matrix multiplication
void dgemm(int N, double **A, double **B, double **C) {
  for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
        // for row i, col j of C
        double sum = 0;
        for (int k = 0; k < size; k++) {
        sum += A[i][k] * B[k][j];
        }
        C[i][j] = sum;
    }
  }
}

The code above assumes that matrices are stored as two-dimensional arrays, but this will complicate our analysis of memory—particularly as it relates to spatial locality for caches and the memory hierarchy. We therefore describe the row-major order implementation below, which we will use as our benchmark task.

3The DGEMM Code Benchmark¶

The above matrix multiplication task is common called DGEMM, which stands for Double-precision GEneral Matrix Multiply.

The following code represents matrices $A$, $B$, and $C$ as variables A, B, and C, respectively. Further assumptions:

• Heap allocated;

• Square with dimensions $n \times n$; and

• Have elements stored in a double array, in row-major order.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
matrix_d_t *dgemm(matrix_d_t *A, matrix_d_t *B) {
  if (A->ncols != B->nrows) return NULL;
  matrix_d_t *C = init_mat_d(A->nrows, B->ncols);
    
  for (int i = 0; i < A->nrows; i++) {
    for (int j = 0; j < B->ncols; j++) {
      double sum = 0; 
      for (int k = 0; k < A->ncols; k++) {
        sum += A->data[i*A->ncols+k]*B->data[k*B->ncols+j];
      }
      C->data[i*C->ncols+j] = sum;
    } 
  }   
  return C;
}

Figure 4 shows the innermost loop of the DGEMM code to compute $C_{ij}$, where $i = j = 0$.

3.1Matrix Dimensions¶

The P&H textbook assumes that all matrices are square, where N specifies both row and column dimensions. For those less familiar with matrix multiplication, reducing to this one value may make it hard to picture matrix multiplication, even if it makes our code simpler.

For the purposes of our course demo, we will assume DGEMM works on a struct that stores the row and column dimensions along with our matrix data. as shown below. All matrices are pointers to heap-allocated structs; these are allocated and freed with init_mat_d and free_mat_d.

1
2
3
4
5
6
7
8
9
10
typedef struct {
    int nrows; 
    int ncols;
  double* data;
} matrix_d_t;
            
matrix_d_t *init_mat_d(int nrows, int ncols);
matrix_d_t *load_mat_d(FILE* f);
void print_mat_d(matrix_d_t *mat);
void free_mat_d(matrix_d_t *mat);

Going forward, we will use this syntax, but all of our timing benchmarks and subsequent optimizations will assume that all matrices are square, i.e., $n = d = m$. We will further assume $n$ powers of two (or at minimum a multiple of 4). We will test:

• $512 \times 512$

• $1024 \times 1024$

• $2048 \times 2048$

• etc.

In practice, matrix dimensions are passed as parameters to avoid structs and keep indirection to a minimum. We will probably revisit and rewrite this benchmark in future semesters to measure performance more clearly.

3.2Row-major order¶

Assume that matrices $A$, $B$, and $C$ are stored as A, B, and C. In the code, these matrices are stoerd as arrays of double arrays. By convention, these arrays store matrix elements in row-major order.^[1]

In Figure 5, each element of A[i] is the $i$-th row of matrix $A$; furthermore, the zero-th element of the $i+1$-th row immediately follows the last element of the $i$-th row, as shown in Figure 5.

We revisit the inner loop of the DGEMM benchmark code. Flip between the two tabs to see how to compute $C_{ij}$, the element of $C$ in row i and column j if we assume all matrices are stored in row-major order.

double sum = 0; 
for (int k = 0; k < A->ncols; k++) {
  sum += A->data[i*A->ncols+k]*B->data[k*B->ncols+j];
}
C->data[i*C->ncols+j] = sum;

double sum = 0;
for (int k = 0; k < size; k++) {
sum += A[i][k] * B[k][j];
}
C[i][j] = sum;