如何在 Fortran 中有效地计算矩阵内积？

Question

我正在尝试计算类似于 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，因为与 Python 不同，Fortran 中的 reshape 会导致将我的数组复制到一个新数组，从而导致更多的 RAM 使用。

有没有什么办法可以加快 Fortran 的性能，从而接近 Numpy 的性能？

Answer 1

您可能没有按照自己的想法进行计时。这是一个完整的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 计算？

Answer 2

这样可以避免制作任何副本。 （注意 blas ddot 方法仍然需要为元素产品制作副本）

   subroutine dot3(n,a,b,c,result)
   implicit none
   real(kind=..) a(*),b(*),c(*),result
   integer i,n
   result=0
   do i=1,n
    result=result+a(i)*b(i)*c(i)
   enddo
   end

dot3 是外部的，意味着 not 在模块/包含结构中。 kind 显然应该匹配主声明。

在主代码中：

  innerprod=0
  do iVar = 1, nVar 
  call dot3(ni*nj*nk, U1(1,1,1,iVar),U2(1,1,1,iVar),intW,result)
  innerProd=innerProd+result
  enddo

Answer 3

比较 Numpy 和 Fortran 代码时，我有同样的观察结果。

差异原来是 BLAS 的版本，我发现使用 netlib 中的 DGEMM 类似于循环，并且比 OpenBLAS 慢大约三倍（请参阅this 答案中的配置文件）。

对我来说最令人惊讶的是，OpenBLAS 提供的代码比仅编译 Fortran 三重嵌套循环要快得多。看来这就是 GotoBLAS 的重点，它是用汇编代码为处理器架构手写的。

即使定时正确、正确排序循环、避免复制并使用每个优化标志（在 gfortran 中），性能仍然比 OpenBLAS 慢大约三倍。我没有尝试过 ifort 或 pgi，但我想知道这是否解释了@kvantour 的赞成评论“我的循环在 0.6 秒内完成”（注意，在某些实现中，内在的 matmul 被 BLAS 取代）。