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代码为
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 函数则结果正确)。对比源码如下:


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 *, "DgemmTime:", 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",...)
```


多谢大咖指点:-handshake   药到病除,立竿见影!!

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. 向量化矩阵运算的差距

            -- 测试源码如下:


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
!


测试结果:   骇人听闻,惨不忍睹。。








页: [1] 2
查看完整版本: 关于几个数学工具,矩阵运算速度的初步测试