fortran - 大多数编译器是否优化 MATMUL(TRANSPOSE(A),B)?
问题描述
在 Fortran 程序中,我需要计算几个表达式,如M · v、M T · v、M T · M、M · M T等......这里,M和v是小尺寸的 2D 和 1D 数组(小于超过 100,通常在 2-10 左右)
我想知道写作MATMUL(TRANSPOSE(M),v)
是否会在编译时展开成一些与 一样有效的代码MATMUL(N,v)
,其中N
显式存储为N=TRANSPOSE(M)
. 我对具有“强”优化标志(例如 -O2、-O3 或 -Ofast)的 gnu 和 ifort 编译器特别感兴趣。
解决方案
您可以在下面找到各种方法的几个执行时间。
系统:
- Intel(R) Core(TM) i5-6500T CPU @ 2.50GHz
- 缓存大小:6144 KB
- 内存:16MB
- GNU Fortran (GCC) 6.3.1 20170216(红帽 6.3.1-3)
- 伊福特 (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
为了使结果尽可能中立,我用完成一组等效操作的平均时间重新调整了答案。我忽略了线程。
矩阵时间向量
比较了六种不同的实现:
- 手动的:
do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do
- 矩阵:
matmul(P,v)
- 布拉斯N:
dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)
- matmul转置:
matmul(transpose(P),v)
- matmul-转置-tmp:
Q=transpose(P); w=matmul(Q,v)
- 爆破:
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。
矩阵乘以矩阵
比较了六种不同的实现:
- 手动的:
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(P,P)
- 布拉斯N:
dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
- matmul转置:
matmul(transpose(P),P)
- matmul-转置-tmp:
Q=transpose(P); matmul(Q,P)
- 爆破:
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
推荐阅读
- compilation - 使用 LINK.EXE 和 C2.EXE 手动编译 VB6 项目
- python - 使用 Python / Apache Beam 进行 Google Cloud Storage 并发控制?
- excel - 代码片段不处理 MS Excel VBA
- mysql - 如何使用触发器更新多个表,从正在更新的表中的行中获取特定数据?
- php - 用当天的结果计算日期
- php - 通过唯一 ID 获取活动资源
- scala - Clojure clojure.core.reducers/fold 和 Scala fold 有什么区别?
- c++ - 不同原始类型的乘法速度
- python - 处理主要损失的缺失数据,辅助损失存在
- eclipse - eclipse没有使用键盘快捷键运行maven构建