关于几个数学工具,矩阵运算速度的初步测试
本帖最后由 gps99 于 2022-4-20 19:49 编辑今天闲的蛋疼,测试了几个数学工具, 关于矩阵运算的大致效率。
测试内容: 随机生成一个大矩阵A(分别取1k*1k,2k*2k,3k*3k,4k*4k);取Ai=A的逆阵;再 A * Ai 应该 =E单位阵。
软件/语言: Maple,Matlab,numpy,Julia, Gfortran.
源码:已表中列出(除fortran外,均为两三行的简单代码)
- Fortran代码为
module inv
contains
! decomposition.Depends on LAPACK.
function inverse(A) result(Ainv)
real(kind=8), dimension(:,:), intent(in) :: A
real(kind=8), dimension(size(A,1),size(A,2)) :: Ainv
real(kind=8), dimension(size(A,1)) :: work! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info
!external DGETRF
!external DGETRI
Ainv = A
n = size(A,1)
call DGETRF(n, n, Ainv, n, ipiv, info)
if (info /= 0) then
stop 'Matrix is numerically singular!'
end if
call DGETRI(n, Ainv, n, ipiv, work, n, info)
if (info /= 0) then
stop 'Matrix inversion failed!'
end if
end function inverse
end module
program MatrixSpeed
use inv
implicit none
real(kind=8), allocatable :: z(:,:), re(:,:)
integer :: i, j, n, p
real(kind=8) :: u
real(kind=4) :: t1, t2
n = 4000
p = 4000
allocate(z(n,p),re(n,n))
do i=1,n
do j=1,p
call random_number(u)
z(i,j) = u * 3.0
enddo
enddo
! 随机矩阵求逆,再乘以原矩阵,得到单位阵。
!external dgemm
print *, "正在计算... ..."
call cpu_time(t1)
re = matmul(z, inverse(z))
call cpu_time(t2)
print *, re(3,3)
print *, "Matmul Time:", t2-t1, "秒"
end program MatrixSpeed
测试结果: 如下表
本帖最后由 gps99 于 2022-4-20 19:59 编辑
4种动态代码,基本是两三行调用既有函数。 进一步优化的空间不大。
fortran代码,矩阵的求逆--使用lapack 3.9.1源码 在本地编译的<.a>静态库,矩阵随机生成用Gfor的内置,相乘用内置matmul函数。在Win10下编译为 gfortran -O3。我自己能做到的fortran 也就这样了,,大咖看看能如何优化加快速度?
我目前的结论:
1. matlab的向量速度确实可以。在避免大规模for迭代的场景下,matlab绝对不算太慢。
2. maple 专注于符号计算,数值计算尤其是向量化规模计算,直接出局了。可能它根本不是向量化算法。(它的算法很特殊, 对角数据是整数1。 是的,没错,不是浮点数0.99999999999xx的形式,直接是 1。与所有竞品不同,它可能是某种符号算法。。。比较NB)
3. Numpy 的 C实现,确实够快!!可以的
4. Julia当仁不让,与numpy持平,速度相当惊艳! 当然它的确需要一些技巧的加持。如果数组和py一样未定义类型--Array{Float64} /或 float32,julia则无法表现出传说中的优秀,将拖慢20倍左右。。。 [此处 2022-04-20修改]
5. Fortran。。。小有郁闷,,废这么大劲儿写代码,居然没速度优势? 纳尼? 应该是我太水。。。
请大咖 就这个程序,给出您的优化。
优化代码也好; 调用更NB的库也好;用其他编译器也好。。。 fortran不应该比matlab慢才对。
gps99 发表于 2022-4-3 22:26
4种动态代码,基本是两三行调用既有函数。 进一步优化的空间不大。
fortran代码,矩阵的求逆--使用lapack...
矩阵运算中,Matlab和Numpy的运算内核用的都是blas, Lapack,可能是C也可能是Fortran。如果所用的库函数是一样的,那么最可能出问题的是内置matmul函数。Gfortran的matmul函数是比较慢的。改用GEMM或许有不一样的结果。 做矩阵整体运算,fortran 不一定比matlab快,matlab在大矩阵的时候是可以自动并行的。另外numpy和matlab内部用的都是mkl或者openblas的库,相比之下优势更为明显。所以总体来说,速度的差距都在调库上。 风平老涡 发表于 2022-4-4 00:03
矩阵运算中,Matlab和Numpy的运算内核用的都是blas, Lapack,可能是C也可能是Fortran。如果所用的库函数 ...
试一下dgemm,好像速度比matmul还慢一点。。(另外dgemm结果有点异常,可能是我没用对)
Transpose 发表于 2022-4-4 09:52
做矩阵整体运算,fortran 不一定比matlab快,matlab在大矩阵的时候是可以自动并行的。另外numpy和matlab内 ...
嗯嗯。某些动态语言除了一贯的便利性上,性能也在逐步改善。numpy,numba,julia。。
据说Ifortran也在预备中,可Jupyter中交互运行。很好奇,这将开启一个fortran的什么姿势。。 本帖最后由 gps99 于 2022-4-4 11:47 编辑
刚试了下Lapack的 dgemm计算矩阵乘法,有异常。请教大家。
dgemm的逆阵*原阵 ≠ 单位阵E。(内置的Matmul 函数则结果正确)。对比源码如下:
module inv
contains
! decomposition.Depends on LAPACK.
function inverse(A) result(Ainv)
real(kind=8), dimension(:,:), intent(in) :: A
real(kind=8), dimension(size(A,1),size(A,2)) :: Ainv
real(kind=8), dimension(size(A,1)) :: work! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info
!external DGETRF
!external DGETRI
Ainv = A
n = size(A,1)
call DGETRF(n, n, Ainv, n, ipiv, info)
if (info /= 0) then
stop 'Matrix is numerically singular!'
end if
call DGETRI(n, Ainv, n, ipiv, work, n, info)
if (info /= 0) then
stop 'Matrix inversion failed!'
end if
end function inverse
end module
program MatrixSpeed
use inv
implicit none
real(kind=8), allocatable :: z(:,:), re(:,:)
integer :: i, j, n, p
real(kind=8) :: u
real(kind=8) :: RealOne = 1.d0, RealZero = 0.d0
real(kind=4) :: t1, t2 ,t3
n = 1000
p = 1000
allocate(z(n,p),re(n,n))
do i=1,n
do j=1,p
call random_number(u)
z(i,j) = u * 3.0
enddo
enddo
! 随机矩阵求逆,再乘以原矩阵,得到单位阵。
print *, "正在计算... ..."
! 1. Matmul 实现乘法
call cpu_time(t1)
re = matmul(z, inverse(z))
call cpu_time(t2)
print *, re(3,3)
do i=1,min(10,n)
write(*,'(10(f10.3,x))') re(i, 1:min(10,n))
enddo
print *, "Matmul Time:", t2-t1, "秒"
! 2. Lapack-dgemm 实现乘法
!external dgemm
call dgemm('n','t' ,n,n,p,RealOne ,Z,n ,inverse(Z),n ,RealZero,RE,n)
call cpu_time(t3)
print *, re(3,3)
do i=1,min(10,n)
write(*,'(10(f10.3,x))') re(i, 1:min(10,n))
enddo
print *, "DgemmTime:", t3-t2, "秒"
end program MatrixSpeed
! gfortran -O3 01_array_inv_mul.f90 -L./lapack -llapack -lrefblas -o array_inv_mul.exe
matmul 得到单位阵; dgemm按netlib.org的介绍 (Dgemm - 说明页面 ),, 查了好几遍也没弄对结果。。
请教大咖,问题出在哪里??
```
call dgemm("N","N",...)
```
不需要转置。另外,lapack编译出来的refblas 速度上不占优势,高度优化的库函数基本都是要用汇编的。 Transpose 发表于 2022-4-4 11:44
```
call dgemm("N","N",...)
```
多谢大咖指点:-handshake 药到病除,立竿见影!!
本帖最后由 gps99 于 2022-4-20 20:01 编辑
比较一轮下来,弄个彻底的,加一项: “原生态的理论代码”(取自 徐士良 fortran常用程序代码)。先说结果吧。。。比Numpy慢150~200倍,比Lapack慢15~20倍。
.
代码:矩阵求逆部分代码,选自 徐士良 ”Forrtan常用算法程序集“。矩阵相乘部分 仍使用内置matmul保持一致。编译 Gfortran -O3
说明:逻辑部分没有任何改动。将原书的代码改为f90的Modern结构。 与Lapack的15倍差距,基本可以认为是 for循环矩阵 vs. 向量化矩阵运算的差距
-- 测试源码如下:
module inv
contains
subroutine inverse(a, n ,l)
implicit none
real(kind=8) ::a(n,n)
integer(kind=4) ::is(n), js(n)
real(kind=8) ::t, d
integer ::i, j, k, n, l
l=1
do k=1, n
d=0.0
do i=k, n
do j=k, n
if (abs(a(i,j)).gt.d) then
d=abs(a(i,j))
is(k)=i
js(k)=j
end if
end do
end do
if (d+1.0.eq.1.0) then
l=0
write(*, '(1x, "Err ** not inv")')
return
end if
do j=1, n
t = a(k, j)
a(k, j) = a(is(k),j)
a(is(k),j) = t
end do
do i=1, n
t=a(i,k)
a(i,k)=a(i,js(k))
a(i,js(k))=t
end do
a(k,k)=1/a(k,k)
do j=1,n
if (j.ne.k) then
a(k,j)=a(k,j)*a(k,k)
end if
end do
do i = 1, n
if (i.ne.k) then
do j=1, n
if (j.ne.k) then
a(i,j)=a(i,j)-a(i,k)*a(k,j)
end if
end do
end if
end do
do i=1,n
if (i.ne.k) then
a(i,k)=-a(i,k)*a(k,k)
end if
end do
end do
do k=n,1,-1
do j=1,n
t=a(k,j)
a(k,j)=a(js(k),j)
a(js(k),j)=t
end do
do i=1,n
t=a(i,k)
a(i,k)=a(i,is(k))
a(i,is(k))=t
end do
end do
return
end subroutine inverse
end module inv
program main
use inv
implicit none
real(kind=8), allocatable :: a(:,:),b(:,:),re(:,:)
integer :: i, j, n, l
real(kind=8) :: u
real(kind=4) :: t1, t2
n = 4000
allocate(a(n,n),b(n,n))
do i=1,n
do j=1,n
call random_number(u)
a(i,j) = u * 3.0
enddo
enddo
b=a
! 随机矩阵求逆,再乘以原矩阵,得到单位阵。
print *, "N=", n
print *, "正在计算... ..."
call cpu_time(t1)
call inverse(a, n, l)
re = matmul(a, b)
call cpu_time(t2)
do i=1,min(10,n)
write(*,'(10(f10.3,x))') re(i, 1:min(10,n))
enddo
print *, "Matmul Time:", t2-t1, "秒"
end program main
! gfortran -O3 array_inv_mul_original.f90 -L./lapack -llapack -lrefblas -o array_inv_mul_original.exe
!
测试结果: 骇人听闻,惨不忍睹。。
页:
[1]
2