优化一维存储对称带状矩阵与通用二维矩阵的加法运算

针对对称带状矩阵B（一维数组存储）与通用矩阵A相加的性能瓶颈问题，结合Fortran95和gfortran-8编译器，提供以下优化方案：

1. 减少循环内重复计算与内存访问次数 #

循环内重复计算索引、多次读取同一数组元素是性能损耗的核心原因之一。提前计算固定值、用临时变量缓存数组元素，能有效降低开销：

subroutine iter(n,ccc,matrix_A,array_B,const)
    implicit none
    integer,intent(in)::n,ccc
    complex*16,intent(inout)::matrix_A(n,n)
    complex*16,intent(in)::array_B(n*ccc), const
    integer::i,j,start_idx,offset
    complex*16::b_val

    do i=1,n
        start_idx = (i-1)*ccc + 1  ! 提前计算当前行的起始索引，避免重复计算
        ! 处理对角线元素(i,i)
        matrix_A(i,i) = matrix_A(i,i) + array_B(start_idx)*const
        ! 处理上三角及对称下三角元素
        do j=i+1,min(i+ccc-1,n)
            offset = j - i
            b_val = array_B(start_idx + offset)*const  ! 一次读取，多次使用
            matrix_A(i,j) = matrix_A(i,j) + b_val
            matrix_A(j,i) = matrix_A(j,i) + b_val
        end do
    end do
end subroutine iter

2. 开启编译器最高级优化选项 #

gfortran-8的优化选项能大幅提升代码性能，尤其针对循环和SIMD指令的利用，推荐编译时添加：

gfortran-8 -O3 -march=native -ffast-math -funroll-loops your_code.f90

• -O3：启用全级别优化，包括循环展开、函数内联等
• -march=native：针对当前CPU架构生成最优指令，充分利用AVX/SSE等SIMD扩展
• -ffast-math：允许数学运算的近似简化，对复数运算性能提升显著（需确保精度符合需求）
• -funroll-loops：强制展开内层小循环（本场景内层循环最多ccc-1次，非常适合）

3. 重构循环结构消除分支判断 #

原代码中的min(i+ccc-1,n)会引入条件分支，阻碍编译器向量化。将循环拆分为无边界判断的中间段和边界段，减少分支干扰：

subroutine iter(n,ccc,matrix_A,array_B,const)
    implicit none
    integer,intent(in)::n,ccc
    complex*16,intent(inout)::matrix_A(n,n)
    complex*16,intent(in)::array_B(n*ccc), const
    integer::i,j,mid_start,mid_end,start_idx,offset
    complex*16::b_val

    mid_start = ccc
    mid_end = n - ccc + 1

    ! 单独处理对角线元素
    do i=1,n
        matrix_A(i,i) = matrix_A(i,i) + array_B((i-1)*ccc +1)*const
    end do

    ! 处理前ccc-1行（存在上边界）
    do i=1,ccc-1
        do j=i+1,min(i+ccc-1,n)
            offset = j - i
            b_val = array_B((i-1)*ccc +1 + offset)*const
            matrix_A(i,j) = matrix_A(i,j) + b_val
            matrix_A(j,i) = matrix_A(j,i) + b_val
        end do
    end do

    ! 处理中间行（无边界判断，适合向量化）
    do i=mid_start,mid_end
        start_idx = (i-1)*ccc +1
        do j=i+1,i+ccc-1
            offset = j - i
            b_val = array_B(start_idx + offset)*const
            matrix_A(i,j) = matrix_A(i,j) + b_val
            matrix_A(j,i) = matrix_A(j,i) + b_val
        end do
    end do

    ! 处理后ccc-1行（存在下边界）
    do i=mid_end+1,n
        do j=i+1,n
            offset = j - i
            b_val = array_B((i-1)*ccc +1 + offset)*const
            matrix_A(i,j) = matrix_A(i,j) + b_val
            matrix_A(j,i) = matrix_A(j,i) + b_val
        end do
    end do
end subroutine iter

4. 按带状偏移量重构循环 #

将循环按带状的偏移量k（从1到ccc-1）拆分，让内存访问模式更连续，便于编译器优化：

subroutine iter(n,ccc,matrix_A,array_B,const)
    implicit none
    integer,intent(in)::n,ccc
    complex*16,intent(inout)::matrix_A(n,n)
    complex*16,intent(in)::array_B(n*ccc), const
    integer::i,k,idx_b
    complex*16::b_val

    ! 处理对角线（k=0）
    do i=1,n
        matrix_A(i,i) = matrix_A(i,i) + array_B((i-1)*ccc +1)*const
    end do

    ! 处理每个偏移量k对应的带状元素
    do k=1,ccc-1
        do i=1,n-k
            idx_b = (i-1)*ccc +1 + k
            b_val = array_B(idx_b)*const
            matrix_A(i,i+k) = matrix_A(i,i+k) + b_val
            matrix_A(i+k,i) = matrix_A(i+k,i) + b_val
        end do
    end do
end subroutine iter

这种方式下，matrix_A(i+k,i)访问的是列i的连续元素（Fortran列优先存储），内存访问更高效，编译器更容易实现向量化。

5. 显式启用SIMD向量化 #

如果编译器自动向量化效果不佳，可通过!$OMP SIMD指令强制对循环进行向量化优化（需确保编译时开启OpenMP支持：-fopenmp）：

! 处理中间行时添加SIMD指令
do i=mid_start,mid_end
    start_idx = (i-1)*ccc +1
    !$OMP SIMD
    do j=i+1,i+ccc-1
        offset = j - i
        b_val = array_B(start_idx + offset)*const
        matrix_A(i,j) = matrix_A(i,j) + b_val
        matrix_A(j,i) = matrix_A(j,i) + b_val
    end do
end do

内容的提问来源于stack exchange，提问作者ZoiloN