performance - 优化循环 : huge arrays operations

Question

我正在对不适合缓存的数组进行大量计算(这里是导数，但看起来类似于图像操作)，这意味着 CPU 必须在缓存中加载部分，计算，然后加载另一部分，等等。但是因为在计算的形状中，一些数据被加载、卸载和重新加载多次。我想知道是否有办法优化这一点。我已经在使用编译器优化(GCC 和 Intel)的 SIMD 指令。

这是 Fortran 计算，但它类似于 C/C++，内存顺序只是倒置，数组使用 () 而不是 []。 for 替换为 do。

在 x 轴上:

   do k=1,N(3)
      do j=1,N(2)
         do i=3,N(1)
            DF(i,j,k)=(F(i+1,j,k)-F(i-1,j,k))*B+(F(i-2,j,k)-F(i+2,j,k))*C
         end do
      end do
   end do

在 y 轴上:

   do k=1,N(3)
      do j=1,N(2)
         do i=3,N(1)
            DF(i,j,k)=(F(i,j+1,k)-F(i,j-1,k))*B+(F(i,j-2,k)-F(i,j+2,k))*C
         end do
      end do
   end do

在 z 轴上:

   do k=1,N(3)
      do j=1,N(2)
         do i=3,N(1)
            DF(i,j,k)=(F(i,j,k+1)-F(i,j,k-1))*B+(F(i,j,k-2)-F(i,j,k+2))*C
         end do
      end do
   end do

ax x 的一阶导数是可以的，因为内存是连续读取的。 y 轴和 z 轴上的导数不连续。

我结合了所有轴的最差计算(这是一个拉普拉斯算子):

 do k=1,N(3)
    do j=1,N(2)
       do i=3,N(1)
          V(i,j,k) = M(i,j,k,1) * p(i,j,k) &
                 & + M(i,j,k,2) * p(i-1,j,k) &
                 & + M(i,j,k,3) * p(i+1,j,k) &
                 & + M(i,j,k,4) * p(i,j-1,k) &
                 & + M(i,j,k,5) * p(i,j+1,k) &
                 & + M(i,j,k,6) * p(i,j,k-1) &
                 & + M(i,j,k,7) * p(i,j,k+1)
       end do
    end do
 end do

请注意，编译器不理解最后一个操作(拉普拉斯算子)。要使用 SIMD(矢量化计算)，我需要像这样拆分操作，这提供了 2.5 倍的加速:

 do k=1,N(3)
    do j=1,N(2)
       do i=3,N(1)
          V(i,j,k) = M(i,j,k,1) * p(i,j,k) &
                 & + M(i,j,k,2) * p(i-1,j,k) &
                 & + M(i,j,k,3) * p(i+1,j,k)
       end do
    end do
 end do
 do k=1,N(3)
    do j=1,N(2)
       do i=3,N(1)
          V(i,j,k) = V(i,j,k) + &
                 & + M(i,j,k,4) * p(i,j-1,k) &
                 & + M(i,j,k,5) * p(i,j+1,k)
       end do
    end do
 end do
 do k=1,N(3)
    do j=1,N(2)
       do i=3,N(1)
          V(i,j,k) = V(i,j,k) + &
                 & + M(i,j,k,6) * p(i,j,k-1) &
                 & + M(i,j,k,7) * p(i,j,k+1)
       end do
    end do
 end do

也许使用 SIMD 我已经达到了最大速度，但是因为这些计算需要几天时间，即使使用 MPI 和超过 1024 个 CPU，减少计算时间，即使是 20% 也是一个很大的进步!你们中有人对如何优化它有想法吗？

Answer 1

当您使用 3D 模板并引用 i,j,k-1 等元素时，i,j,k+1 是您通过的线性顺序该阵列将不是最佳的。缓存效率可以提高 loop tiling .

在我的代码中我使用

!$omp parallel private(i,j,k,bi,bj,bk)
!$omp do schedule(runtime) collapse(3)
do bk = 1, Unz, tnz
 do bj = 1, Uny, tny
  do bi = 1, Unx, tnx
   do k = bk, min(bk+tnz-1,Unz)
    do j = bj, min(bj+tny-1,Uny)
     do i = bi, min(bi+tnx-1,Unx)
         U2 (i,j,k) = U2(i,j,k) + &
           (U(i+1,j,k)-U(i,j,k)) * ...
         U2(i,j,k) = U2(i,j,k) - &
            (U(i,j,k)-U(i-1,j,k)) * ...
         U2(i,j,k) = U2(i,j,k) + &
            (U(i,j+1,k)-U(i,j,k)) * ...
         U2(i,j,k) = U2(i,j,k) - &
            (U(i,j,k)-U(i,j-1,k)) * ...
         U2(i,j,k) = U2(i,j,k) + &
            (U(i,j,k+1)-U(i,j,k)) * ...
         U2(i,j,k) = U2(i,j,k) - &
            (U(i,j,k)-U(i,j,k-1)) * ...
     end do
    end do
   end do
  end do
 end do
end do
!$omp end do

其中 tnx、tny、tnz 是您在 i,j,k 中迭代的图 block 的大小 顺序。大小必须设置为接近 L1 缓存。这将增加加载到缓存中的内容的重用。

如果您需要分开方向，您当然可以这样做并仍然保持平铺。