gps99 发表于 2022-4-4 18:49:57

本帖最后由 gps99 于 2022-4-4 19:00 编辑

喔噻,,我居然有耐心等这个4k的结果出来 溜溜的跑了十几分钟。、。

是的,十几分钟一致在想,这也许是今后 自己手写数值计算的下场。

著名前辈尚且如此,我等千万不要自己造轮子!!!千万千万!!

zoziha 发表于 2022-4-4 20:18:57

本帖最后由 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:10

zoziha 发表于 2022-4-4 20:18
你是使用纯Fortran的Lapack吗?从我的认知上看,Numpy后端矩阵求解使用的好像是OpenBLAS,在Fortran中也能 ...

是的,纯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:33

Julia那个应该是把时间花在类型转换上了。类似的“benchmark”网上有很多,但是大多数都是疑点重重的,关键在于怎么用不同的语言写出等效的代码。比如说你的julia代码,M乘以Mi的时候是一个整数矩阵乘以一个浮点数矩阵,每一个元素相乘的时候都要做一次类型转换。但是如果你提前把M转化成浮点数的矩阵,速度就上来了。python或者是numpy应该是把这一步给包起来了。我测试了一下

(1)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的高亮)
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:16

本帖最后由 gps99 于 2022-4-20 19:28 编辑

唐汉 发表于 2022-4-17 07:49
Julia那个应该是把时间花在类型转换上了。类似的“benchmark”网上有很多,但是大多数都是疑点重重的,关键 ...
感谢您对Julia的指导。
前面的Jl代码没有指定类型,确实影响了julia的测试速度。指定Array{Float64}后,Julia速度提高近20倍。核Numpy基本一致,甚至略快。。 果然Julia是需要注入技巧,才有速度。
同是 float64格式下:
Julia         3.3 s
Numpy   3.4 s


Julia 代码如下:

@time begin
global M = rand(1:99, 4000, 4000)
global M = Array{Float64}(M)
Mi = inv(M)
Re = M * Mi
println(Re)
end






gps99 发表于 2022-4-20 20:16:28

本帖最后由 gps99 于 2022-4-20 20:20 编辑

Numpy 和 Julia在纯向量的优化程度上,如此的接近,应该不是巧合。 也许是目前商业库能达到的某种特征水平。

过几天 再把openBLAS试试,看看是一个什么结果。。。估计达不到4 s以内。

论坛内大咖云集,哪位有时间的话 c++或fortran调一下openBLAS,并优化一下。。。。。。。最后放一个EXE的附件。

测试目的只有一个:即便是一个看似很简单的调用, 我们能否实现商业库的效率??尽力追,,能追到一个什么程度?
页: 1 [2]
查看完整版本: 关于几个数学工具,矩阵运算速度的初步测试