我正在阅读“使用 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
最佳答案
循环遍历数组元素时应该注意几件事:
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。
关于performance - 不同优化下的 Fortran 矩阵乘法性能,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/15580572/