2

我正在尝试在 Fortran 中编写一个函数，它将多个具有不同权重的矩阵相乘，然后将它们加在一起形成一个矩阵。我已经确定这个过程是我程序中的瓶颈（这个权重将在程序的单次运行中进行多次，具有不同的权重）。现在我正试图通过从 Matlab 切换到 Fortran 来使其运行得更快。我是 Fortran 的新手，所以我感谢所有帮助。

在 Matlab 中，我发现进行此类计算的最快方法如下所示：

function B = weight_matrices()
n = 46;
m = 1800;
A = rand(n,m,m);
w = rand(n,1);
tic;
B = squeeze(sum(bsxfun(@times,w,A),1));
toc;

在我的机器（Matlab R2012b，MacBook Pro 13" Retina，2.5 GHz Intel Core i5，8 GB 1600 MHz DDR3）上分配的行在B大约 0.9 秒内运行。应该注意的是，对于我的问题，张量A将是整个程序运行（初始化后）相同（常量），但 w 可以取任何值。此外，这里使用n和的典型值m，这意味着张量A在内存中的大小约为 1 GB。

我能想到的用 Fortran 写这个的最清晰的方法是这样的：

pure function weight_matrices(w,A) result(B)
    implicit none
    integer, parameter :: n = 46
    integer, parameter :: m = 1800
    double precision, dimension(num_sizes), intent(in) :: w
    double precision, dimension(num_sizes,msize,msize), intent(in) :: A
    double precision, dimension(msize,msize) :: B
    integer :: i
    B = 0
    do i = 1,n
        B = B + w(i)*A(i,:,:)
    end do
end function weight_matrices

当使用 gfortran 4.7.2 编译时，使用 -O3（使用“call cpu_time(t)”计时的函数调用），此函数在大约 1.4 秒内运行。如果我手动将循环解包到

B = w(1)*A(1,:,:)+w(2)*A(2,:,:)+ ... + w(46)*A(46,:,:)

该函数需要大约 0.11 秒才能运行。这很棒，意味着与 Matlab 版本相比，我得到了大约 8 倍的加速。但是，我仍然对可读性和性能有一些疑问。

首先，我想知道是否有更快的方法来执行矩阵的加权和求和。我查看了 BLAS 和 LAPACK，但找不到任何似乎适合的功能。我还尝试将A枚举矩阵的维度作为最后一个维度（即从元素切换(i,j,k)到(k,i,j)元素），但这会导致代码变慢。

其次，这个快速版本不是很灵活，而且实际上看起来很丑，因为对于这样一个简单的计算来说，它包含了太多的文本。对于我正在运行的测试，我想尝试使用不同数量的权重，以便 w 的长度会有所不同，以了解它如何影响我的算法的其余部分。但是，这意味着我B每次重写作业都相当繁琐。有什么方法可以让这更灵活，同时保持性能相同（或更好）？

第三，A如前所述，张量将在程序运行期间保持不变。我在我的程序中使用他们自己模块中的“参数”属性设置了常量标量值，并使用“使用”表达式将它们导入到需要它们的函数/子例程中。为张量做等效事情的最佳方法是什么A？我想告诉编译器这个张量在 init. 之后将是常量，以便可以进行任何相应的优化。请注意，A它的大小通常约为 1 GB，因此直接在源文件中输入它是不切实际的。

提前感谢您的任何意见！:)

4

4 回答 4

3

也许你可以尝试类似

    do k=1,m
       do j=1,m
          B(j,k)=sum( [ ( (w(i)*A(i,j,k)), i=1,n) ])
       enddo
    enddo

方括号是 (/ /) 的一种更新形式，即一维矩阵（向量）。insum是一个维度矩阵，对(n)所有sum这些元素求和。这正是您展开的代码所做的（并且不完全等于do您拥有的循环）。

于 2013-04-20T21:05:35.793 回答
1

我不会隐藏任何循环，因为这通常比较慢。您可以显式地编写它，然后您会看到内部循环访问超过了最后一个索引，从而使其效率低下。因此，您应该通过存储 A 来确保您的n维度是最后一个A(m,m,n)：

B = 0
do i = 1,n
    w_tmp = w(i)
    do j = 1,m
        do k = 1,m
            B(k,j) = B(k,j) + w_tmp*A(k,j,i)
        end do
    end do
