如何基于现有优化矩阵乘法函数进一步提升运算速度？

Nice work on the initial optimized implementation for GF(2) matrix multiplication (where * maps to bitwise AND and + maps to XOR) — let's break down concrete, actionable optimizations to push performance further, tailored to both CPU cache behavior and instruction-level parallelism.

1. Fix Loop Order for Cache Locality (Biggest Win!) #

Your current loop order (i → j → k) leads to strided memory access for matrix A, which kills cache hit rates. Since all matrices are stored in column-major order (element (row, col) is at row + col * n), we need to rearrange loops to prioritize contiguous memory access:

• For B[k + j*n]: When j is fixed, incrementing k accesses contiguous memory (entire j-th column of B).
• For A[i + k*n]: When k is fixed, incrementing i accesses contiguous memory (entire k-th column of A).
• For C[i + j*n]: When j is fixed, incrementing i accesses contiguous memory (entire j-th column of C).

Optimized loop order (j → k → i):

void matmul_optimized_cache(int n, int *A, int *B, int *C) {
    int i, j, k;
    int b_val;
    for (j = 0; j < n; ++j) {          // Traverse columns of C/B
        for (k = 0; k < n; ++k) {      // Traverse columns of A / rows of B
            b_val = B[k + j * n];      // Load B's (k,j) once (contiguous access)
            int *a_col = &A[k * n];    // Pointer to start of A's k-th column
            int *c_col = &C[j * n];    // Pointer to start of C's j-th column
            for (i = 0; i < n; ++i) {  // Contiguous access to A/C columns
                c_col[i] ^= a_col[i] & b_val;
            }
        }
    }
}

This change alone can yield 2-5x speedups depending on matrix size, as it eliminates most cache misses.

2. Loop Unrolling to Reduce Overhead #

Loop control (incrementing counters, checking conditions) adds overhead. Unroll the innermost loop to batch operations and leverage CPU instruction-level parallelism:

void matmul_optimized_unrolled(int n, int *A, int *B, int *C) {
    int i, j, k;
    int b_val;
    const int UNROLL_FACTOR = 4;
    for (j = 0; j < n; ++j) {
        for (k = 0; k < n; ++k) {
            b_val = B[k + j * n];
            int *a_col = &A[k * n];
            int *c_col = &C[j * n];
            // Unroll main loop
            for (i = 0; i < n - UNROLL_FACTOR + 1; i += UNROLL_FACTOR) {
                c_col[i] ^= a_col[i] & b_val;
                c_col[i+1] ^= a_col[i+1] & b_val;
                c_col[i+2] ^= a_col[i+2] & b_val;
                c_col[i+3] ^= a_col[i+3] & b_val;
            }
            // Handle remaining elements
            for (; i < n; ++i) {
                c_col[i] ^= a_col[i] & b_val;
            }
        }
    }
}

Adjust UNROLL_FACTOR based on your CPU (8 works well for wider CPUs).

3. SIMD Instructions for Parallel Bitwise Operations #

Modern CPUs support SIMD (Single Instruction, Multiple Data) which can process 8-16 32-bit integers per instruction. For x86/AVX2, use built-in functions to accelerate the bitwise operations:

#include <immintrin.h>

void matmul_optimized_simd(int n, int *A, int *B, int *C) {
    int i, j, k;
    const int SIMD_WIDTH = 8; // AVX2 handles 8x32-bit ints per instruction
    for (j = 0; j < n; ++j) {
        for (k = 0; k < n; ++k) {
            int b_val = B[k + j * n];
            // Broadcast b_val to all 8 lanes of a SIMD register
            __m256i b_vec = _mm256_set1_epi32(b_val);
            int *a_col = &A[k * n];
            int *c_col = &C[j * n];
            // Process SIMD-aligned blocks
            for (i = 0; i < n - SIMD_WIDTH + 1; i += SIMD_WIDTH) {
                __m256i a_vec = _mm256_loadu_si256((__m256i*)&a_col[i]);
                __m256i c_vec = _mm256_loadu_si256((__m256i*)&c_col[i]);
                __m256i and_result = _mm256_and_si256(a_vec, b_vec);
                __m256i xor_result = _mm256_xor_si256(c_vec, and_result);
                _mm256_storeu_si256((__m256i*)&c_col[i], xor_result);
            }
            // Handle leftover elements
            for (; i < n; ++i) {
                c_col[i] ^= a_col[i] & b_val;
            }
        }
    }
}

Pro Tip: Align Memory #

For maximum SIMD performance, allocate memory aligned to 32 bytes (AVX2's requirement):

// Replace malloc with aligned allocation
A = _mm_malloc(n * n * sizeof(int), 32);
B = _mm_malloc(n * n * sizeof(int), 32);
C1 = _mm_malloc(n * n * sizeof(int), 32);
C2 = _mm_malloc(n * n * sizeof(int), 32);

// Use _mm_free instead of free
_mm_free(A);
_mm_free(B);
_mm_free(C1);
_mm_free(C2);

Then replace _mm256_loadu_si256/_mm256_storeu_si256 with _mm256_load_si256/_mm256_store_si256.

4. Blocking (Tiling) for Multi-Level Cache Utilization #

Split the matrix into smaller tiles that fit entirely into L1/L2 cache, reducing trips to slower main memory. Combine this with the above optimizations for even better results:

#define TILE_SIZE 64 // Adjust based on your CPU's cache size (e.g., 32/64/128)

void matmul_optimized_tiled(int n, int *A, int *B, int *C) {
    int i, j, k, ii, jj, kk;
    // Traverse matrix in tiles
    for (jj = 0; jj < n; jj += TILE_SIZE) {
        for (kk = 0; kk < n; kk += TILE_SIZE) {
            for (ii = 0; ii < n; ii += TILE_SIZE) {
                // Process each tile
                for (j = jj; j < jj + TILE_SIZE && j < n; ++j) {
                    for (k = kk; k < kk + TILE_SIZE && k < n; ++k) {
                        int b_val = B[k + j * n];
                        int *a_col = &A[k * n];
                        int *c_col = &C[j * n];
                        for (i = ii; i < ii + TILE_SIZE && i < n; ++i) {
                            c_col[i] ^= a_col[i] & b_val;
                        }
                    }
                }
            }
        }
    }
}

5. Compiler Optimization Flags #

Don't forget to enable compiler optimizations to let the compiler handle instruction scheduling and additional micro-optimizations:

• GCC/Clang: -O3 -mavx2 -march=native
• MSVC: /O2 /arch:AVX2

Final Notes on Testing #

• Ensure your reference implementation (matmul_reference) matches the GF(2) logic (AND + XOR) to validate correctness.
• For large matrices (2048x2048), run multiple iterations (like your current REPEAT=5) to get stable cycle counts.

内容的提问来源于stack exchange，提问作者연승현