Fortran Coder
标题: 关于几个数学工具，矩阵运算速度的初步测试 [打印本页]
作者: gps99 时间: 2022-4-3 21:37
标题: 关于几个数学工具，矩阵运算速度的初步测试
本帖最后由 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代码为[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-3 22:26
本帖最后由 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慢才对。
作者: 风平老涡 时间: 2022-4-4 00:03
矩阵运算中，Matlab和Numpy的运算内核用的都是blas, Lapack，可能是C也可能是Fortran。如果所用的库函数是一样的，那么最可能出问题的是内置matmul函数。Gfortran的matmul函数是比较慢的。改用GEMM或许有不一样的结果。
作者: Transpose 时间: 2022-4-4 09:52
做矩阵整体运算，fortran 不一定比matlab快，matlab在大矩阵的时候是可以自动并行的。另外numpy和matlab内部用的都是mkl或者openblas的库，相比之下优势更为明显。所以总体来说，速度的差距都在调库上。
作者: gps99 时间: 2022-4-4 10:42
试一下dgemm，好像速度比matmul还慢一点。。（另外dgemm结果有点异常，可能是我没用对）
作者: gps99 时间: 2022-4-4 11:17
嗯嗯。 某些动态语言除了一贯的便利性上，性能也在逐步改善。numpy，numba，julia。。
据说Ifortran也在预备中，可Jupyter中交互运行。很好奇，这将开启一个fortran的什么姿势。。
作者: gps99 时间: 2022-4-4 11:35
本帖最后由 gps99 于 2022-4-4 11:47 编辑
刚试了下Lapack的 dgemm计算矩阵乘法，有异常。请教大家。
dgemm的逆阵*原阵 ≠ 单位阵E。（内置的Matmul 函数则结果正确）。对比源码如下：
[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=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 *, "Dgemm Time: ", 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 - 说明页面 )，, 查了好几遍也没弄对结果。。
请教大咖，问题出在哪里？？
作者: Transpose 时间: 2022-4-4 11:44
```
call dgemm("N","N",...)
```
不需要转置。另外，lapack编译出来的refblas 速度上不占优势，高度优化的库函数基本都是要用汇编的。
作者: gps99 时间: 2022-4-4 11:54
多谢大咖指点
药到病除，立竿见影！！
作者: gps99 时间: 2022-4-4 18:38
本帖最后由 gps99 于 2022-4-20 20:01 编辑
比较一轮下来，弄个彻底的，加一项： “原生态的理论代码”（取自 徐士良 fortran常用程序代码）。先说结果吧。。。比Numpy慢150~200倍，比Lapack慢15~20倍。
.
代码： 矩阵求逆部分代码，选自 徐士良 ”Forrtan常用算法程序集“。矩阵相乘部分 仍使用内置matmul保持一致。编译 Gfortran -O3
说明： 逻辑部分没有任何改动。将原书的代码改为f90的Modern结构。 与Lapack的15倍差距，基本可以认为是 for循环矩阵 vs. 向量化矩阵运算的差距
-- 测试源码如下：
[Fortran] 纯文本查看 复制代码
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
!
测试结果： 骇人听闻，惨不忍睹。。
作者: gps99 时间: 2022-4-4 18:49
本帖最后由 gps99 于 2022-4-4 19:00 编辑
喔噻，，我居然有耐心等这个4k的结果出来 溜溜的跑了十几分钟。、。
是的，十几分钟一致在想，这也许是今后 自己手写数值计算的下场。
著名前辈尚且如此，我等千万不要自己造轮子！！！千万千万！！
作者: zoziha 时间: 2022-4-4 20:18
本帖最后由 zoziha 于 2022-4-4 20:28 编辑
你是使用纯Fortran的Lapack吗？从我的认知上看，Numpy后端矩阵求解使用的好像是OpenBLAS，在Fortran中也能方便地使用OpenBLAS，理论上Fortran写法(OpenBLAS)和Numpy写法性能表现应该一致。
Python代码越少，OpenBLAS调用相关代码越多，Python的Numpy写法效率应该更高。
同样的Intel的Fortran编译器使用自家的MKL性能应该也是飞起。（据我所知，matlab使用的是intel家的MKL库，这谁也打不过啊，intel家软硬件结合，cpu届应该无敌）
自实现算法与使用著名函数库（如OpenBLAS、MKL、IMSL、Petsc、Numpy、fftw等），后者往往是许多人集体工作的结晶，也普遍存在大量优化，是典型的商业或社区驱动“产品”，而自实现算法往往只能在代码层面优化，难以完美掌控底层编译器行为，性能一般是不如著名函数库。但是，这也是跟着具体的业务和需求走的，不需要长时间大量计算的代码，也不一定极致追求调用性能最好的库，盲目引入外部库也可能是一种不好的实践。
互联网时代，Python自带优势，借着这东风，已经很强了，我能看到很多c/cpp/fortran代码被Python社区（从scipy/numpy源码中可见）很好地利用起来了，但编译型语言始终代表性能top1行列。希望Fortran能迎头赶上。
作者: gps99 时间: 2022-4-4 23:12
是的，纯fortran的Lapack。 在 netlib-Lapack 3.9.1 下载源码，按包内的makefile本地编译得到三个.a静态库.
从几个软件所各自代表的LinearAlgabra包角度看，实测的差距感觉比预想的大。 纯矩阵效率上 自实现代码 <- 专用包Lapack... <- 商业专业包MKL Matlab包... 三个台阶泾渭分明，居然能有数量级的差距 的确很意外。当然实际场景中综合各种运算，不会有这类白手套测试的差距明显，正如您所说 盲目引入外部包不一定是明智的选择，结构需综合考虑。
另外今天还大致跑了其他几个测试，包括纯for迭代的整形、浮点速度，fibnaci的递归效率， 。。 有一个明显的感觉就是：编译语言的各项性能比较稳定（gcc，gfortran .VS. matlab， python+numba） 。在递归、for迭代、整形、浮点的效率上，c/fortran比动态语言有几十倍的压倒性优势。综合性能上 尤其是C，绝对实至名归，一个字 稳。
作者: 唐汉 时间: 2022-4-17 07:49
Julia那个应该是把时间花在类型转换上了。类似的“benchmark”网上有很多，但是大多数都是疑点重重的，关键在于怎么用不同的语言写出等效的代码。比如说你的julia代码，M乘以Mi的时候是一个整数矩阵乘以一个浮点数矩阵，每一个元素相乘的时候都要做一次类型转换。但是如果你提前把M转化成浮点数的矩阵，速度就上来了。python或者是numpy应该是把这一步给包起来了。我测试了一下
（1）Python：
[Python] 纯文本查看 复制代码
import numpy as np
matrix_size = 4000
int_range = (1, 99)
M = np.random.randint(*int_range, (matrix_size, matrix_size))
Mi = np.linalg.inv(M)
Re = np.dot(M, Mi)
（2）Julia：（论坛里貌似没有Julia的高亮）
[Python] 纯文本查看 复制代码
matrix_size = (4000, 4000)
int_range = 1:99
M = rand(int_range, matrix_size...)
Mᵣ = Array{Float32}(M)
Mᵢ = inv(Mᵣ)
Rₑ = Mᵣ * Mᵢ
结果是差不多的julia稍微快一点：
$ time python test.py
real 0m3.473s
user 0m22.469s
sys 0m2.250s
$time julia test.jl
real 0m2.202s
user 0m8.245s
sys 0m2.246s
如果换成Float64 julia速度会降下来一点但是也就是4秒多一点，绝对不至于跌出一个数量级来。然后调用OpenBLAS的时候矩阵超过一定大小应该就会自动超线程了，Julia和python在这个层面上是否做了一样的事情我也不太清楚。
大部分情况下如果代码耗时在同一个硬件上出现了数量级的差距，那基本上都是算法上不等效（当然你如果用纯python不调用任何库那肯定会慢出天际……）。
Fortran那个代码我暂时还没测试，有空了我再来玩玩~
作者: gps99 时间: 2022-4-20 19:25
本帖最后由 gps99 于 2022-4-20 19:28 编辑
感谢您对Julia的指导。
前面的Jl代码没有指定类型，确实影响了julia的测试速度。指定Array{Float64}后，Julia速度提高近20倍。核Numpy基本一致，甚至略快。。 果然Julia是需要注入技巧，才有速度。
同是 float64格式下：
Julia 3.3 s
Numpy 3.4 s
Julia 代码如下：
[Fortran] 纯文本查看 复制代码
@time begin
global M = rand(1:99, 4000, 4000)
global M = Array{Float64}(M)
Mi = inv(M)
Re = M * Mi
println(Re[11, 11])
end
作者: gps99 时间: 2022-4-20 20:16
本帖最后由 gps99 于 2022-4-20 20:20 编辑
Numpy 和 Julia在纯向量的优化程度上，如此的接近，应该不是巧合。 也许是目前商业库能达到的某种特征水平。
过几天 再把openBLAS试试，看看是一个什么结果。。。估计达不到4 s以内。
论坛内大咖云集，哪位有时间的话 c++或fortran调一下openBLAS，并优化一下。。。 。。。。最后放一个EXE的附件。
测试目的只有一个： 即便是一个看似很简单的调用， 我们能否实现商业库的效率？？尽力追，，能追到一个什么程度？
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) |
Powered by Discuz! X3.2 |