关于几个数学工具，矩阵运算速度的初步测试

gps99

刚试了下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 - 说明页面 )，, 查了好几遍也没弄对结果。。

请教大咖，问题出在哪里？？

gps99

Transpose 发表于 2022-4-4 09:52
做矩阵整体运算，fortran 不一定比matlab快，matlab在大矩阵的时候是可以自动并行的。另外numpy和matlab内 ...

嗯嗯。  某些动态语言除了一贯的便利性上，性能也在逐步改善。numpy，numba，julia。。

据说Ifortran也在预备中，可Jupyter中交互运行。很好奇，这将开启一个fortran的什么姿势。。

gps99

风平老涡 发表于 2022-4-4 00:03
矩阵运算中，Matlab和Numpy的运算内核用的都是blas, Lapack，可能是C也可能是Fortran。如果所用的库函数 ...

试一下dgemm，好像速度比matmul还慢一点。。（另外dgemm结果有点异常，可能是我没用对）

Transpose

做矩阵整体运算，fortran 不一定比matlab快，matlab在大矩阵的时候是可以自动并行的。另外numpy和matlab内部用的都是mkl或者openblas的库，相比之下优势更为明显。所以总体来说，速度的差距都在调库上。

风平老涡

gps99 发表于 2022-4-3 22:26
4种动态代码，基本是两三行调用既有函数。 进一步优化的空间不大。
fortran代码，矩阵的求逆--使用lapack  ...

矩阵运算中，Matlab和Numpy的运算内核用的都是blas, Lapack，可能是C也可能是Fortran。如果所用的库函数是一样的，那么最可能出问题的是内置matmul函数。Gfortran的matmul函数是比较慢的。改用GEMM或许有不一样的结果。

gps99

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慢才对。