Fortran cumsum累计和分块优化性能提升原因咨询

云原生阿岚

2026-10-5

cumsum分段优化性能提升原因咨询

cumsum函数接收向量输入，输出对应元素的累计求和向量。我需要了解自己尝试的cumsum优化方案的性能提升逻辑。

我采用的优化思路如下：不再对整个向量执行单一遍历循环，而是编写循环同时对向量的4个分段分别做局部累加，之后再根据前序分段的总和调整后序分段的数值，该方案带来的微小误差在可接受范围内。

完整测试程序

module cumsum_mod
    implicit none
    integer, parameter, private :: dp = kind(1d0)
contains
    ! cumsum in one straight loop
    subroutine cumsum1(n, a, b)
        integer :: n, i
        real(dp) :: a(n), b(n)
        
        b(1) = a(1)
        do i = 2, n
            b(i) = a(i) + b(i-1)
        end do
    end subroutine
    
    subroutine cumsum2(n, a, b)
        integer :: n, i, m
        real(dp) :: a(n), b(n)
        
        m = n/4
        
        ! Loop over the four parts
        b(1) = a(1)
        b(1+m) = a(1+m)
        b(1+2*m) = a(1+2*m)
        b(1+3*m) = a(1+3*m)
        do i = 2, m
            b(i) = a(i) + b(i-1)
            b(i+m) = a(i+m) + b(i+m-1)
            b(i+2*m) = a(i+2*m) + b(i+2*m-1)
            b(i+3*m) = a(i+3*m) + b(i+3*m-1)
        end do
        
        ! Adjusting
        b(m+1:2*m) = b(m+1:2*m) + b(m)
        b(2*m+1:3*m) = b(2*m+1:3*m) + b(2*m)
        b(3*m+1:4*m) = b(3*m+1:4*m) + b(3*m)
        
        do i = 4*m+1, n
            b(i) = a(i) + b(i-1)
        end do
    end subroutine
    
    subroutine cumsum3(n, a, b)
        integer :: n, i, m
        real(dp) :: a(n), b(n)
        real(dp) :: k1, k2, k3
        
        m = n/4
        
        ! Loop over the four parts
        b(1) = a(1)
        b(1+m) = a(1+m)
        b(1+2*m) = a(1+2*m)
        b(1+3*m) = a(1+3*m)
        do i = 2, m
            b(i) = a(i) + b(i-1)
            b(i+m) = a(i+m) + b(i+m-1)
            b(i+2*m) = a(i+2*m) + b(i+2*m-1)
            b(i+3*m) = a(i+3*m) + b(i+3*m-1)
        end do
        
        ! Adjusting
        k1 = b(m)
        k2 = b(2*m) + k1
        k3 = b(3*m) + k2
        do i = 1, m
            b(i+m) = b(i+m) + k1
            b(i+2*m) = b(i+2*m) + k2
            b(i+3*m) = b(i+3*m) + k3
        end do
        
        do i = 4*m+1, n
            b(i) = a(i) + b(i-1)
        end do
    end subroutine
end module

program cumsum_test
    use cumsum_mod
    implicit none
    integer, parameter :: dp = kind(1d0)
    real(dp), allocatable :: a(:), b(:), b1(:), b2(:), b3(:)
    integer :: n, m, i
    real(dp) :: t1, t2
    
    read *, n, m
    
    allocate (a(n), b(n), b1(n), b2(n), b3(n))
    call random_number(a)
    
    ! Heating up
    do i = 1, 20
        call cumsum1(n, a, b)
        call cumsum2(n, a, b)
        call cumsum3(n, a, b)
    end do
    
    call cpu_time(t1)
    do i = 1, m
        call cumsum1(n, a, b)
    end do
    call cpu_time(t2)
    print *, t2-t1
    
    call cpu_time(t1)
    do i = 1, m
        call cumsum2(n, a, b)
    end do
    call cpu_time(t2)
    print *, t2-t1

    call cpu_time(t1)
    do i = 1, m
        call cumsum3(n, a, b)
    end do
    call cpu_time(t2)
    print *, t2-t1
    
    ! Quick check of the difference
    call cumsum1(n, a, b1)
    call cumsum2(n, a, b2)
    call cumsum3(n, a, b3)
    print *, maxval(abs(b2-b1))
    print *, maxval(abs(b3-b1))
    deallocate (a, b, b1, b2, b3)
end program

