我正在尝试计算类似于 Fortran 中的加权矩阵内积的东西。我当前用于计算内积的脚本如下
! --> In
real(kind=8), intent(in), dimension(ni, nj, nk, nVar) :: U1, U2
real(kind=8), intent(in), dimension(ni, nj, nk) :: intW
! --> Out
real(kind=8), intent(out) :: innerProd
! --> Local
integer :: ni, nj, nk, nVar, iVar
! --> Computing inner product
do iVar = 1, nVar
innerProd = innerProd + sum(U1(:,:,:,iVar)*U2(:,:,:,iVar)*intW)
enddo
但是我发现我当前使用的上述脚本效率不是很高。可以使用 NumPy 在 Python 中执行相同的操作,如下所示,
import numpy as np
import os
# --> Preventing numpy from multi-threading
os.environ['OPENBLAS_NUM_THREADS'] = '1'
os.environ['MKL_NUM_THREADS'] = '1'
innerProd = 0
# --> Toy matrices
U1 = np.random.random((ni,nj,nk,nVar))
U2 = np.random.random((ni,nj,nk,nVar))
intW = np.random.random((ni,nj,nk))
# --> Reshaping
U1 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
U2 = np.reshape(np.ravel(U1), (ni*nj*nk, nVar))
intW = np.reshape(np.ravel(intW), (ni*nj*nk))
# --> Calculating inner product
for iVar in range(nVar):
innerProd = innerProd + np.dot(U1[:, iVar], U2[:, iVar]*intW)
使用 Numpy 的第二种方法似乎比使用 Fortran 的方法快得多。对于ni = nj = nk = nVar = 130
的具体情况,两种方法所花费的时间如下
fortran_time = 25.8641 s
numpy_time = 6.8924 s
我尝试使用 BLAS 的 ddot
改进我的 Fortran 代码,如下所示,
do iVar = 1, nVar
do k = 1, nk
do j = 1, nj
innerProd = innerProd + ddot(ni, U1(:,j,k,iVar), 1, U2(:,j,k,iVar)*intW(:,j,k), 1)
enddo
enddo
enddo
但时间上并没有明显的改善。对于ni = nj = nk = nVar = 130
的情况,上述方法所花费的时间为~24s
。 (我忘了提及,我使用“-O2”选项编译了 Fortran 代码以优化性能)。
不幸的是,Fortran 中没有用于逐元素矩阵乘法的 BLAS 函数。而且我不想在 Fortran 中使用 reshape,因为与 Fortran 中的 python 不同,reshape 会导致将我的数组复制到新数组,从而导致更多的 RAM 使用。
有没有办法可以加快Fortran的性能,从而接近Numpy的性能?
最佳答案
您可能没有按照您认为的时间进行计时。这是一个完整的 Fortran 示例
program test
use iso_fortran_env, r8 => real64
implicit none
integer, parameter :: ni = 130, nj = 130, nk = 130, nvar = 130
real(r8), allocatable :: u1(:,:,:,:), u2(:,:,:,:), w(:,:,:)
real(r8) :: sum, t0, t1
integer :: i,j,k,n
call cpu_time(t0)
allocate(u1(ni,nj,nk,nvar))
allocate(u2(ni,nj,nk,nvar))
allocate(w(ni,nj,nk))
call cpu_time(t1)
write(*,'("allocation time(s):",es15.5)') t1-t0
call cpu_time(t0)
call random_seed()
call random_number(u1)
call random_number(u2)
call random_number(w)
call cpu_time(t1)
write(*,'("random init time (s):",es15.5)') t1-t0
sum = 0.0_r8
call cpu_time(t0)
do n = 1, nvar
do k = 1, nk
do j = 1, nj
do i = 1, ni
sum = sum + u1(i,j,k,n)*u2(i,j,k,n)*w(i,j,k)
end do
end do
end do
end do
call cpu_time(t1)
write(*,'("Sum:",es15.5," time(s):",es15.5)') sum, t1-t0
end program
输出:
$ gfortran -O2 -o inner_product inner_product.f90
$ time ./inner_product
allocation time(s): 3.00000E-05
random init time (s): 5.73293E+00
Sum: 3.57050E+07 time(s): 5.69066E-01
real 0m6.465s
user 0m4.634s
sys 0m1.798s
在此 Fortran 代码中计算内积不到运行时间的 10%。你如何/什么时间安排非常重要。你确定你在 fortran 和 python 版本中计时相同的事情吗?您确定只对inner_product 计算进行计时吗?
关于fortran - 如何在Fortran中高效计算矩阵内积?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/49718256/