确实，Numpy 和 Scipy 有效地调用 LAPACK 例程来执行 numpy.linalg.inv 和 scipy.linalg.inv。
要对一般矩阵求逆，numpy.linalg.inv 求解 A.x=np.eye((n,n))。函数inv() 调用ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)，其中calls call_@lapack_func@(&params); 其中params.B 是恒等矩阵，@lapack_func@ 是sgesv, dgesv, cgesv, zgesv 之一，它们是一般矩阵的线性求解器。
另一方面，scipy.linalg.invcallsgetri，定义为get_lapack_funcs(('getri'),(a1,))。它对应于 lapack 的 DGETRI() 函数，旨在使用 DGETRF() 计算的 LU 分解来计算矩阵的逆。 因此，如果您在 Fortran 中使用 DGETRI()，在 python 中使用 scipy.linalg.inv() 可能会获得类似的性能和结果。
大部分 Lapack 函数都可以使用 scipy.linalg.lapack 调用。这是一个在 cython 模块中使用 scipy.linalg.cython_lapack.dgetri() 的示例：How to compile C extension for Python where C function uses LAPACK library? 这是一个示例代码，比较 scipy.linalg.cython_lapack.dgetrf()+scipy.linalg.cython_lapack.dgetri() 、numpy 和 scipy.linalg .inv() 在 1000x1000 矩阵上：
import numpy as np
from scipy import linalg
import time
import myinverse
n=1000
A=np.random.rand(n,n)
start= time.time()
Am,info,string=myinverse.invert(A.copy())
end= time.time()
print 'DGETRF+DGETRI, ', end-start, ' seconds'
if info==0:
print 'residual ',np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)
else :
print "inversion failed, info=",info, string
start= time.time()
Am=np.linalg.inv(A.copy())
end= time.time()
print 'np.linalg.inv ', end-start, ' seconds'
print 'residual ', np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)
start= time.time()
Am=linalg.inv(A.copy())
end= time.time()
print 'scipy.linalg.inv ', end-start, ' seconds'
print 'residual ',np.linalg.norm(A.dot(Am)-np.identity(n), np.inf)
输出是：
DGETRF+DGETRI, 0.22541308403 seconds
residual 4.2155882951089296e-11
np.linalg.inv 0.29932808876 seconds
residual 4.371813154546711e-11
scipy.linalg.inv 0.298856973648 seconds
residual 9.110997546690758e-11
对于 2000x2000 矩阵：
DGETRF+DGETRI, 1.64830899239 seconds
residual 8.541625644634121e-10
np.linalg.inv 2.02795410156 seconds
residual 7.448244269611659e-10
scipy.linalg.inv 1.61937093735 seconds
residual 1.6453560233026243e-09
LAPACK inversion routine strangely mixes up all variables 中提供了 Fortran 代码链接 DGETRF()+DGETRI()
进行一些更改后，让我们运行：
PROGRAM solvelinear
implicit none
REAL(8), dimension(1000,1000) :: A,Ainv,M,LU
REAL(8),allocatable :: work(:)
REAL(8) :: wwork
INTEGER :: info,lwork
INTEGER,dimension(1000) :: ipiv
INTEGER :: i,j
real :: start, finish
! put code to test here
info=0
!work=0
ipiv=0
call RANDOM_NUMBER(A)
call cpu_time(start)
!-- LU factorisation
LU = A
CALL DGETRF(1000,1000,LU,1000,ipiv,info)
!-- Inversion of matrix A using the LU
Ainv=LU
lwork=-1
CALL DGETRI(1000,Ainv,1000,Ipiv,wwork,lwork,info)
lwork =INT( wwork+0.1)
allocate(work(lwork))
CALL DGETRI(1000,Ainv,1000,Ipiv,work,lwork,info)
deallocate(work)
call cpu_time(finish)
print '("Time = ",f6.3," seconds.")',finish-start
!-- computation of A^-1 * A to check the inverse
M = matmul(Ainv,A)
print*,"M = "
do i=1,3
do j=1,3
print*,M(i,j)
enddo
end do
END PROGRAM solvelinear
使用gfortran main2.f90 -o main2 -llapack -lblas -lm -Wall编译后，1000x1000矩阵需要0.42s，2000x2000矩阵需要3s。
最后，如果 Fortran 代码和 python 代码未链接到相同的 Blas/Lapack 库，则可能会出现不同的性能。 要调查此问题，请键入 np.__config__.show() 之类的命令，如 Link ATLAS/MKL to an installed Numpy 所示或How to check BLAS/LAPACK linkage in NumPy and SciPy? 中报告的命令。
为了进一步利用分布式计算，petsc 不鼓励反转完整矩阵，因为它很少需要。还写到MatMatSolve(A,B,X)，其中B 和X 是稠密矩阵，可用于执行此操作。此外，此函数在 python 接口petsc4py 中作为对象petsc4py.PETSc.Mat 的方法matSolve(self, Mat B, Mat X) 提供。 no-maintained-anymore Elemental library 被列为实现密集矩阵的直接求解器。虽然 Elemental 库支持 python 接口，但它的 fork Hydrogen 不再支持它。
尽管如此，Elemental 页面列出了一些与分布式密集线性代数相关的开源项目。 ScaLapack 提供了例程 PDGETRI()/PZGETRI() 以使用 LU 分解来反转分布式密集矩阵。这可能会为更快的反转留下一些空间。