优化方案分为两个版本：调整阶段要么使用4个数组切片语句实现，例如b(m+1:2*m) = b(m+1:2*m) + b(m)，要么使用一个循环同步处理三个后序分段。

测试环境

编译器：Intel Fortran，编译参数为/O3 /QxHost
运行硬件：支持AVX2指令集的Intel Core i7 9700F

汇编核心片段对比

我不贴出全部250KB的汇编输出，仅展示第一个累加循环的核心部分：

cumsum1 核心汇编

L1:
    inc       rcx
    vaddsd    xmm0, xmm0, QWORD PTR [8+rdx+r10]
    vmovsd    QWORD PTR [8+rdx+r8], xmm0
    vaddsd    xmm0, xmm0, QWORD PTR [16+rdx+r10]
    vmovsd    QWORD PTR [16+rdx+r8], xmm0
    add       rdx, 16
    cmp       rcx, r11
    jb        L1

cumsum2 与 cumsum3 核心汇编

L1:
    vaddsd    xmm0, xmm0, QWORD PTR [8+rdi+r14*8]
    vaddsd    xmm1, xmm1, QWORD PTR [8+r12+r14*8]
    vaddsd    xmm2, xmm2, QWORD PTR [8+r11+r14*8]
    vaddsd    xmm3, xmm3, QWORD PTR [8+r9+r14*8]
    vmovsd    QWORD PTR [8+r8+r14*8], xmm0
    vmovsd    QWORD PTR [8+rcx+r14*8], xmm1
    vmovsd    QWORD PTR [8+rbp+r14*8], xmm2
    vmovsd    QWORD PTR [8+rsi+r14*8], xmm3
    inc       r14
    cmp       r14, r10
    jb        L1

因此两个版本的累加循环本质上都是「标量」循环，未用到并行SIMD操作，仅调整阶段如预期用到了并行指令。

问题核心

为什么cumsum2/cumsum3的运行速度远快于cumsum1？实测性能提升数据如下：长度1000时快3倍，长度10000时快2倍，长度100000时仍快30%~60%：

length=1000 loops=100000000
cumsum1   85.8593750000000
cumsum2   27.2187500000000
cumsum3   28.0937500000000

length=10000 loops=1000000
cumsum1   8.78125000000000
cumsum2   4.51562500000000
cumsum3   4.56250000000000
  
length=100000 loops=1000000
cumsum1   87.8281250000000
cumsum2   52.9687500000000
cumsum3   63.3281250000000

我最初猜测是加法并行执行带来的提升，但所有版本的累加阶段都是串行执行，且优化版还有额外的数值调整步骤，性能仍远高于原生版本。我怀疑该现象和内存缓存有关，但未找到明确依据。

恳请各位提供相关解释，也欢迎给出进一步优化的可行方案。

补充说明（回应评论疑问）

cumsum3中第二个调整循环的汇编如下：

L1:
    vaddpd    ymm6, ymm5, YMMWORD PTR [rcx+r10*8]
    vmovupd   YMMWORD PTR [rcx+r10*8], ymm6
    vaddpd    ymm6, ymm4, YMMWORD PTR [rbp+r10*8]
    vmovupd   YMMWORD PTR [rbp+r10*8], ymm6
    vaddpd    ymm6, ymm3, YMMWORD PTR [rsi+r10*8]
    vmovupd   YMMWORD PTR [rsi+r10*8], ymm6
    vaddpd    ymm6, ymm5, YMMWORD PTR [32+rcx+r10*8]
    vmovupd   YMMWORD PTR [32+rcx+r10*8], ymm6
    vaddpd    ymm6, ymm4, YMMWORD PTR [32+rbp+r10*8]
    vmovupd   YMMWORD PTR [32+rbp+r10*8], ymm6
    vaddpd    ymm6, ymm3, YMMWORD PTR [32+rsi+r10*8]
    vmovupd   YMMWORD PTR [32+rsi+r10*8], ymm6
    vaddpd    ymm6, ymm5, YMMWORD PTR [64+rcx+r10*8]
    vmovupd   YMMWORD PTR [64+rcx+r10*8], ymm6
    vaddpd    ymm6, ymm4, YMMWORD PTR [64+rbp+r10*8]
    vmovupd   YMMWORD PTR [64+rbp+r10*8], ymm6
    vaddpd    ymm6, ymm3, YMMWORD PTR [64+rsi+r10*8]
    vmovupd   YMMWORD PTR [64+rsi+r10*8], ymm6
    vaddpd    ymm6, ymm5, YMMWORD PTR [96+rcx+r10*8]
    vmovupd   YMMWORD PTR [96+rcx+r10*8], ymm6
    vaddpd    ymm6, ymm4, YMMWORD PTR [96+rbp+r10*8]
    vmovupd   YMMWORD PTR [96+rbp+r10*8], ymm6
    vaddpd    ymm6, ymm3, YMMWORD PTR [96+rsi+r10*8]
    vmovupd   YMMWORD PTR [96+rsi+r10*8], ymm6
    add       r10, 16
    cmp       r10, r9
    jb        L1

