fortran - 大多数编译器会优化MATMUL(TRANSPOSE(A)，B)吗？

在Fortran程序中，我需要计算几个表达式，例如 M · v ， M T· v ， M T· M 等。。
在这里， M 和 v 是尺寸较小(小于100，通常为2-10)的2D和1D数组

我想知道是否可以在编译时将MATMUL(TRANSPOSE(M),v)编写为与MATMUL(N,v)一样有效的某些代码，其中N明确存储为N=TRANSPOSE(M)。我对带有“强”优化标志(例如-O2，-O3或-Ofast)的gnu和ifort编译器特别感兴趣。

最佳答案

在下面，您可以找到几种不同方法的执行时间。

为了使结果尽可能中性，我用一组等效操作的平均时间重新调整了答案。
我忽略了线程。

矩阵时间向量

比较了六个不同的实现:

手册: do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do

matmul: matmul(P,v)

blas N: dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)

matmul-transpose: matmul(transpose(P),v)

matmul-transpose-tmp: Q=transpose(P); w=matmul(Q,v)

blas T: dgemv('T',n,n,1.0D0,P,n,v,1,0,w,1)

在图1和图2中，您可以比较上述情况的时序结果。总的来说，我们不建议在gfortran和ifort中同时使用临时选项。两种编译器都可以更好地优化MATMUL(TRANSPOSE(P),v)。在gfortran中，MATMUL的实现比默认的BLAS更快，ifort清楚地表明mkl-blas更快。

fortran - 大多数编译器会优化MATMUL(TRANSPOSE(A)，B)吗？-LMLPHP

图1:矩阵向量乘法。各种实现的比较都在gfortran上进行。左面板显示绝对时间除以手动计算的大小为1000的系统的总时间。右面板显示绝对时间除以n2×δ。此处δ是大小为1000的手动计算的平均时间除以1000×1000。

图2:矩阵向量乘法。各种实现的比较是在单线程ifort编译上进行的。左面板显示绝对时间除以手动计算的大小为1000的系统的总时间。右面板显示绝对时间除以n2×δ。此处δ是大小为1000的手动计算的平均时间除以1000×1000。

矩阵乘以矩阵

比较了六个不同的实现:

手册: do l=1,n; do j=1,n; do k=1,n; Q(j,l) = P(j,k)*P(k,l); end do; end do; end do

matmul: matmul(P,P)

blas N: dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)

matmul-transpose: matmul(transpose(P),P)

matmul-transpose-tmp: Q=transpose(P); matmul(Q,P)

blas T: dgemm('T','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)

在图3和图4中，您可以比较上述情况的时序结果。与向量情况相反，仅建议将临时情况用于gfortran。在gfortran中，MATMUL的实现比默认的BLAS更快，ifort清楚地表明mkl-blas更快。值得注意的是，手动实现与mkl-blas相当。

图3:矩阵矩阵乘法。各种实现的比较都在gfortran上进行。左面板显示绝对时间除以手动计算的总时间(对于大小为1000的系统)。右面板显示绝对时间除以n3×δ。此处δ是大小为1000的手动计算的平均时间除以1000×1000×1000。

图4:矩阵矩阵乘法。各种实现的比较是在单线程ifort编译上进行的。左面板显示绝对时间除以手动计算的总时间(对于大小为1000的系统)。右面板显示绝对时间除以n3×δ。此处δ是大小为1000的手动计算的平均时间除以1000×1000×1000。

使用的代码:

program matmul_test

  implicit none

  double precision, dimension(:,:), allocatable :: P,Q,R
  double precision, dimension(:), allocatable :: v,w

  integer :: n,i,j,k,l
  double precision,dimension(12) :: t1,t2

  do n = 1,1000
     allocate(P(n,n),Q(n,n), R(n,n), v(n),w(n))
     call random_number(P)
     call random_number(v)

     i=0

     i=i+1
     call cpu_time(t1(i))
     do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     w=matmul(P,v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     w=matmul(transpose(P),v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=transpose(P)
     w=matmul(Q,v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemv('T',n,n,1.0D0,P,n,v,1,0,w,1)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     do l=1,n; do j=1,n; do k=1,n; Q(j,l) = P(j,k)*P(k,l); end do; end do; end do
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=matmul(P,P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=matmul(transpose(P),P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=transpose(P)
     R=matmul(Q,P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemm('T','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
     call cpu_time(t2(i))

     write(*,'(I6,12D25.17)') n, t2-t1
     deallocate(P,Q,R,v,w)
  end do

end program matmul_test

Matmul

fortran - 大多数编译器会优化MATMUL(TRANSPOSE(A)，B)吗？