语言, Fortran专题

Fortran+MKL和Python+NumPy矩阵求逆的时间对比

发布时间：

最近修改时间：2025年4月4日

这是之前关于 Fortran 的几篇博文：

本篇以矩阵求逆为例，测试 Fortran+MKL 和 Python+NumPy 的计算时间。

在测试时，需要注意设备的内存大小。如果内存不够大，可以适当减小测试的矩阵维度范围，参考这篇：不同大小的矩阵的内存占用（8B至1TB完整版 ）。

主要结论为：

  • 由于 NumPy 的 numpy.linalg.inv 的底层代码同样是调用 LAPACK 库，因此在计算速度上和 Fortran+MKL 还是比较相近的。在不追求极致计算速度的情况下，Python 完全可以代替 Fortran。
  • Fortran+MKL 的矩阵求逆速度一般会比 Python+NumPy 稍微快一些，提升约 20% 左右，大概在 1~2 倍区间，这可能是因为：（1）MKL 的 LAPACK 库会更高效一些，MKL 是 intel 开发的，在 intel 芯片会有更多的性能调优，但 Anaconda 的 NumPy 可能也已经默认链接到 MKL 版本的 LAPACK 库，所以这个情况不一定；（2）Python 是解释型语言，尽管 NumPy 的核心操作由 C 实现，但 Python 层仍需额外开销，导致速度会慢一些。
  • 在有些情况下，Python+NumPy 的计算速度可能会反超 Fortran+MKL 的计算速度，这可能是由于测试样本数量较少、内存带宽限制等原因。另外，也可以在更多不同型号的 CPU 和内存上进行测试等，这里不做更多讨论。

Fortran+MKL 矩阵求逆时间测试代码：

! This code is supported by the website: https://www.guanjihuan.com
! The newest version of this code is on the web page: https://www.guanjihuan.com/archives/45966

module random_matrix_mod
    implicit none
    contains
        subroutine generate_random_matrix(n, A)
            integer, intent(in) :: n
            double precision, allocatable, intent(out) :: A(:,:)
            integer :: ierr
            
            allocate(A(n, n), stat=ierr)
            if (ierr /= 0) stop "内存分配失败"
            
            call init_random_seed()
            call random_number(A)
        end subroutine

        subroutine init_random_seed()
            integer :: i, n, clock, ierr
            integer, allocatable :: seed(:)
            
            call random_seed(size = n)
            allocate(seed(n), stat=ierr)
            if (ierr /= 0) stop "种子分配失败"
            
            call system_clock(count=clock)
            seed = clock + 37 * [(i - 1, i = 1, n)]
            call random_seed(put=seed)
            
            deallocate(seed)
        end subroutine
end module

program main
    use random_matrix_mod
    use f95_precision
    use blas95
    use lapack95
    implicit none

    integer, allocatable :: index1(:)
    integer n, i, j, info, ierr, stage, start, end_val, step, count_start, count_end, count_rate, test_0, test_times
    double precision, allocatable :: A(:,:)
    double precision time_used

    test_times = 20

    ! 定义不同阶段的参数
    do stage = 1, 3
        select case(stage)
        case(1)  ! 第一阶段：100-1000，步长100
            start = 100
            end_val = 1000
            step = 100
        case(2)  ! 第二阶段：2000-10000，步长1000
            start = 2000
            end_val = 10000
            step = 1000
        case(3)  ! 第三阶段：20000-30000，步长10000
            start = 20000
            end_val = 30000
            step = 10000
        end select

        n = start
        do while (n <= end_val)

            allocate(index1(n), stat=ierr)
            
            call system_clock(count_start, count_rate)
            test_0 = 1
            do while (test_0 <= test_times)
                call generate_random_matrix(n, A)
                call getrf(A, index1, info);  call getri(A, index1, info)  ! 使用 getrf 和 getri 对矩阵求逆。这时候 A 不再是原来的矩阵了，而是求逆后的矩阵。
                deallocate(A, stat=ierr)
                test_0 = test_0 + 1
            end do
            call system_clock(count_end)

            ! 打印计算时间
            if (count_rate > 0) then
                time_used = real(count_end - count_start) / real(count_rate) / test_times
                write(*, '(a, I6, a, f12.6, a)') 'n = ', n, ' 的计算时间: ', time_used, ' 秒'
            else
                write(*,*) "无法获取计算时间"
            endif

            deallocate(index1, stat=ierr)

            n = n + step
        end do
    end do

end program

Python+NumPy 矩阵求逆时间测试代码：

