在 OpenMP 中使用嵌套循环计算协方差

Question

我有一段核心代码用于计算模型状态和观测值之间的协方差，该协方差将运行数十万次，我想使用 OpenMP 来加快速度。目前的实施似乎并没有像我希望的那样加快速度。有一个计算协方差的嵌套循环。 N很大（10^4左右），M比较小（10^2左右）。内部循环计算两个总和并将一个总和除以另一个。关于如何加快速度的任何建议？

  xa     = ...
  hxens  = ...
  istate = ...

  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(iupdt, iobs, mem_count, iens)                                            
  !$OMP DO
  do iupdt = 1,N                                        
     mem_count = 0

     do iens = 1,M
           cij(iupdt) = cij(iupdt) + xa(istate,iens) * hxens(iupdt,iens)   
           mem_count = mem_count + 1
     end do
     cij(iupdt) = cij(iupdt)/(mem_count-1)          
  end do
  !$OMP END DO          
  !$OMP END PARALLEL

非常感谢任何帮助！

Answer 1

1) 您可以删除 mem_count ：总和始终为 M.

2) 您可以展开外循环：如果您有 Haswell CPU，您将能够进行矢量化 FMA。在 32 字节边界上对齐数组将提高使用编译器选项或指令的性能。

3) 如果你可以转到单精度，你可以获得额外的两倍。可以单精度累加，双精度累加 divide（注意 xa * hxens 总是同号，会很不精确...）

4) 您只能计算一次 1.d0/(dble(M)-1.d0) 以避免执行昂贵的 div 操作。

如果你有英特尔编译器，你在哪里分配你的数组：

  !DIR$ ATTRIBUTES ALIGN : 32 :: cij, hxens

然后，

  n8 = (n/8)*8                             
  f = 1.d0/(dble(M)-1.d0)

  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(iupdt, iobs, iens, l)       
  !$OMP DO
  do iupdt = 1,n8,8                                        

     do iens = 1,M 
       do l=0,7
         cij(iupdt+l) = cij(iupdt+l) + xa(istate,iens) * hxens(iupdt+l,iens)   
       end do
     end do

     do l=0,7
       cij(iupdt+l) = cij(iupdt+l) * f
     end do

  end do
  !$OMP END DO

  !$OMP END PARALLEL

  do iupdt = n8+1,n                                        
     do iens = 1,M
        cij(iupdt) = cij(iupdt) + xa(istate,iens) * hxens(iupdt,iens)   
     end do
     cij(iupdt) = cij(iupdt)*f          
  end do

Answer 2

代码中的主要问题是数组 xa 和 hxens 的访问顺序。但是，如果交换 iupdt 和 iens 循环，代码会表现得更好。如果 cij 的最后一个除法从 iupdt 循环中移除，并设置在它自己的新循环中，这就很容易做到。另外，必须注意 mem_count 的累加是无用的，因为结果总是 M.

为了能够编译代码，我将代码片段放在了一个子例程中。这是最终的样子：

subroutine covariance( xa, hxens, cij, istate, N, M )
    implicit none
    integer, intent(in)             :: istate, N, M 
    double precision, intent(in)    :: xa(N,M), hxens(N,M)
    double precision, intent(inout) :: cij(N)

    integer          :: iupdt, iens
    double precision :: ratio, xaiens

    ratio = 1.d0 / (M-1)
    !$omp parallel private( xaiens, iens, iupdt )
    !$omp do
    do iens = 1,M
        xaiens = xa(istate,iens)
        !$omp simd
        do iupdt = 1,N
             cij(iupdt) = cij(iupdt) + xaiens * hxens(iupdt,iens)   
        end do
        !$omp end simd
    end do
    !$omp end do
    !$omp do simd
    do iupdt = 1,N
        cij(iupdt) = cij(iupdt) * ratio          
    end do
    !$omp end do simd
    !$omp end parallel
end subroutine covariance

几点说明：

所有内存访问现在都是线性的，这很好
然而，cij 以前可能只被访问过一次（因为在 iens 循环期间缓存），现在被访问 iens 次。这可能非常糟糕。但是，由于您说 N 大约为 10000，这意味着 cij 大约为 80KB，这将适合缓存。所以这不应该有任何后果。
最里面的循环现在可以简单地向量化。此外，它非常适合使用融合的乘法和加法运算，从而可以充分发挥处理器的潜力。我已经放置了一个 omp simd 指令来加强这一点，但如果您的编译器不支持它，您可以将其删除。无论如何，矢量化案例是如此明显，以至于我怀疑编译器是否需要它。最后一个循环也是如此。
计算如此有限且实际数据适合缓存，我不确定 OpenMP 并行化是否如此 effective/necessary。我鼓励您在使用和不使用它的情况下尝试代码，并删除全部或部分代码，如果它被证明是适得其反的。我已经把所有可以合理设置的东西都放在这里了，但不要犹豫，删掉需要的东西。

在 OpenMP 中使用嵌套循环计算协方差

Covariance calculation with a nested loop in OpenMP

fortran

openmp