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

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

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

测试目的只有一个：  即便是一个看似很简单的调用， 我们能否实现商业库的效率？？尽力追，，能追到一个什么程度？

唐汉 发表于 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 代码如下：

[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

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

（1）Python：

[Python] 纯文本查看 复制代码

@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

（2）Julia：（论坛里貌似没有Julia的高亮）

[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)

结果是差不多的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那个代码我暂时还没测试，有空了我再来玩玩~

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，绝对实至名归，一个字 稳。

你是使用纯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能迎头赶上。

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

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

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

```
    call dgemm("N","N",...)
```
不需要转置。另外，lapack编译出来的refblas 速度上不占优势，高度优化的库函数基本都是要用汇编的。