performance - 不同优化下的 Fortran 矩阵乘法性能

Question

我正在阅读“使用 Fortran 进行科学软件开发”一书，其中有一个我认为非常有趣的练习:

“创建一个名为 MatrixMultiplyModule 的 Fortran 模块。向其中添加三个名为 LoopMatrixMultiply、IntrinsicMatrixMultiply 和 MixMatrixMultiply 的子例程。每个例程应将两个实矩阵作为参数，执行矩阵乘法，并通过第三个参数返回结果。LoopMatrixMultiply 应完全编写有do循环，没有数组操作或内在过程； IntrinsicMatrixMultiply应该使用matmul内在函数编写； MixMatrixMultiply应该使用一些do循环和内在函数dot_product编写。编写一个小程序来测试这三种不同方式的性能对不同大小的矩阵执行矩阵乘法。”

我对两个 2 阶矩阵的乘法进行了一些测试，以下是不同优化标志下的结果:

compiler:ifort version 13.0.0 on Mac 

这是我的问题:

为什么在 -O0 下它们具有大致相同的性能，但使用 -O3 时 matmul 具有巨大的性能提升，而显式循环和点积的性能提升较小？另外，为什么 dot_product 似乎与显式循环相比具有相同的性能？

我使用的代码如下:

module MatrixMultiplyModule

contains
    subroutine LoopMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))
        mtx3=0.

        do i=1,m
            do j=1,n
                do k=1,size(mtx1,dim=2)
                    mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
                end do      
            end do
        end do
    end subroutine

    subroutine IntrinsicMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))
        mtx3=matmul(mtx1,mtx2)

    end subroutine

    subroutine MixMatrixMultiply(mtx1,mtx2,mtx3)
        real,intent(in)                 ::  mtx1(:,:),mtx2(:,:)
        real,intent(out),allocatable    ::  mtx3(:,:)
        integer                         ::  m,n
        integer                         ::  i,j
        if(size(mtx1,dim=2) /= size(mtx2,dim=1)) stop "input array size not match"
        m=size(mtx1,dim=1)
        n=size(mtx2,dim=2)
        allocate(mtx3(m,n))

        do i=1,m
            do j=1,n
                    mtx3(i,j)=dot_product(mtx1(i,:),mtx2(:,j))
            end do
        end do

    end subroutine

end module




program main
use MatrixMultiplyModule
implicit none

real,allocatable        ::  a(:,:),b(:,:)
real,allocatable    ::  c1(:,:),c2(:,:),c3(:,:)
integer ::  n
integer ::  count, rate
real    ::  timeAtStart, timeAtEnd
real    ::  time(3,10)
do n=100,1000,100
    allocate(a(n,n),b(n,n))

    call random_number(a)
    call random_number(b)

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call LoopMatrixMultiply(a,b,c1)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(1,n/100)=timeAtEnd-timeAtStart

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call IntrinsicMatrixMultiply(a,b,c2)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(2,n/100)=timeAtEnd-timeAtStart

    call system_clock(count = count, count_rate = rate)
    timeAtStart = count / real(rate)
    call MixMatrixMultiply(a,b,c3)
    call system_clock(count = count, count_rate = rate)
    timeAtEnd = count / real(rate)
    time(3,n/100)=timeAtEnd-timeAtStart


    deallocate(a,b)

end do

open(1,file="time.txt")
do n=1,10
    write(1,*) time(:,n)
end do
close(1)
deallocate(c1,c2,c3)
end program

Answer 1

循环遍历数组元素时应该注意几件事:

确保内部循环位于内存中的连续元素之上。在您当前的“循环”算法中，内部循环超过索引 k。由于矩阵在内存中作为列布局(第一个索引在通过内存时变化最快)，访问新的 k 值可能需要将新页面加载到缓存中。在这种情况下，您可以通过将循环重新排序为:

do j=1,n
    do k=1,size(mtx1,dim=2)
        do i=1,m
            mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
        end do      
    end do
end do

现在，内部循环在内存中的连续元素上(mtx2(k,j) 值可能只会在内部循环之前被编译器获取一次，如果不是，您可以将它存储在循环之前的临时变量中)

确保整个循环可以放入缓存中，以避免过多的缓存未命中。这可以通过阻塞算法来完成。在这种情况下，解决方案可能是例如:

l=size(mtx1,dim=2)
ichunk=512 ! I have a 3MB cache size (real*4)
do jj=1,n,ichunk
    do kk=1,l,ichunk

do j=jj,min(jj+ichunk-1,n)
    do k=kk,min(kk+ichunk-1,l)
        do i=1,m
            mtx3(i,j)=mtx3(i,j)+mtx1(i,k)*mtx2(k,j)
        end do      
    end do
end do

    end do
end do

在这种情况下，性能将取决于 ichunk 的大小，尤其是对于足够大的矩阵(您甚至可以阻止内部循环，这只是一个示例)。

确保执行循环所需的工作远小于循环内部的工作。这可以通过“循环展开”来解决，即在循环的一次迭代中组合多个语句。通常编译器可以通过提供标志 -funroll-loops 来做到这一点。 .

如果我使用上面的代码并使用标志 -O3 -funroll-loops 进行编译, 我的性能比 matmul 稍微好一点.

这三个要记住的重要一点是关于循环排序的第一点，因为这会影响其他用例的性能，而编译器通常无法修复它。循环展开，您可以留给编译器(但对其进行测试，因为这并不总能提高性能)。至于第二点，由于这取决于硬件，因此您不应该(通常)尝试自己实现一个非常有效的矩阵乘法，而是考虑使用诸如 e.g. 之类的库。 atlas，可以优化缓存大小，或供应商库，如 MKL 或 ACML。