代码

#define idx(i, j) ((i) * n + (j))

void matrix_multiply(int n, const double *A, const double *B, double *C) {
	const int BLOCK = 2;
	for (int ii = 0; ii < n; ii += BLOCK)
		for (int jj = 0; jj < n; jj += BLOCK)
			for (int i = ii; i < ii + BLOCK; i ++)
				for (int j = jj; j < jj + BLOCK; j ++)
				{
					double a_ij = A[idx(i, j)];
					register double *c_i = C + i * n;
                                        register const double *b_j = B + j * n;
					const int UNROLL = 4;
					for (int k = 0; k < n; k += UNROLL)
					{
						c_i[0] += a_ij * b_j[0];
						c_i[1] += a_ij * b_j[1];
						c_i[2] += a_ij * b_j[2];
						c_i[3] += a_ij * b_j[3];
                                                c_i += UNROLL;
                                                b_j += UNROLL;
					}
				}
}

评测结果