C与Fortran矩阵乘法性能差异显著的原因及优化咨询

我对比了Fortran与C语言的N×N方阵乘法性能，两者均采用循环分块优化内存访问，但Fortran性能显著更优：当N=1024时，Fortran耗时约0.18秒，C耗时约0.57秒，前者速度约为后者3倍。

测试环境与编译命令 #

• 系统：Windows 10
• CPU：Intel i7-6700
• 内存：16GB
• 编译器：MinGW 14.2

Fortran编译命令：

gfortran -O3 -march=native -funroll-all-loops matmul.f90

C编译命令：

gcc -O3 -march=native -funroll-all-loops matmul.c

测试代码 #

Fortran代码（matmul.f90） #

program matmul_fortran
    implicit none
    integer, parameter :: N = 1024   ! matrix size 
    integer, parameter :: BLOCK_SIZE = 32  ! block size 
    real(8), dimension(N, N) :: A, B, C
    integer :: i, j, k, i_block, j_block, k_block
    real(8) :: start, finish, temp

    ! initialize matrices A and B with random values
    call random_seed()
    call random_number(A)
    call random_number(B)
    C = 0.0  ! set matrix C to zero values

    call cpu_time(start)

    ! multiplication
    do i_block = 1, N, BLOCK_SIZE
        do j_block = 1, N, BLOCK_SIZE
            do k_block = 1, N, BLOCK_SIZE
                do i = i_block, min(i_block + BLOCK_SIZE - 1, N)
                    do j = j_block, min(j_block + BLOCK_SIZE - 1, N)
                        do k = k_block, min(k_block + BLOCK_SIZE - 1, N)
                            C(k, i) = C(k, i) + A(k, j)*B(j, i)
                        end do
                    end do
                end do
            end do
        end do
    end do

    call cpu_time(finish)

    print *, "Fortran Matrix Multiplication Time: ", finish - start, " seconds"

end program matmul_fortran

C代码（matmul.c） #

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 1024  // matrix size
#define BLOCK_SIZE 32  // block size

// function to initialize matrices with random values
void initialize_matrix(double *matrix) {
    for (int i = 0; i < N * N; i++) {
        matrix[i] = (double)rand() / RAND_MAX;  // Random values between 0 and 1
    }
}

int main() {
    double *A, *B, *C;
    clock_t start, end;

    A = (double *)malloc(N * N * sizeof(double));
    B = (double *)malloc(N * N * sizeof(double));
    C = (double *)malloc(N * N * sizeof(double));

    // set matrix C to zero values
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            C[i * N + j] = 0.0; 

    // Initialize matrices
    srand(time(NULL));
    initialize_matrix(A);
    initialize_matrix(B);

    start = clock();

    // multiplication
    for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE) {
        for (int j_block = 0; j_block < N; j_block += BLOCK_SIZE) {
            for (int k_block = 0; k_block < N; k_block += BLOCK_SIZE) {
                for (int i = i_block; i < i_block + BLOCK_SIZE && i < N; i++) {
                    for (int j = j_block; j < j_block + BLOCK_SIZE && j < N; j++) {
                        for (int k = k_block; k < k_block + BLOCK_SIZE && k < N; k++) {
                            C[i*N + k] += A[j*N + k]*B[i*N + j];
                        }
                    }
                }
            }
        }
    }
    
    end = clock();

    printf("C Matrix Multiplication Time: %.6f seconds\n", ((double)(end - start)) / CLOCKS_PER_SEC);

    free(A);
    free(B);
    free(C);

    return 0;
}

1. 内存布局不匹配：Fortran默认采用列优先存储，而C采用行优先。Fortran代码的访问模式完全契合列优先的连续内存读取；但C代码中A[j*N + k]是跨列的非连续访问，触发大量缓存失效，严重拖慢性能。
2. 循环顺序与缓存利用率：Fortran的循环顺序适配列优先的缓存命中逻辑，C代码的循环顺序未针对行优先优化，缓存命中率极低。

1. 调整循环顺序与内存访问模式 #

修改C代码的循环结构和索引，适配行优先存储，让内存访问变为连续模式：

// 优化后的核心乘法循环
for (int i_block = 0; i_block < N; i_block += BLOCK_SIZE) {
    for (int k_block = 0; k_block < N; k_block += BLOCK_SIZE) {
        for (int j_block = 0; j_block < N; j_block += BLOCK_SIZE) {
            for (int i = i_block; i < i_block + BLOCK_SIZE && i < N; i++) {
                for (int k = k_block; k < k_block + BLOCK_SIZE && k < N; k++) {
                    double temp = 0.0;
                    // 内层循环遍历连续内存
                    for (int j = j_block; j < j_block + BLOCK_SIZE && j < N; j++) {
                        temp += A[i*N + j] * B[j*N + k];
                    }
                    C[i*N + k] += temp;
                }
            }
        }
    }
}

• 调整循环块顺序为i→k→j，内层循环遍历j，保证A[i*N + j]是连续行访问。
• 引入临时变量减少对C的重复写操作，降低内存开销。

2. 使用对齐内存分配 #

malloc分配的内存不一定满足CPU缓存对齐要求，改用Windows下的_aligned_malloc分配64字节对齐内存（适配多数CPU缓存行大小）：

// 替换原malloc代码
A = (double*)_aligned_malloc(N*N*sizeof(double), 64);
B = (double*)_aligned_malloc(N*N*sizeof(double), 64);
C = (double*)_aligned_malloc(N*N*sizeof(double), 64);

// 释放内存时对应使用_aligned_free
_aligned_free(A);
_aligned_free(B);
_aligned_free(C);

3. 启用更激进的编译器优化 #

在编译命令中添加-ffast-math（允许数学优化，不影响矩阵乘法结果正确性）：

gcc -O3 -march=native -funroll-all-loops -ffast-math matmul.c

4. 手动SIMD向量化（可选） #

若编译器自动向量化效果不佳，可手动使用AVX指令优化内层循环：

#include <immintrin.h>

// 内层循环替换为AVX向量化代码
for (int j = j_block; j < j_block + BLOCK_SIZE; j += 4) {
    __m256d a_vec = _mm256_loadu_pd(&A[i*N + j]);
    __m256d b_vec = _mm256_loadu_pd(&B[j*N + k]);
    __m256d temp_vec = _mm256_add_pd(temp_vec, _mm256_mul_pd(a_vec, b_vec));
}
// 将向量结果累加后写入C
double temp[4];
_mm256_storeu_pd(temp, temp_vec);
C[i*N + k] += temp[0] + temp[1] + temp[2] + temp[3];

调整循环顺序和内存访问模式后，C代码性能会大幅提升，与Fortran的差距可缩小至10%-20%；配合对齐内存分配和-ffast-math，甚至可以追平Fortran的性能。

内容的提问来源于stack exchange，提问作者Ante Jurčević