代码

#pragma GCC optimize("Ofast,no-stack-protector")
#pragma GCC target("avx2,fma")
#include <string.h>
#include <x86intrin.h>

#define n 1024
#define idx(i, j) (((i) << 10) + (j))

static void gao(int s, long x, long y, long z, int dx, int dy, int dz, int dx2, int dy2, int dz2, int dx3, int dy3, int dz3, const double* __restrict__ A, const double* __restrict__ B, double* __restrict__ C) {
    if (s == 4) {
        const long m = 16;
        for (long j = 0; j < m; j += 4) {
            __m256d b00 = _mm256_load_pd(&B[idx(y + j + 0, z + 4 * 0)]);
            __m256d b01 = _mm256_load_pd(&B[idx(y + j + 0, z + 4 * 1)]);
            __m256d b02 = _mm256_load_pd(&B[idx(y + j + 0, z + 4 * 2)]);
            __m256d b03 = _mm256_load_pd(&B[idx(y + j + 0, z + 4 * 3)]);
            __m256d b10 = _mm256_load_pd(&B[idx(y + j + 1, z + 4 * 0)]);
            __m256d b11 = _mm256_load_pd(&B[idx(y + j + 1, z + 4 * 1)]);
            __m256d b12 = _mm256_load_pd(&B[idx(y + j + 1, z + 4 * 2)]);
            __m256d b13 = _mm256_load_pd(&B[idx(y + j + 1, z + 4 * 3)]);
            __m256d b20 = _mm256_load_pd(&B[idx(y + j + 2, z + 4 * 0)]);
            __m256d b21 = _mm256_load_pd(&B[idx(y + j + 2, z + 4 * 1)]);
            __m256d b22 = _mm256_load_pd(&B[idx(y + j + 2, z + 4 * 2)]);
            __m256d b23 = _mm256_load_pd(&B[idx(y + j + 2, z + 4 * 3)]);
            __m256d b30 = _mm256_load_pd(&B[idx(y + j + 3, z + 4 * 0)]);
            __m256d b31 = _mm256_load_pd(&B[idx(y + j + 3, z + 4 * 1)]);
            __m256d b32 = _mm256_load_pd(&B[idx(y + j + 3, z + 4 * 2)]);
            __m256d b33 = _mm256_load_pd(&B[idx(y + j + 3, z + 4 * 3)]);
            for (long i = 0; i < m; ++i) {
                __m256d c0 = _mm256_load_pd(&C[idx(x + i, z + 4 * 0)]);
                __m256d c1 = _mm256_load_pd(&C[idx(x + i, z + 4 * 1)]);
                __m256d c2 = _mm256_load_pd(&C[idx(x + i, z + 4 * 2)]);
                __m256d c3 = _mm256_load_pd(&C[idx(x + i, z + 4 * 3)]);
                __m256d a;
                a = _mm256_set1_pd(A[idx(x + i, y + j + 0)]);
                c0 = _mm256_fmadd_pd(a, b00, c0);
                c1 = _mm256_fmadd_pd(a, b01, c1);
                c2 = _mm256_fmadd_pd(a, b02, c2);
                c3 = _mm256_fmadd_pd(a, b03, c3);
                a = _mm256_set1_pd(A[idx(x + i, y + j + 1)]);
                c0 = _mm256_fmadd_pd(a, b10, c0);
                c1 = _mm256_fmadd_pd(a, b11, c1);
                c2 = _mm256_fmadd_pd(a, b12, c2);
                c3 = _mm256_fmadd_pd(a, b13, c3);
                a = _mm256_set1_pd(A[idx(x + i, y + j + 2)]);
                c0 = _mm256_fmadd_pd(a, b20, c0);
                c1 = _mm256_fmadd_pd(a, b21, c1);
                c2 = _mm256_fmadd_pd(a, b22, c2);
                c3 = _mm256_fmadd_pd(a, b23, c3);
                a = _mm256_set1_pd(A[idx(x + i, y + j + 3)]);
                c0 = _mm256_fmadd_pd(a, b30, c0);
                c1 = _mm256_fmadd_pd(a, b31, c1);
                c2 = _mm256_fmadd_pd(a, b32, c2);
                c3 = _mm256_fmadd_pd(a, b33, c3);
                _mm256_store_pd(&C[idx(x + i, z + 4 * 0)], c0);
                _mm256_store_pd(&C[idx(x + i, z + 4 * 1)], c1);
                _mm256_store_pd(&C[idx(x + i, z + 4 * 2)], c2);
                _mm256_store_pd(&C[idx(x + i, z + 4 * 3)], c3);
            }
        }
        return;
    }
    --s;
    if (dx < 0) x -= dx << s;
    if (dy < 0) y -= dy << s;
    if (dz < 0) z -= dz << s;
    if (dx2 < 0) x -= dx2 << s;
    if (dy2 < 0) y -= dy2 << s;
    if (dz2 < 0) z -= dz2 << s;
    if (dx3 < 0) x -= dx3 << s;
    if (dy3 < 0) y -= dy3 << s;
    if (dz3 < 0) z -= dz3 << s;
    gao(s, x, y, z, dx2, dy2, dz2, dx3, dy3, dz3, dx, dy, dz, A, B, C);
    gao(s, x + (dx << s), y + (dy << s), z + (dz << s), dx3, dy3, dz3, dx, dy, dz, dx2, dy2, dz2, A, B, C);
    gao(s, x + (dx << s) + (dx2 << s), y + (dy << s) + (dy2 << s), z + (dz << s) + (dz2 << s), dx3, dy3, dz3, dx, dy, dz, dx2, dy2, dz2, A, B, C);
    gao(s, x + (dx2 << s), y + (dy2 << s), z + (dz2 << s), -dx, -dy, -dz, -dx2, -dy2, -dz2, dx3, dy3, dz3, A, B, C);
    gao(s, x + (dx2 << s) + (dx3 << s), y + (dy2 << s) + (dy3 << s), z + (dz2 << s) + (dz3 << s), -dx, -dy, -dz, -dx2, -dy2, -dz2, dx3, dy3, dz3, A, B, C);
    gao(s, x + (dx << s) + (dx2 << s) + (dx3 << s), y + (dy << s) + (dy2 << s) + (dy3 << s), z + (dz << s) + (dz2 << s) + (dz3 << s), -dx3, -dy3, -dz3, dx, dy, dz, -dx2, -dy2, -dz2, A, B, C);
    gao(s, x + (dx << s) + (dx3 << s), y + (dy << s) + (dy3 << s), z + (dz << s) + (dz3 << s), -dx3, -dy3, -dz3, dx, dy, dz, -dx2, -dy2, -dz2, A, B, C);
    gao(s, x + (dx3 << s), y + (dy3 << s), z + (dz3 << s), dx2, dy2, dz2, -dx3, -dy3, -dz3, -dx, -dy, -dz, A, B, C);
}
void matrix_multiply(int, const double* A, const double* B, double* C) {
    memset(C, 0, 1024 * 1024 * sizeof(double));
    gao(10, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, A, B, C);
}

评测结果