Fortran中OpenMP循环与数组操作的效率对比及技术疑问

我正在对比Fortran中OpenMP在不同循环索引组织方式及是否手动拆分索引情况下的执行效率，测试使用ifort编译器，编译选项为-qopenmp、-O3。

测试代码 #

Program test
    Use, Intrinsic :: iso_fortran_env, Only :  wp => real64, li => int64
    use omp_lib

    integer, parameter :: dp = selected_real_kind(15, 307)
    
    Real( dp ), Dimension( :, : ), Allocatable :: a
    Real( dp ), Dimension( :, : ), Allocatable :: c

    Integer :: na
    Integer :: i, j, m
    Integer( li ) :: start, finish, rate
    Integer :: numthreads
    real(dp) :: sum_time1
    
    Write( *, * ) 'numthreads'
    Read( *, * ) numthreads

    call omp_set_num_threads(numthreads)  

    Write( *, * ) 'na'
    Read( *, * ) na
    Allocate( a ( 1:na, 1:na ) ) 
    Allocate( c ( 1:na, 1:na ) ) 

    Call Random_number( a )


    sum_time1 = 0.0  
    Call System_clock( start, rate )
    do m = 1, m_iter    
      do i = 1, na
          c(i,:) = a(i,:)
      end do  
    end do  
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate  
    Write( *, * ) 'Time for loop first index', sum_time1 

    sum_time1 = 0.0  
    Call System_clock( start, rate )
      do i = 1, na
          c(:,i) = a(:,i)
      end do  
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate  
    Write( *, * ) 'Time for loop last index', sum_time1   


    sum_time1 = 0.0  
    Call System_clock( start, rate )
      do i = 1, na
        do j = 1, na
          c(i,j) = a(i,j)
        end do  
      end do  
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate  
    Write( *, * ) 'Time for loop two index first one', sum_time1


    sum_time1 = 0.0  
    Call System_clock( start, rate )
      do i = 1, na
        do j = 1, na
          c(j,i) = a(j,i)
        end do  
      end do  
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate  
    Write( *, * ) 'Time for loop two index inner most', sum_time1    

    sum_time1 = 0.0  
    Call System_clock( start, rate )
    !$omp parallel 
    !$omp do private(i)  
      do i = 1, na
        c(i,:) = a(i,:)
      end do  
    !$omp end do
    !$omp end parallel
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate 

    Write( *, * ) 'Time for loop-omp first index', sum_time1  


    sum_time1 = 0.0  
    Call System_clock( start, rate )
    !$omp parallel 
    !$omp do private(i)  
      do i = 1, na
        c(:,i) = a(:,i)
      end do  
    !$omp end do
    !$omp end parallel
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate 

    Write( *, * ) 'Time for loop-omp last index', sum_time1  


    sum_time1 = 0.0  
    Call System_clock( start, rate )
    !$omp parallel private(nthr, nstart, nend)
      nthr   = omp_get_thread_num()
      nstart = nthr * na / numthreads
      nend   = (nthr + 1) *  na / numthreads
      c(:, nstart+1:nend) = a(:, nstart+1:nend)
    !$omp end parallel
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate 

    Write( *, * ) 'Time for omp split last index', sum_time1   

    sum_time1 = 0.0  
    Call System_clock( start, rate )
    !$omp parallel private(nthr, nstart, nend)
      nthr   = omp_get_thread_num()
      nstart = nthr * na / numthreads
      nend   = (nthr + 1) *  na / numthreads
      c(nstart+1:nend, :) = a(nstart+1:nend, :)
    !$omp end parallel
    Call System_clock( finish, rate )
    sum_time1 = sum_time1 + Real( finish - start, dp ) / rate 

    Write( *, * ) 'Time for omp split first index', sum_time1   
  End

测试结果 #

numthreads
8
 na
5000
 Time for loop first index  0.116499780000000
 Time for loop last index  3.983250000000000E-002
 Time for loop two index first one  8.187200000000000E-003
 Time for loop two index inner most  8.229439999999999E-003
 Time for loop-omp first index  3.069090000000000E-003
 Time for loop-omp last index  3.028410000000000E-003
 Time for omp split last index  1.025315000000000E-002
 Time for omp split first index  9.189559999999999E-003

疑问解答 #

(i) 为何双层索引循环比(:,i)数组切片更快？ #

核心原因和Fortran的列优先存储特性、编译器优化能力有关：

• Fortran数组按列连续存储，双层循环中内层遍历列索引j时，是连续访问内存，缓存命中率极高；而(:,i)切片操作看似整列赋值，但单循环结构下，编译器可能无法完全消除数组描述符处理、边界检查等额外开销，也难最大化利用连续内存的缓存优势。
• 双层循环的逻辑更直白，编译器能轻松识别连续内存访问模式，生成向量化、预取等高效机器码；切片操作的抽象层级更高，编译器优化空间反而受限。

(ii) 为何手动拆分索引的OpenMP实现比!$omp do更慢？ #

主要有三个关键点：

• 负载不均：手动拆分nstart = nthr * na / numthreads的方式，当na无法被线程数整除时，最后一个线程会处理更多元素，导致部分线程闲置等待；而!$omp do会自动做负载均衡，把迭代更均匀地分配给线程。
• 优化受限：手动拆分后的数组切片操作，编译器难以对跨线程的切片做进一步向量化、内存预取优化；而!$omp do配合循环结构，编译器能结合OpenMP运行时库识别并行模式，生成更高效的代码。
• 额外开销：手动计算线程号、起止索引的操作会带来微小但累积的开销，而!$omp do的调度由成熟的OpenMP运行时库处理，开销更低。

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