#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 kernel32(long x, long y, long z, const double* __restrict__ A, const double* __restrict__ B, double* __restrict__ C) {
const long m = 32;
for (long k = 0; k < 2; ++k, z += 16) {
for (long j = 0; j < 2; ++j, y += 16) {
__m256d b00 = _mm256_load_pd(&B[idx(y + 0, z + 4 * 0)]);
__m256d b01 = _mm256_load_pd(&B[idx(y + 0, z + 4 * 1)]);
__m256d b02 = _mm256_load_pd(&B[idx(y + 0, z + 4 * 2)]);
__m256d b03 = _mm256_load_pd(&B[idx(y + 0, z + 4 * 3)]);
__m256d b10 = _mm256_load_pd(&B[idx(y + 1, z + 4 * 0)]);
__m256d b11 = _mm256_load_pd(&B[idx(y + 1, z + 4 * 1)]);
__m256d b12 = _mm256_load_pd(&B[idx(y + 1, z + 4 * 2)]);
__m256d b13 = _mm256_load_pd(&B[idx(y + 1, z + 4 * 3)]);
__m256d b20 = _mm256_load_pd(&B[idx(y + 2, z + 4 * 0)]);
__m256d b21 = _mm256_load_pd(&B[idx(y + 2, z + 4 * 1)]);
__m256d b22 = _mm256_load_pd(&B[idx(y + 2, z + 4 * 2)]);
__m256d b23 = _mm256_load_pd(&B[idx(y + 2, z + 4 * 3)]);
__m256d b30 = _mm256_load_pd(&B[idx(y + 3, z + 4 * 0)]);
__m256d b31 = _mm256_load_pd(&B[idx(y + 3, z + 4 * 1)]);
__m256d b32 = _mm256_load_pd(&B[idx(y + 3, z + 4 * 2)]);
__m256d b33 = _mm256_load_pd(&B[idx(y + 3, z + 4 * 3)]);
__m256d b40 = _mm256_load_pd(&B[idx(y + 4, z + 4 * 0)]);
__m256d b41 = _mm256_load_pd(&B[idx(y + 4, z + 4 * 1)]);
__m256d b42 = _mm256_load_pd(&B[idx(y + 4, z + 4 * 2)]);
__m256d b43 = _mm256_load_pd(&B[idx(y + 4, z + 4 * 3)]);
__m256d b50 = _mm256_load_pd(&B[idx(y + 5, z + 4 * 0)]);
__m256d b51 = _mm256_load_pd(&B[idx(y + 5, z + 4 * 1)]);
__m256d b52 = _mm256_load_pd(&B[idx(y + 5, z + 4 * 2)]);
__m256d b53 = _mm256_load_pd(&B[idx(y + 5, z + 4 * 3)]);
__m256d b60 = _mm256_load_pd(&B[idx(y + 6, z + 4 * 0)]);
__m256d b61 = _mm256_load_pd(&B[idx(y + 6, z + 4 * 1)]);
__m256d b62 = _mm256_load_pd(&B[idx(y + 6, z + 4 * 2)]);
__m256d b63 = _mm256_load_pd(&B[idx(y + 6, z + 4 * 3)]);
__m256d b70 = _mm256_load_pd(&B[idx(y + 7, z + 4 * 0)]);
__m256d b71 = _mm256_load_pd(&B[idx(y + 7, z + 4 * 1)]);
__m256d b72 = _mm256_load_pd(&B[idx(y + 7, z + 4 * 2)]);
__m256d b73 = _mm256_load_pd(&B[idx(y + 7, z + 4 * 3)]);
__m256d b80 = _mm256_load_pd(&B[idx(y + 8, z + 4 * 0)]);
__m256d b81 = _mm256_load_pd(&B[idx(y + 8, z + 4 * 1)]);
__m256d b82 = _mm256_load_pd(&B[idx(y + 8, z + 4 * 2)]);
__m256d b83 = _mm256_load_pd(&B[idx(y + 8, z + 4 * 3)]);
__m256d b90 = _mm256_load_pd(&B[idx(y + 9, z + 4 * 0)]);
__m256d b91 = _mm256_load_pd(&B[idx(y + 9, z + 4 * 1)]);
__m256d b92 = _mm256_load_pd(&B[idx(y + 9, z + 4 * 2)]);
__m256d b93 = _mm256_load_pd(&B[idx(y + 9, z + 4 * 3)]);
__m256d ba0 = _mm256_load_pd(&B[idx(y + 10, z + 4 * 0)]);
__m256d ba1 = _mm256_load_pd(&B[idx(y + 10, z + 4 * 1)]);
__m256d ba2 = _mm256_load_pd(&B[idx(y + 10, z + 4 * 2)]);
__m256d ba3 = _mm256_load_pd(&B[idx(y + 10, z + 4 * 3)]);
__m256d bb0 = _mm256_load_pd(&B[idx(y + 11, z + 4 * 0)]);
__m256d bb1 = _mm256_load_pd(&B[idx(y + 11, z + 4 * 1)]);
__m256d bb2 = _mm256_load_pd(&B[idx(y + 11, z + 4 * 2)]);
__m256d bb3 = _mm256_load_pd(&B[idx(y + 11, z + 4 * 3)]);
__m256d bc0 = _mm256_load_pd(&B[idx(y + 12, z + 4 * 0)]);
__m256d bc1 = _mm256_load_pd(&B[idx(y + 12, z + 4 * 1)]);
__m256d bc2 = _mm256_load_pd(&B[idx(y + 12, z + 4 * 2)]);
__m256d bc3 = _mm256_load_pd(&B[idx(y + 12, z + 4 * 3)]);
__m256d bd0 = _mm256_load_pd(&B[idx(y + 13, z + 4 * 0)]);
__m256d bd1 = _mm256_load_pd(&B[idx(y + 13, z + 4 * 1)]);
__m256d bd2 = _mm256_load_pd(&B[idx(y + 13, z + 4 * 2)]);
__m256d bd3 = _mm256_load_pd(&B[idx(y + 13, z + 4 * 3)]);
__m256d be0 = _mm256_load_pd(&B[idx(y + 14, z + 4 * 0)]);
__m256d be1 = _mm256_load_pd(&B[idx(y + 14, z + 4 * 1)]);
__m256d be2 = _mm256_load_pd(&B[idx(y + 14, z + 4 * 2)]);
__m256d be3 = _mm256_load_pd(&B[idx(y + 14, z + 4 * 3)]);
__m256d bf0 = _mm256_load_pd(&B[idx(y + 15, z + 4 * 0)]);
__m256d bf1 = _mm256_load_pd(&B[idx(y + 15, z + 4 * 1)]);
__m256d bf2 = _mm256_load_pd(&B[idx(y + 15, z + 4 * 2)]);
__m256d bf3 = _mm256_load_pd(&B[idx(y + 15, 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 + 0)]);
c0 += a * b00;
c1 += a * b01;
c2 += a * b02;
c3 += a * b03;
a = _mm256_set1_pd(A[idx(x + i, y + 1)]);
c0 += a * b10;
c1 += a * b11;
c2 += a * b12;
c3 += a * b13;
a = _mm256_set1_pd(A[idx(x + i, y + 2)]);
c0 += a * b20;
c1 += a * b21;
c2 += a * b22;
c3 += a * b23;
a = _mm256_set1_pd(A[idx(x + i, y + 3)]);
c0 += a * b30;
c1 += a * b31;
c2 += a * b32;
c3 += a * b33;
a = _mm256_set1_pd(A[idx(x + i, y + 4)]);
c0 += a * b40;
c1 += a * b41;
c2 += a * b42;
c3 += a * b43;
a = _mm256_set1_pd(A[idx(x + i, y + 5)]);
c0 += a * b50;
c1 += a * b51;
c2 += a * b52;
c3 += a * b53;
a = _mm256_set1_pd(A[idx(x + i, y + 6)]);
c0 += a * b60;
c1 += a * b61;
c2 += a * b62;
c3 += a * b63;
a = _mm256_set1_pd(A[idx(x + i, y + 7)]);
c0 += a * b70;
c1 += a * b71;
c2 += a * b72;
c3 += a * b73;
a = _mm256_set1_pd(A[idx(x + i, y + 8)]);
c0 += a * b80;
c1 += a * b81;
c2 += a * b82;
c3 += a * b83;
a = _mm256_set1_pd(A[idx(x + i, y + 9)]);
c0 += a * b90;
c1 += a * b91;
c2 += a * b92;
c3 += a * b93;
a = _mm256_set1_pd(A[idx(x + i, y + 10)]);
c0 += a * ba0;
c1 += a * ba1;
c2 += a * ba2;
c3 += a * ba3;
a = _mm256_set1_pd(A[idx(x + i, y + 11)]);
c0 += a * bb0;
c1 += a * bb1;
c2 += a * bb2;
c3 += a * bb3;
a = _mm256_set1_pd(A[idx(x + i, y + 12)]);
c0 += a * bc0;
c1 += a * bc1;
c2 += a * bc2;
c3 += a * bc3;
a = _mm256_set1_pd(A[idx(x + i, y + 13)]);
c0 += a * bd0;
c1 += a * bd1;
c2 += a * bd2;
c3 += a * bd3;
a = _mm256_set1_pd(A[idx(x + i, y + 14)]);
c0 += a * be0;
c1 += a * be1;
c2 += a * be2;
c3 += a * be3;
a = _mm256_set1_pd(A[idx(x + i, y + 15)]);
c0 += a * bf0;
c1 += a * bf1;
c2 += a * bf2;
c3 += a * bf3;
_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);
}
}
y -= 32;
}
}
static void gao(int s, int x, int y, int 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 == 5) {
kernel32(x, y, z, A, B, C);
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);
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 110.421 ms | 8 MB + 20 KB | Accepted | Score: 100 | 显示更多 |