代码

// Naive-opt-2

#pragma GCC target("avx2,fma")

#include <string.h>
#include <x86intrin.h>

// dst += src * k
static inline void vector_mul_add(int n, double *dst, const double *src, double k) {
	__m256d K = _mm256_set1_pd(k);
	
	#define DST(i) (* (__m256d *) (dst + (i)))
	#define SRC(i) (* (__m256d *) (src + (i)))
	
	// assume 4096-aligned
	for (int i = 0; i < n; i += 16) {
		DST(i + 4 * 0) = _mm256_fmadd_pd(SRC(i + 4 * 0), K, DST(i + 4 * 0));
		DST(i + 4 * 1) = _mm256_fmadd_pd(SRC(i + 4 * 1), K, DST(i + 4 * 1));
		DST(i + 4 * 2) = _mm256_fmadd_pd(SRC(i + 4 * 2), K, DST(i + 4 * 2));
		DST(i + 4 * 3) = _mm256_fmadd_pd(SRC(i + 4 * 3), K, DST(i + 4 * 3));
	}
}

void matrix_multiply(int n, const double *A, const double *B, double *C) {
	for (int i = 0; i < n; i++) {
		const double *ai = A + i * n;
		double *ci = C + i * n;
		memset(ci, 0, n * sizeof(double));
		
		for (int k = 0; k < n; k++) {
			const double *bk = B + k * n;
			const double aik = ai[k];
			
			vector_mul_add(n, ci, bk, aik);
		}
	}
}

评测结果