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

gps99 · 发表于 2022-4-3 21:37:58

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

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

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

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

Transpose · 发表于 2022-4-4 09:52:33

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

gps99 · 发表于 2022-4-4 10:42:29

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

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

gps99 · 发表于 2022-4-4 11:17:10

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

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

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

gps99 · 发表于 2022-4-4 11:35:33

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

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

gps99 · 发表于 2022-4-4 11:54:27

Transpose 发表于 2022-4-4 11:44
```
call dgemm("N","N",...)
```

多谢大咖指点

药到病除，立竿见影！！

gps99 · 发表于 2022-4-4 18:38:59

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

测试结果： 骇人听闻，惨不忍睹。。

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

规矩勋章