I need help to make matrix multiplication in C run as fast as possible. On my AMD Phenom(tm) II X6 1090T, my program multiplies two square singe precision 4096x4096 matrices in about 6.9 seconds. NumPy, which is based on OpenBLAS, multiplies the same sized matrices in about 1.5 seconds. So I think it should be possible to double the speed of my program. Since my processor is old, only 128-bit SIMD instructions are available.
Original program
#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define SIZE 4096
#define N_ELS (SIZE * SIZE)
#define N_BYTES (sizeof(float) * N_ELS)
#define MIN(a, b) ((a > b) ? (b) : (a))
#define N_THREADS 6
#define TILE_SIZE 32
static void
mul_tile(int i0, int i1, int j0, int j1, int k0, int k1,
float * restrict A, float * restrict B, float * restrict C) {
for (int i = i0; i < i1; i++) {
for (int k = k0; k < k1; k++) {
float a = A[SIZE * i + k];
__builtin_prefetch(&C[SIZE * i + j0]);
__builtin_prefetch(&B[SIZE * (k + 1) + j0]);
__builtin_prefetch(&B[SIZE * (k + 2) + j0]);
for (int j = j0; j < j1; j++) {
float b = B[SIZE * k + j];
C[SIZE * i + j] += a * b;
}
}
}
}
typedef struct {
float *A, *B, *C;
int start_i, end_i;
} mul_job_t;
void *mul_thread(void *arg) {
mul_job_t *job = (mul_job_t *)arg;
int start_i = job->start_i;
int end_i = job->end_i;
for (int i = start_i; i < end_i; i += TILE_SIZE) {
int imax = MIN(i + TILE_SIZE, SIZE);
for (int j = 0; j < SIZE; j += TILE_SIZE) {
int jmax = MIN(j + TILE_SIZE, SIZE);
for (int k = 0; k < SIZE; k += TILE_SIZE) {
int kmax = MIN(k + TILE_SIZE, SIZE);
mul_tile(i, imax, j, jmax, k, kmax,
job->A, job->B, job->C);
}
}
}
return 0;
}
static void
mul(float * restrict A, float * restrict B, float * restrict C) {
pthread_t threads[N_THREADS];
mul_job_t jobs[N_THREADS];
memset(C, 0, N_BYTES);
int n_i_tiles = (int)ceil((float)SIZE / (float)TILE_SIZE);
int tiles_per_thread = (int)ceil((float)n_i_tiles / (float)N_THREADS);
for (int i = 0; i < N_THREADS; i++) {
int start = TILE_SIZE * i * tiles_per_thread;
int end = MIN(start + TILE_SIZE * tiles_per_thread, SIZE);
jobs[i] = (mul_job_t){A, B, C, start, end};
pthread_create(&threads[i], NULL, mul_thread, &jobs[i]);
}
for (int i = 0; i < N_THREADS; i++) {
pthread_join(threads[i], NULL);
}
}
int
main(int argc, char *argv[]) {
float *A = (float *)malloc(N_BYTES);
float *B = (float *)malloc(N_BYTES);
float *C = (float *)malloc(N_BYTES);
for (int i = 0; i < SIZE * SIZE; i++) {
A[i] = ((float)rand() / (float)RAND_MAX) * 10;
B[i] = ((float)rand() / (float)RAND_MAX) * 10;
}
struct timespec begin, end;
clock_gettime(CLOCK_MONOTONIC_RAW, &begin);
mul(A, B, C);
clock_gettime(CLOCK_MONOTONIC_RAW, &end);
double delta = (end.tv_nsec - begin.tv_nsec) / 1000000000.0 +
(end.tv_sec - begin.tv_sec);
printf("%.6lfs\n", delta);
free(A);
free(B);
free(C);
}
For best results, compile with clang -o mul mul.c -O3 -fomit-frame-pointer -march=native -mtune=native -lpthread.
Here is my cache hierarchy:
Any answer that significantly improves the performance is acceptable to me. I'm only interested in performance and not other issues the code might have.
Notes and additions:
- gcc 12.2 is about 0.4 seconds slower than clang 14.0 at 7.2 seconds.
- The machine has a Radeon HD 4200 card, so if someone could optimize the code using GPU instructions that would be cool too.
- The code is veeery cache sensitive; setting
TILE_SIZEto 64 quadruples the execution time, for example.
Optimized program
#include <math.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define SIZE 4096
#define N_ELS (SIZE * SIZE)
#define N_BYTES (sizeof(float) * N_ELS)
#define MIN(a, b) ((a > b) ? (b) : (a))
#define N_THREADS 6
#define TILE_SIZE 32
typedef unsigned int uint_t;
static void
mul_tile(uint_t i0, uint_t i1, uint_t j0, uint_t j1, uint_t k0, uint_t k1,
float * restrict A, float * restrict B, float * restrict C
) {
for (uint_t i = i0; i < i1; i ++) {
for (uint_t k = k0; k < k1; k++) {
float a = A[SIZE * i + k];
for (uint_t j = j0; j < j1; j++) {
float b = B[SIZE * k + j];
C[SIZE * i + j] += a * b;
}
}
}
}
typedef struct {
float *A, *B, *C;
uint_t start_i, end_i;
} mul_job_t;
void *mul_thread(void *arg) {
mul_job_t *job = (mul_job_t *)arg;
uint_t start_i = job->start_i;
uint_t end_i = job->end_i;
for (uint_t i = start_i; i < end_i; i += TILE_SIZE) {
uint_t imax = MIN(i + TILE_SIZE, SIZE);
for (uint_t j = 0; j < SIZE; j += TILE_SIZE) {
uint_t jmax = MIN(j + TILE_SIZE, SIZE);
for (uint_t k = 0; k < SIZE; k += TILE_SIZE) {
uint_t kmax = MIN(k + TILE_SIZE, SIZE);
mul_tile(i, imax, j, jmax, k, kmax,
job->A, job->B, job->C);
}
}
}
return 0;
}
static void
mul(float * restrict A, float * restrict B, float * restrict C) {
pthread_t threads[N_THREADS];
mul_job_t jobs[N_THREADS];
memset(C, 0, N_BYTES);
int n_i_tiles = (int)ceil((float)SIZE / (float)TILE_SIZE);
int tiles_per_thread = (int)ceil((float)n_i_tiles / (float)N_THREADS);
for (int i = 0; i < N_THREADS; i++) {
int start = TILE_SIZE * i * tiles_per_thread;
int end = MIN(start + TILE_SIZE * tiles_per_thread, SIZE);
jobs[i] = (mul_job_t){A, B, C, start, end};
pthread_create(&threads[i], NULL, mul_thread, &jobs[i]);
}
for (int i = 0; i < N_THREADS; i++) {
pthread_join(threads[i], NULL);
}
}
int
main(int argc, char *argv[]) {
float *A = (float *)malloc(N_BYTES);
float *B = (float *)malloc(N_BYTES);
float *C = (float *)malloc(N_BYTES);
for (int i = 0; i < SIZE * SIZE; i++) {
A[i] = ((float)rand() / (float)RAND_MAX) * 10;
B[i] = ((float)rand() / (float)RAND_MAX) * 10;
}
struct timespec begin, end;
clock_gettime(CLOCK_MONOTONIC_RAW, &begin);
mul(A, B, C);
clock_gettime(CLOCK_MONOTONIC_RAW, &end);
double delta = (end.tv_nsec - begin.tv_nsec) / 1000000000.0 +
(end.tv_sec - begin.tv_sec);
printf("%.6lfs\n", delta);
free(A);
free(B);
free(C);
}
