这是之前关于 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 programPython+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.exePBS 文件示例（Python）：
#!/bin/sh
#PBS -N python
#PBS -l nodes=1:ppn=24
python a.py以下分别是在三个设备中用 24 核 CPU 运行的结果，仅供参考。
Fortran+MKL 24 核运行结果（设备 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 核运行结果（设备 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 核运行结果（设备 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 核运行结果（设备 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 核运行结果（设备 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 核运行结果（设备 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 单核运行结果（设备 H）：
n = 100 的计算时间: 0.037915 秒
n = 200 的计算时间: 0.001645 秒
n = 300 的计算时间: 0.003420 秒
n = 400 的计算时间: 0.006320 秒
n = 500 的计算时间: 0.011255 秒
n = 600 的计算时间: 0.017210 秒
n = 700 的计算时间: 0.027165 秒
n = 800 的计算时间: 0.039520 秒
n = 900 的计算时间: 0.055745 秒
n = 1000 的计算时间: 0.078775 秒
n = 2000 的计算时间: 0.522050 秒
n = 3000 的计算时间: 1.459600 秒
n = 4000 的计算时间: 3.026015 秒
n = 5000 的计算时间: 5.068755 秒
n = 6000 的计算时间: 8.457271 秒
n = 7000 的计算时间: 13.082525 秒
n = 8000 的计算时间: 19.187775 秒
n = 9000 的计算时间: 26.694241 秒
n = 10000 的计算时间: 36.200172 秒
因为设备 H 对单核的内存有限制，1 核绑定为 4G 内存，所以这里跳过 2w 和 3w 维度的矩阵的测试。Fortran+MKL 单核运行结果（设备 S）：
n = 100 的计算时间: 0.023145 秒
n = 200 的计算时间: 0.046640 秒
n = 300 的计算时间: 0.069955 秒
n = 400 的计算时间: 0.134025 秒
n = 500 的计算时间: 0.148670 秒
n = 600 的计算时间: 0.307725 秒
n = 700 的计算时间: 0.992025 秒
n = 800 的计算时间: 1.187035 秒
n = 900 的计算时间: 1.524415 秒
n = 1000 的计算时间: 1.321075 秒
n = 2000 的计算时间: 1.948630 秒
n = 3000 的计算时间: 4.112365 秒
n = 4000 的计算时间: 6.514895 秒
n = 5000 的计算时间: 6.258010 秒
n = 6000 的计算时间: 8.316165 秒
n = 7000 的计算时间: 11.349990 秒
n = 8000 的计算时间: 15.349805 秒
n = 9000 的计算时间: 19.495274 秒
n = 10000 的计算时间: 25.048605 秒
n = 20000 的计算时间: 126.816452 秒
n = 30000 的计算时间: 554.398315 秒Fortran+MKL 单核运行结果（设备 P）：
n = 100 的计算时间: 0.001630 秒
n = 200 的计算时间: 0.001505 秒
n = 300 的计算时间: 0.003160 秒
n = 400 的计算时间: 0.005850 秒
n = 500 的计算时间: 0.010475 秒
n = 600 的计算时间: 0.016065 秒
n = 700 的计算时间: 0.024410 秒
n = 800 的计算时间: 0.034770 秒
n = 900 的计算时间: 0.048545 秒
n = 1000 的计算时间: 0.062055 秒
n = 2000 的计算时间: 0.407050 秒
n = 3000 的计算时间: 1.346525 秒
n = 4000 的计算时间: 3.166775 秒
n = 5000 的计算时间: 6.064680 秒
n = 6000 的计算时间: 10.295170 秒
n = 7000 的计算时间: 16.081350 秒
n = 8000 的计算时间: 23.533381 秒
n = 9000 的计算时间: 33.309010 秒
n = 10000 的计算时间: 45.399574 秒
n = 20000 的计算时间: 450.167816 秒
n = 30000 的计算时间: 1491.506836 秒Python+NumPy 单核运行结果（设备 H）：
n = 100 的计算时间: 0.020551 秒
n = 200 的计算时间: 0.001476 秒
n = 300 的计算时间: 0.005053 秒
n = 400 的计算时间: 0.006819 秒
n = 500 的计算时间: 0.013365 秒
n = 600 的计算时间: 0.021095 秒
n = 700 的计算时间: 0.030566 秒
n = 800 的计算时间: 0.043758 秒
n = 900 的计算时间: 0.058170 秒
n = 1000 的计算时间: 0.074706 秒
n = 2000 的计算时间: 0.550421 秒
n = 3000 的计算时间: 1.589639 秒
n = 4000 的计算时间: 3.597673 秒
n = 5000 的计算时间: 6.017714 秒
n = 6000 的计算时间: 10.068739 秒
n = 7000 的计算时间: 16.116104 秒
n = 8000 的计算时间: 23.762686 秒
n = 9000 的计算时间: 32.594450 秒
n = 10000 的计算时间: 44.582610 秒
因为设备 H 对单核的内存有限制，1 核绑定为 4G 内存，所以这里跳过 2w 和 3w 维度的矩阵的测试。Python+NumPy 单核运行结果（设备 S）：
n = 100 的计算时间: 0.027154 秒
n = 200 的计算时间: 0.082743 秒
n = 300 的计算时间: 0.081564 秒
n = 400 的计算时间: 0.109224 秒
n = 500 的计算时间: 0.109470 秒
n = 600 的计算时间: 0.428341 秒
n = 700 的计算时间: 1.012723 秒
n = 800 的计算时间: 1.089292 秒
n = 900 的计算时间: 1.360963 秒
n = 1000 的计算时间: 1.302347 秒
n = 2000 的计算时间: 2.407817 秒
n = 3000 的计算时间: 4.414838 秒
n = 4000 的计算时间: 5.896137 秒
n = 5000 的计算时间: 6.105954 秒
n = 6000 的计算时间: 8.344055 秒
n = 7000 的计算时间: 11.416265 秒
n = 8000 的计算时间: 15.638111 秒
n = 9000 的计算时间: 20.103910 秒
n = 10000 的计算时间: 25.833898 秒
n = 20000 的计算时间: 157.841558 秒
n = 30000 的计算时间: 807.751309 秒Python+NumPy 单核运行结果（设备 P）：
n = 100 的计算时间: 0.001823 秒
n = 200 的计算时间: 0.002341 秒
n = 300 的计算时间: 0.004830 秒
n = 400 的计算时间: 0.009922 秒
n = 500 的计算时间: 0.018024 秒
n = 600 的计算时间: 0.020881 秒
n = 700 的计算时间: 0.031204 秒
n = 800 的计算时间: 0.044493 秒
n = 900 的计算时间: 0.061796 秒
n = 1000 的计算时间: 0.079052 秒
n = 2000 的计算时间: 0.598854 秒
n = 3000 的计算时间: 1.876777 秒
n = 4000 的计算时间: 4.259775 秒
n = 5000 的计算时间: 8.068309 秒
n = 6000 的计算时间: 13.616758 秒
n = 7000 的计算时间: 21.326297 秒
n = 8000 的计算时间: 31.155525 秒
n = 9000 的计算时间: 43.874915 秒
n = 10000 的计算时间: 59.834241 秒
n = 20000 的计算时间: 597.668836 秒
n = 30000 的计算时间: 1861.923084 秒【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】