而cumsum2的相同逻辑通过三个循环实现，每个循环结构如下：

L1:
    vaddpd    ymm2, ymm1, YMMWORD PTR [rcx+r11*8]
    vaddpd    ymm3, ymm1, YMMWORD PTR [32+rcx+r11*8]
    vaddpd    ymm4, ymm1, YMMWORD PTR [64+rcx+r11*8]
    vaddpd    ymm5, ymm1, YMMWORD PTR [96+rcx+r11*8]
    vmovupd   YMMWORD PTR [rcx+r11*8], ymm2
    vmovupd   YMMWORD PTR [32+rcx+r11*8], ymm3
    vmovupd   YMMWORD PTR [64+rcx+r11*8], ymm4
    vmovupd   YMMWORD PTR [96+rcx+r11*8], ymm5
    add       r11, 16
    cmp       r11, r12
    jb        L1

回答

性能差异核心原因

1. 指令级并行（ILP）利用率差异

这是最核心的原因：
cumsum1的累加循环只有1条依赖链：b(i) = a(i) + b(i-1)，每一步计算都必须等上一步的结果输出才能执行，现代CPU的乱序执行能力完全发挥不出来，执行端口长期闲置，流水线里充满等待气泡。
而cumsum2/3的累加循环有4条完全独立的依赖链（四个分段的局部累加互不影响），i7-9700F处理器有3个浮点加法执行端口，乱序执行逻辑可以同时调度这4条独立的加法链，把执行端口的利用率拉满，相当于在标量层面实现了并行计算。这也是为什么小数据量（完全在L1缓存内，无内存瓶颈）下性能提升接近3~4倍，刚好匹配CPU的并行执行能力。

2. 调整阶段的SIMD加速抵消了额外开销

调整操作是给整个分段的所有元素加同一个固定值，没有数据依赖，编译器生成了256位的vaddpd指令，一次可以处理4个双精度数，这部分的执行效率非常高，额外的调整开销远低于累加阶段ILP提升带来的收益。

3. 性能随数据长度衰减的原因

当向量长度达到10万时，数据已经超出L1/L2缓存容量，瓶颈从计算侧转移到内存带宽侧，此时ILP的收益被内存访问速度限制，所以性能提升幅度降到30%~60%。
另外cumsum2比cumsum3性能更好，是因为cumsum2的调整阶段是连续访问单个数组切片，缓存命中率和访存预取效率远高于cumsum3循环里交替访问三个分段的跳址访问模式。

进一步优化方案

调整分段数：如果后续升级到支持AVX512的处理器，可以把分段数改成8，进一步提升ILP利用率，同时调整阶段可以用512位SIMD指令获得更高吞吐量。
内存对齐优化：给数组分配时添加对齐属性，让编译器可以生成对齐的vmovapd指令代替非对齐的vmovupd，访存效率可以提升10%~15%。
多级分块优化：对于超大规模向量，可以用二级分块策略：先做16个分段的局部累加，再做两级全局调整，进一步压榨CPU的乱序执行能力。
边界处理优化：当前代码里最后剩余的4*m+1 ~ n的元素还是用单链累加，可以把这部分也合并到分段逻辑里，消除剩余的串行瓶颈。

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

Fortran cumsum累计和分块优化性能提升原因咨询

完整测试程序 #

测试环境 #

汇编核心片段对比 #

cumsum1 核心汇编 #

cumsum2 与 cumsum3 核心汇编 #

问题核心 #

补充说明（回应评论疑问） #

性能差异核心原因 #

1. 指令级并行（ILP）利用率差异 #

2. 调整阶段的SIMD加速抵消了额外开销 #

3. 性能随数据长度衰减的原因 #

进一步优化方案 #