优化七个循环

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

这些事情的诀窍是寻找常见的模式,并使用现有的有效例程来加速它们.

像往常一样,MSB是完全正确的,只要翻转你的指数就会给你带来可观的加速,尽管具有高优化性能的英特尔fortran编译器已经为你带来了一些好处.

但是让我们剥离m索引一秒钟(这很容易做,正如MSB指出的那样,这是最慢的移动指数)并且只看一下乘法:

A _i,j,k,l =ΣBk _,l,r,s ×P _i,j,r,s
A _i,j,k,l =ΣPi _,j,r,s ×B _{k,l ,R,S}

重塑数组:

A _ij,kl = _ΣPij,rs ×B _kl,rs
A _ij,kl = _ΣPij,rs ×B ^T _rs,kl
A = P×B ^T

我们现在有矩阵乘法,其中存在非常有效的例程.因此,如果我们重塑P和B矩阵并转置B,我们可以进行简单的矩阵乘法并重塑结果; 在这种情况下,这种重塑甚至不一定需要任何副本.所以改变这样的事情:

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调用中进行重新整形时它变得更快,如果你-mkl=parallel用于线程例程,它会更快(对于大矩阵).如果您还没有MKL,您可以链接到其他一些快速LAPACK _GEMM例程.