"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/45966
"""

import numpy as np
import time

n_array = np.concatenate((np.arange(100, 1000, 100),
                          np.arange(1000, 10000, 1000),
                          np.arange(10000, 40000, 10000)))

for n in n_array:
    test_times = 20
    start_time = time.time()
    for _ in range(test_times):
        A = np.random.rand(n, n)
        inv_A = np.linalg.inv(A)
    inv_time = (time.time() - start_time)/test_times
    print(f"n = {n} 的计算时间: {inv_time:.6f} 秒")

PBS 文件示例（Fortran）：

#!/bin/sh
#PBS -N fortran
#PBS -l nodes=1:ppn=24
./a.exe

PBS 文件示例（Python）：

#!/bin/sh
#PBS -N python
#PBS -l nodes=1:ppn=24
python a.py

以下分别是在三个设备中用 24 核 CPU 运行的结果，仅供参考。

Fortran+MKL 24 核运行结果（设备 1-H）：

n =    100 的计算时间:     0.055525 秒
n =    200 的计算时间:     0.001935 秒
n =    300 的计算时间:     0.006570 秒
n =    400 的计算时间:     0.004825 秒
n =    500 的计算时间:     0.006890 秒
n =    600 的计算时间:     0.009320 秒
n =    700 的计算时间:     0.012285 秒
n =    800 的计算时间:     0.015865 秒
n =    900 的计算时间:     0.019650 秒
n =   1000 的计算时间:     0.024070 秒
n =   2000 的计算时间:     0.096740 秒
n =   3000 的计算时间:     0.239795 秒
n =   4000 的计算时间:     0.448425 秒
n =   5000 的计算时间:     0.706470 秒
n =   6000 的计算时间:     1.066410 秒
n =   7000 的计算时间:     1.518030 秒
n =   8000 的计算时间:     2.048155 秒
n =   9000 的计算时间:     2.782685 秒
n =  10000 的计算时间:     3.637680 秒
n =  20000 的计算时间:    21.123285 秒
n =  30000 的计算时间:    63.932243 秒

Fortran+MKL 24 核运行结果（设备 2-S）：

n =    100 的计算时间:     0.003510 秒
n =    200 的计算时间:     0.002340 秒
n =    300 的计算时间:     0.003930 秒
n =    400 的计算时间:     0.005770 秒
n =    500 的计算时间:     0.008620 秒
n =    600 的计算时间:     0.011105 秒
n =    700 的计算时间:     0.015280 秒
n =    800 的计算时间:     0.019305 秒
n =    900 的计算时间:     0.025700 秒
n =   1000 的计算时间:     0.030405 秒
n =   2000 的计算时间:     0.115055 秒
n =   3000 的计算时间:     0.300335 秒
n =   4000 的计算时间:     0.620170 秒
n =   5000 的计算时间:     1.014275 秒
n =   6000 的计算时间:     1.582420 秒
n =   7000 的计算时间:     2.364375 秒
n =   8000 的计算时间:     3.117580 秒
n =   9000 的计算时间:     4.269075 秒
n =  10000 的计算时间:     5.594840 秒
n =  20000 的计算时间:    34.387581 秒
n =  30000 的计算时间:   102.449661 秒

Fortran+MKL 24 核运行结果（设备 3-P）：

n =    100 的计算时间:     0.016710 秒
n =    200 的计算时间:     0.011890 秒
n =    300 的计算时间:     0.004345 秒
n =    400 的计算时间:     0.005825 秒
n =    500 的计算时间:     0.008130 秒
n =    600 的计算时间:     0.010765 秒
n =    700 的计算时间:     0.014010 秒
n =    800 的计算时间:     0.017120 秒
n =    900 的计算时间:     0.022685 秒
n =   1000 的计算时间:     0.028345 秒
n =   2000 的计算时间:     0.112030 秒
n =   3000 的计算时间:     0.297470 秒
n =   4000 的计算时间:     0.592875 秒
n =   5000 的计算时间:     0.984750 秒
n =   6000 的计算时间:     1.459420 秒
n =   7000 的计算时间:     2.141540 秒
n =   8000 的计算时间:     2.881255 秒
n =   9000 的计算时间:     4.117225 秒
n =  10000 的计算时间:     5.585145 秒
n =  20000 的计算时间:    38.037880 秒
n =  30000 的计算时间:   118.826248 秒

Python+NumPy 24 核运行结果（设备 1-H）：

n = 100 的计算时间: 0.020844 秒
n = 200 的计算时间: 0.015869 秒
n = 300 的计算时间: 0.002238 秒
n = 400 的计算时间: 0.014731 秒
n = 500 的计算时间: 0.053861 秒
n = 600 的计算时间: 0.015732 秒
n = 700 的计算时间: 0.010169 秒
n = 800 的计算时间: 0.013698 秒
n = 900 的计算时间: 0.016606 秒
n = 1000 的计算时间: 0.019759 秒
n = 2000 的计算时间: 0.098382 秒
n = 3000 的计算时间: 0.258957 秒
n = 4000 的计算时间: 0.585829 秒
n = 5000 的计算时间: 0.717663 秒
n = 6000 的计算时间: 1.061896 秒
n = 7000 的计算时间: 1.469291 秒
n = 8000 的计算时间: 2.057066 秒
n = 9000 的计算时间: 2.614083 秒
n = 10000 的计算时间: 3.406435 秒
n = 20000 的计算时间: 20.181018 秒
n = 30000 的计算时间: 64.642707 秒

Python+NumPy 24 核运行结果（设备 2-S）：

n = 100 的计算时间: 0.031161 秒
n = 200 的计算时间: 0.015650 秒
n = 300 的计算时间: 0.013643 秒
n = 400 的计算时间: 0.019081 秒
n = 500 的计算时间: 0.031401 秒
n = 600 的计算时间: 0.027563 秒
n = 700 的计算时间: 0.034102 秒
n = 800 的计算时间: 0.049760 秒
n = 900 的计算时间: 0.061175 秒
n = 1000 的计算时间: 0.072876 秒
n = 2000 的计算时间: 0.242256 秒
n = 3000 的计算时间: 0.557108 秒
n = 4000 的计算时间: 1.090208 秒
n = 5000 的计算时间: 1.459641 秒
n = 6000 的计算时间: 2.157938 秒
n = 7000 的计算时间: 2.701496 秒
n = 8000 的计算时间: 3.478990 秒
n = 9000 的计算时间: 4.626030 秒
n = 10000 的计算时间: 5.833200 秒
n = 20000 的计算时间: 34.834619 秒
n = 30000 的计算时间: 79.427341 秒

Python+NumPy 24 核运行结果（设备 3-P）：

n = 100 的计算时间: 0.003126 秒
n = 200 的计算时间: 0.007001 秒
n = 300 的计算时间: 0.002940 秒
n = 400 的计算时间: 0.004251 秒
n = 500 的计算时间: 0.006071 秒
n = 600 的计算时间: 0.008462 秒
n = 700 的计算时间: 0.012643 秒
n = 800 的计算时间: 0.017243 秒
n = 900 的计算时间: 0.021306 秒
n = 1000 的计算时间: 0.027384 秒
n = 2000 的计算时间: 0.149329 秒
n = 3000 的计算时间: 0.383520 秒
n = 4000 的计算时间: 0.790218 秒
n = 5000 的计算时间: 1.415115 秒
n = 6000 的计算时间: 2.301454 秒
n = 7000 的计算时间: 3.199709 秒
n = 8000 的计算时间: 4.400124 秒
n = 9000 的计算时间: 5.673701 秒
n = 10000 的计算时间: 7.702747 秒
n = 20000 的计算时间: 45.646189 秒
n = 30000 的计算时间: 139.430721 秒

以下是在一个设备中用单核 CPU 运行的结果，仅供参考。

Fortran+MKL 单核运行结果（设备 2-S）：

n =    100 的计算时间:     0.003510 秒
n =    200 的计算时间:     0.004425 秒
n =    300 的计算时间:     0.072140 秒
n =    400 的计算时间:     0.069145 秒
n =    500 的计算时间:     0.039185 秒
n =    600 的计算时间:     0.068060 秒
n =    700 的计算时间:     0.235925 秒
n =    800 的计算时间:     0.254295 秒
n =    900 的计算时间:     0.031515 秒
n =   1000 的计算时间:     0.256685 秒
n =   2000 的计算时间:     0.796440 秒
n =   3000 的计算时间:     1.319735 秒
n =   4000 的计算时间:     1.896235 秒
n =   5000 的计算时间:     2.551130 秒
n =   6000 的计算时间:     3.527985 秒
n =   7000 的计算时间:     4.837920 秒
n =   8000 的计算时间:     5.670895 秒
n =   9000 的计算时间:     8.859410 秒
n =  10000 的计算时间:    11.079590 秒
n =  20000 的计算时间:    60.444641 秒
n =  30000 的计算时间:   175.640289 秒

Python+NumPy 单核运行结果（设备 2-S）：

n = 100 的计算时间: 0.027093 秒
n = 200 的计算时间: 0.064048 秒
n = 300 的计算时间: 0.077069 秒
n = 400 的计算时间: 0.101299 秒
n = 500 的计算时间: 0.081949 秒
n = 600 的计算时间: 0.420397 秒
n = 700 的计算时间: 1.059398 秒
n = 800 的计算时间: 1.278059 秒
n = 900 的计算时间: 1.401319 秒
n = 1000 的计算时间: 1.400217 秒
n = 2000 的计算时间: 2.247926 秒
n = 3000 的计算时间: 4.093585 秒
n = 4000 的计算时间: 5.348192 秒
n = 5000 的计算时间: 7.043575 秒
n = 6000 的计算时间: 11.013543 秒
n = 7000 的计算时间: 15.523776 秒
n = 8000 的计算时间: 21.623775 秒
n = 9000 的计算时间: 28.986586 秒
n = 10000 的计算时间: 40.401018 秒
n = 20000 的计算时间: 268.765458 秒
n = 30000 的计算时间: 611.031390 秒
6 次浏览

【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】

头像

Published by guanjihuan

发表评论

您的邮箱地址不会被公开。 必填项已用 * 标注

Captcha Code