I have been playing around with Creel's video on optimising matrix multiplicationn and I don't get the speedups he does. What is the reason for this? Below is the program I used to benchmark. There are three functions: naive multiplication, in-place transpose of B, and in-place transpose of B + blocking. I ran this with n = 4000 and block sizes 1, 10, 20, 50, 100, 200. My caches are 32KB L1D, 256KB L2, 4MB L3 shared, so block size 10 should be 20 * 20 * 8 * 2 = 6.4KB, and fit comfortably in L1 cache. No matter the block size, it takes 50s, which is the same as for only transposing. I compiled with gcc -O3 -mavx2.
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
void matmul(size_t n, double A[n][n], double B[n][n], double result[n][n])
{
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < n; j++) {
double acc = 0;
for (size_t k = 0; k < n; k++) {
acc += A[i][k] * B[k][j];
}
result[i][j] = acc;
}
}
}
void transpose(size_t n, double matrix[n][n])
{
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < i; j++) {
double temp = matrix[i][j];
matrix[i][j] = matrix[j][i];
matrix[j][i] = temp;
}
}
}
void matmulTrans(size_t n, double A[n][n], double B[n][n], double result[n][n])
{
transpose(n, B);
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < n; j++) {
double acc = 0;
for (size_t k = 0; k < n; k++) {
acc += A[i][k] * B[j][k];
}
result[i][j] = acc;
}
}
}
void matmulBlock(size_t n, double A[n][n], double B[n][n],
double result[n][n], size_t blockSize)
{
transpose(n, B);
for (size_t i = 0; i < n; i += blockSize) {
for (size_t j = 0; j < n; j += blockSize) {
for (size_t iBlock = i; iBlock < i + blockSize; iBlock++) {
for (size_t jBlock = j; jBlock < j + blockSize; jBlock++) {
double acc = 0;
for (size_t k = 0; k < n; k++) {
acc += A[iBlock][k] * B[jBlock][k];
}
result[iBlock][jBlock] = acc;
}
}
}
}
}
int main(int argc, char **argv)
{
if (argc != 3) {
printf("Provide two arguments!\n");
return 1;
}
int n = atoi(argv[1]);
int blockSize = atoi(argv[2]);
double (*A)[n] = malloc(n * n * sizeof(double));
double (*B)[n] = malloc(n * n * sizeof(double));
double (*result)[n] = malloc(n * n * sizeof(double));
clock_t time1 = clock();
matmulBlock(n, A, B, result, blockSize);
clock_t time2 = clock();
// matmul(n, A, B, result);
clock_t time3 = clock();
matmulTrans(n, A, B, result);
clock_t time4 = clock();
printf("Blocked version: %lfs.\nNaive version: %lfs.\n"
"Transposed version: %lfs.\n",
(double) (time2 - time1) / CLOCKS_PER_SEC,
(double) (time3 - time2) / CLOCKS_PER_SEC,
(double) (time4 - time3) / CLOCKS_PER_SEC);
free(A);
free(B);
free(result);
return 0;
}
The problem is that I had only blocked the i and the j loop. This means that we essentially block A into an blockSize x 1
matrix of (n / blockSize) x n
blocks and B into an 1 x blockSize
matrix of n x (n / blockSize)
blocks. These blocks are far too large to fit in cache. Using
void matmulBlock(size_t n, double A[n][n], double B[n][n],
double result[__restrict__ n][n], size_t blockSize)
{
for (size_t i = 0; i < n; i += blockSize) {
for (size_t j = 0; j < n; j += blockSize) {
for (size_t k = 0; k < n; k += blockSize) {
for (size_t iBlock = i; iBlock < i + blockSize; iBlock++) {
for (size_t kBlock = k; kBlock < k + blockSize; kBlock++) {
for (size_t jBlock = j; jBlock < j + blockSize; jBlock++) {
result[iBlock][jBlock] += A[iBlock][kBlock] * B[kBlock][jBlock];
}
}
}
}
}
}
}
instead does lead to speedups.