在 fortran 中优化数组乘法答案

【问题标题】：Optimizing array multiplication in fortran在 fortran 中优化数组乘法
【发布时间】：2022-01-20 13:52:37
【问题描述】：

我正在尝试在fortran中优化以下类型的数组乘法：

    do i = 1, n
    do j = 1, n
    do k = 1, n
            do i1 = 1, n
            do j1 = 1, n
            do k1 = 1, n
                    B(k1,j1,i1) = B(k1,j1,i1) + &
                    A(k,j,i)*v1(k,k1)*v2(j,j1)*v3(i,i1)
            enddo
            enddo
            enddo
    enddo
    enddo
    enddo

这里 A 和 B 的大小为 (n,n,n)，而 v1、v2 和 v3 的大小为 (n,n)。似乎应该有某种矢量化可以让我加快这个评估。我尝试了不同的索引排序。有没有更快的计算 B 的方法？

编辑：

好的，可以使用 matmul 获得大约 10 倍的加速：

    do i = 1, n
    do i1 = 1, n
    do j = 1, n
    do j1 = 1, n
            B(:,j1,i1) = B(:,j1,i1) + &   
            matmul(A(:,j,i),v1)*v2(j,j1)*v3(i,i1)                                             
    enddo
    enddo
    enddo
    enddo

我们还能做得更好吗？

【问题讨论】：

只是为了让我可以在脑海中梳理一些想法，大概 n 有多大？
n 通常低于 10，但在极端情况下，我预计它甚至会达到数百。目前，我正在尝试 50 左右的 ns。
急着出门，但我怀疑你可以将订单从N^6降低到N^4。注意 Sum_i A(k,j,i)*v(i,i1) 独立于 k1 和 j1 ，因此可以吊出循环 - 我怀疑反过来这样做（但我还没有说服自己）将计算的顺序减少到 3 组四重嵌套循环，而不是 1 组六深。
啊，这就是我看到的@veryreverie 似乎正在走向

标签： loops multidimensional-array fortran

【解决方案1】：

您至少可以使用 matmul 或 BLAS/LAPACK 将两个循环对减少为矩阵乘法，例如

do i1=1,n
  do i=1,n
    B(:,:,i1) = B(:,:,i1) &
            & + matmul(matmul(transpose(v1), A(:,:,i)), v2) * v3(i,i1)
  enddo
enddo

几乎可以肯定，您可以通过缓存双精度 matmul 的结果来进一步加快速度，例如

real :: intermediate(n,n,n)

do i=1,n
  intermediate(:,:,i) = matmul(matmul(transpose(v1), A(:,:,i)), v2)
enddo

do i1=1,n
  do i=1,n
    B(:,:,i1) = B(:,:,i1) + intermediate(:,:,i) * v3(i,i1)
  enddo
enddo

这可能会被更多matmul 加速，例如

real :: intermediate(n,n,n)

do i=1,n
  intermediate(:,:,i) = matmul(matmul(transpose(v1), A(:,:,i)), v2)
enddo

do j1=1,n
  B(:,j1,:) = matmul(intermediate(:,j1,:), v3))
enddo

【讨论】：

非常感谢！这很好用。发布问题后，我立即意识到 matmul 可以帮助我减少其中一个循环。这个技巧可以将它再加速几次。