大多数编译器是否优化了MATMUL（TRANSPOSE（A），B）？

Question

在Fortran程序中，我需要计算几个表达式，如M·v，MT·v，MT·M，M·MT等......这里，M和v是2D和1D小尺寸的数组（小于100，通常在2-10左右）

我想知道编写MATMUL(TRANSPOSE(M),v)是否会在编译时展开到像MATMUL(N,v)一样有效的代码，其中N被明确存储为N=TRANSPOSE(M)。我特别感兴趣的是具有“强”优化标志的gnu和ifort编译器（例如-O2，-O3或-Ofast）。

Answer 1

下面是各种方法的几个执行时间。

系统：

• Intel（R）Core（TM）i5-6500T CPU @ 2.50GHz
• 缓存大小：6144 KB
• 内存：16MB
• GNU Fortran（GCC）6.3.1 20170216（Red Hat 6.3.1-3）
• ifort（IFORT）18.0.5 20180823
• BLAS：对于gnu编译器，使用的blas是默认版本

汇编：

[gnu] $ gfortran -O3 x.f90 -lblas
[intel] $ ifort -O3 -mkl x.f90

执行：

[gnu] $ ./a.out > matmul.gnu.txt
[intel] $ EXPORT MKL_NUM_THREADS=1; ./a.out > matmul.intel.txt

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

matrix times vector

比较了六种不同的实现：

1. 手册：do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do
2. MATMUL：matmul(P,v)
3. blas N：dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)
4. matmul-transpose：matmul(transpose(P),v)
5. matmul-transpose-tmp：Q=transpose(P); w=matmul(Q,v)
6. 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更快。

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

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

matrix times matrix

比较了六种不同的实现：

1. 手册：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
2. MATMUL：matmul(P,P)
3. blas N：dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
4. matmul-transpose：matmul(transpose(P),P)
5. matmul-transpose-tmp：Q=transpose(P); matmul(Q,P)
6. blas T：dgemm('T','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)

在图3和图4中，您可以比较上述情况的时序结果。与vector-case相比，只建议使用临时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