end do

这应该更有效，因为您现在正在内部循环中访问内存中的连续元素。

另一种解决方案是使用 1 级 BLAS 子例程 _AXPY (y = a*x + y)：

B = 0
do i = 1,n
    CALL DAXPY(m*m, w(i), A(1,1,i), 1, B(1,1), 1)
end do

使用英特尔 MKL，这应该更有效，但您应该再次确保最后一个索引是在外部循环中发生变化的索引（在这种情况下是您正在编写的循环）。您可以在此处找到此调用的必要参数：MKL

编辑：您可能还想使用一些并行化？（我不知道 Matlab 是否利用了这一点）

EDIT2：在凯尔的回答中，内部循环超过不同的值w，这比n重新加载的时间更有效，B因为w可以保存在缓存中（使用A(n,m,m)）：

B = 0
do i = 1,m
    do j = 1,m
        B(j,i)=0.0d0
        do k = 1,n
            B(j,i) = B(j,i) + w(k)*A(k,j,i)
        end do
    end do
end do

这种显式循环的性能比使用全数组操作的 Kyle 代码要好 10%。带宽ifort -O3 -xHost约为 6600 MB/s，gfortran -O3约为 6000 MB/s，使用任一编译器的全阵列版本也约为 6000 MB/s。

于 2013-04-21T08:48:06.817 回答
1

我试图完善 Kyle Vanos 的解决方案。

因此，我决定使用sumFortran 的矢量功能。

我不知道结果是否正确，因为我只寻找时间！

版本 1：（用于比较）

B = 0
do i = 1,n
    B = B + w(i)*A(i,:,:)
end do

版本 2：（来自 Kyle Vanos）

do k=1,m
   do j=1,m
      B(j,k)=sum( [ ( (w(i)*A(i,j,k)), i=1,n) ])
   enddo
enddo

版本 3：（混合索引，一次处理一行/列）

do j = 1, m
    B(:,j)=sum( [ ( (w(i)*A(:,i,j)), i=1,n) ], dim=1)
enddo

版本 4：（完整矩阵）

B=sum( [ ( (w(i)*A(:,:,i)), i=1,n) ], dim=1)

定时

如您所见，我不得不混合索引以获得更快的执行时间。第三种解决方案真的很奇怪，因为矩阵的数量是中间索引，但这是内存顺序原因所必需的。

V1: 1.30s
V2: 0.16s
V3: 0.02s
V4: 0.03s

最后，我想说的是，如果您有可能以任意顺序更改矩阵索引的顺序，您可以获得巨大的加速。

于 2013-04-22T09:40:18.557 回答
0

我知道这是一篇旧帖子，但是当我使用大多数已发布的解决方案时，我很乐意带来我的贡献。

通过为权重循环添加本地展开（来自​​ Steabert 的回答），与完整展开版本相比（使用不同大小的矩阵从 10% 到 80%），我可以稍微加快速度。部分展开可以帮助编译器在一个 SSE 调用中向量化 4 个操作。

pure function weight_matrices_partial_unroll_4(w,A) result(B)
  implicit none
  integer, parameter  :: n = 46
  integer, parameter  :: m = 1800
  real(8), intent(in) :: w(n)
  real(8), intent(in) :: A(n,m,m)
  real(8)             :: B(m,m)
  real(8)             :: Btemp(4)
  integer             :: i, j, k, l, ndiv, nmod, roll
  !==================================================
  roll = 4
  ndiv = n / roll
  nmod = mod( n, roll )

  do i = 1,m
    do j = 1,m
        B(j,i)=0.0d0
        k = 1
        do l = 1,ndiv
          Btemp(1) = w(k  )*A(k  ,j,i)
          Btemp(2) = w(k+1)*A(k+1,j,i)
          Btemp(3) = w(k+2)*A(k+2,j,i)
          Btemp(4) = w(k+3)*A(k+3,j,i)
          k = k + roll
          B(j,i) = B(j,i) + sum( Btemp )
        end do

        do l = 1,nmod !---- process the rest of the loop
          B(j,i) = B(j,i) + w(k)*A(k,j,i)
          k = k + 1
        enddo
    end do
  end do
end function
于 2018-05-30T10:06:10.170 回答