performance - 七次循环的优化

Question

我有 3 个数组，我必须做这个求和

实现的代码是

do i=1,320
  do j=1,320
    do k=1,10
     do l=1,10
      do m=1,10
       do r=1,10
        do s=1,10
          sum=sum+B(k,l,r,s,m)*P(i,j,r,s,m)
        end do
       end do
       A(i,j,k,l,m)=sum
     end do 
    end do 
   end do 
 end do
end do

执行代码需要 1 天。 有什么办法可以优化吗？

谢谢。

Answer 1

这些事情的诀窍是寻找共同的模式并使用现有的高效例程来加速它们。

M.S.B 和往常一样完全正确，只需翻转索引即可显着加快速度，尽管具有高度优化的英特尔 Fortran 编译器已经为您带来了一些好处。

但让我们剥掉 m索引一秒钟(这很容易做到，正如 MSB 指出的那样，这是移动最慢的索引)，然后只看乘法:

A_i,j,k,l = ∑ B_k,l,r,s × P_i,j,r,s
A_i,j,k,l = ∑ P_i,j,r,s × B_k,l,r,s

reshape 数组:

A_ij,kl = ∑ P_ij,rs × B_kl,rs
A_ij,kl = ∑ P_ij,rs × B^T_rs,kl
A = P × B^T

我们现在有矩阵乘法，为此存在非常有效的例程。因此，如果我们对 P 和 B 矩阵进行整形，并转置 B，我们可以进行简单的矩阵乘法并对结果进行整形；在这种情况下，这种 reshape 甚至不一定需要任何副本。所以改变这样的东西:

program testpsum
implicit none

integer, dimension(10,10,10,10,10) :: B
integer, dimension(32,32,10,10,10) :: P
integer, dimension(32,32,10,10,10) :: A
integer :: psum
integer :: i, j, k, l, m, r, s

B = 1
P = 2

do i=1,32
  do j=1,32
    do k=1,10
     do l=1,10
      do m=1,10
       do r=1,10
        do s=1,10
          psum=psum+B(k,l,r,s,m)*P(i,j,r,s,m)
        end do
       end do
       A(i,j,k,l,m)=psum
       psum = 0
     end do
    end do
   end do
 end do
end do

print *,minval(A), maxval(A)

end program testpsum

对此:

program testmatmult
implicit none

integer, dimension(10,10,10,10,10) :: B
integer, dimension(32,32,10,10,10) :: P
integer, dimension(10*10,10*10) :: Bmt
integer, dimension(32*32,10*10) :: Pm
integer, dimension(32,32,10,10,10) :: A
integer :: m

B = 1
P = 2

do m=1,10
    Pm  = reshape(P(:,:,:,:,m),[32*32,10*10])
    Bmt = transpose(reshape(B(:,:,:,:,m),[10*10,10*10]))
    A(:,:,:,:,m) = reshape(matmul(Pm,Bmt),[32,32,10,10])
end do

print *,minval(A), maxval(A)

end program testmatmult

给出时间:

$ time ./psum
         200         200

real    0m2.239s
user    0m1.197s
sys 0m0.008s

$ time ./matmult
         200         200

real    0m0.064s
user    0m0.027s
sys 0m0.008s

使用 ifort -O3 -xhost -mkl 编译时所以我们可以使用快速的英特尔 MKL 库。 如果您不创建 Pm，它会变得更快临时的，只需在 matmult 调用中进行 reshape ，如果使用 -mkl=parallel，速度会更快(对于大矩阵)对于线程例程。如果您还没有 MKL，则可以链接到其他一些快速 LAPACK _GEMM 例程。