使用 python 求解 7000x7000 线性系统时的最佳性能方法

2023-12-21

我需要一种有效的方法来反转 python 中的 7000x7000 空气动力学影响系数(密集)矩阵。在使用 FORTRAN 例程之前,我已经开始使用 LAPACK 中的 LU 分解例程来处理问题,我已经看到它在其他相关应用程序中的使用非常有效。不过,我读到,NumPy 和 SciPy 线性系统求解器大多基于直接调用 C 中相同的 LAPACK/BLAS 函数,并且想知道切换到 FORTRAN 是否真的会将计算时间减少到足以证明放弃的水平一种更简单、更高级的语言。

如果有 Python 求解器可以保证该大小(1000 到 10000,平方)的矩阵具有相似的性能,那么它们是哪些?

我确实需要矩阵逆,所以切换到迭代 Ax=b 解决方案不是一个选择。


事实上,Numpy 和 Scipy 有效地调用 LAPACK 例程来执行numpy.linalg.inv https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html and scipy.linalg.inv https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html.

要逆一般矩阵,numpy.linalg.inv https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html solves A.x=np.eye((n,n))。功能inv() https://github.com/numpy/numpy/blob/v1.17.0/numpy/linalg/linalg.py#L486-L552 calls ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj), which calls https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src#L1693 call_@lapack_func@(&params); where params.B是单位矩阵并且@lapack_func@是其中之一sgesv, dgesv, cgesv, zgesv,它们是一般矩阵的线性求解器。

另一方面,scipy.linalg.inv calls https://github.com/scipy/scipy/blob/v1.4.1/scipy/linalg/basic.py#L911-L983 getri, 定义为get_lapack_funcs(('getri'),(a1,))。它对应于DGETRI() http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f_source.htmllapack 函数,旨在使用 LU 分解计算矩阵的逆,计算公式为DGETRF(). 因此,如果您正在使用DGETRI()在 Fortran 中,利用scipy.linalg.inv()在 python 中可能会实现类似的性能和结果。

大多数 Lapack 函数可以使用以下方式调用scipy.linalg.lapack https://docs.scipy.org/doc/scipy/reference/linalg.lapack.html#module-scipy.linalg.lapack。这是一个使用的示例scipy.linalg.cython_lapack.dgetri()在 cython 模块中:如何为Python编译C扩展,其中C函数使用LAPACK库? https://stackoverflow.com/questions/52106978/how-to-compile-c-extension-for-python-where-c-function-uses-lapack-library/52132359#52132359下面是一个示例代码,在 1000x1000 矩阵上比较 scipy.linalg.cython_lapack.dgetrf()+scipy.linalg.cython_lapack.dgetri() 、 numpy 和 scipy.linalg.inv() :

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

中提供了链接 DGETRF()+DGETRI() 的 Fortran 代码LAPACK 反演例程奇怪地混合了所有变量 https://stackoverflow.com/questions/26475987/lapack-inversion-routine-strangely-mixes-up-all-variables进行一些更改后,运行:

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()如图所示将 ATLAS/MKL 链接到已安装的 Numpy https://stackoverflow.com/questions/21671040/link-atlas-mkl-to-an-installed-numpy或报告中的命令如何在 NumPy 和 SciPy 中检查 BLAS/LAPACK 链接? https://stackoverflow.com/questions/9000164/how-to-check-blas-lapack-linkage-in-numpy-and-scipy .

为了进一步利用分布式计算,petsc https://www.mcs.anl.gov/petsc/documentation/faq.html#invertmatrix不鼓励对完整矩阵求逆,因为很少需要这样做。书中还写道MatMatSolve(A,B,X), where B and X可以使用稠密矩阵来做到这一点。此外,这个函数是在python接口中提供的petsc4py https://www.mcs.anl.gov/petsc/petsc4py-current/docs/apiref/index.html作为方法matSolve(self, Mat B, Mat X)对于对象petsc4py.PETSc.Mat。不再维护的元素库 https://github.com/elemental/Elemental被列为实现密集矩阵的直接求解器。虽然 Elemental 库支持 python 接口,但它的分支 Hydrogen 不再支持它。 尽管如此,Elemental 页面还是列出了一些与分布式密集线性代数相关的开源项目。ScaLapack 提供了例程PDGETRI()/PZGETRI() http://www.netlib.org/scalapack/explore-html/d3/df3/pdgetri_8f_source.html使用 LU 分解来反转分布式稠密矩阵。这可能会为更快的反转留下一些空间。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

使用 python 求解 7000x7000 线性系统时的最佳性能方法 的相关文章

随机推荐

  • Excel Range 奇怪行为的 SpecialCells 方法

    我编写了一个宏 使用 Excel 范围对象的 SpecialCells 方法从某个范围中查找空白单元格 当我尝试执行以下代码时 出现 未找到单元格 的异常 Sub test Debug Print Sheet1 Range A1 D4 Sp
  • 无法找到组合@Published - Xcode11 Beta 5(11M382q)

    我正在尝试使用以下内容运行一个简单的项目 Published var currentPlacemark CLPlacemark nil XCode 11 Beta 5 11M382q iOS13 17A5556d 出现以下错误 dyld S
  • azure辅助角色中的异步/等待导致角色回收

    我正在我的 WorkerRole RoleEntryPoint 中使用任务 异步和等待 我有一些无法解释的回收 现在我发现 如果等待调用中的某些内容运行时间过长 则角色会回收 要重现它 只需在 Run 方法中执行 await Task De
  • Jquery UI 可调整大小 - 调整放置在 iframe 上的 div 的大小

    如果你查看这个 jsbin http jsbin com efosed 5 edit http jsbin com efosed 5 edit然后你按 Run with JS 就会出现一个可以用 jquery ui 调整大小的 div 一切
  • Azure Functions 部署时无法运行

    我是新来的 如果帖子不完整 抱歉 我正在尝试在 azure 上部署一个与 blob 交互的 python 脚本 该脚本在本地运行良好 我可以与我的存储帐户交互 上传和下载 blob 但是当我在 azure 上部署我的函数时 它不会运行 日志
  • bash 使用序列号批量重命名文件夹和子文件夹中的文件

    我需要一个 bash 脚本来执行以下操作 对于文件夹及其子文件夹中存在的特定类型的每个文件 它都会在前面添加一个序列号 4 位数字 后跟一个分隔符 例如我有 Queen 1986 A Kind of Magic 01 One vision
  • 流程图 - 动态更改 y 轴

    我是飞行新手 但很快就设置了我的时间图 这是我基于时间的情节 plot placeholder d xaxis mode time minTickSize 1 month min new Date 2008 05 20 getTime ma
  • 为什么我需要将“get”包装在 J“lapply”调用中的虚拟函数中?

    我希望通过类或常见模式匹配等标准来处理列grep 我的第一次尝试没有成功 require data table test table lt data table a 1 10 ab 1 10 b 101 110 this does not
  • 在 Netbeans 内运行时停止 Tomcat

    我使用 NetBeans 运行 Apache Tomcat 6 当我的代码出现故障 例如 NullPointerException 时 tomcat 会失败并且不会运行任何其他请求 我的问题是我无法让 tomcat 停止 我必须重新启动整个
  • 查找 Java 应用程序中的连接泄漏

    我有一个应用程序在一段时间后开始出现内部服务器错误 我询问的一些人告诉我 这可能是因为我的应用程序中的连接泄漏 我开始搜索并发现这个查询来模拟连接泄漏 select LAST CALL ET SQL TEXT username machin
  • 堆积条形图未正确更新 d3js

    In this https plnkr co edit X7JYRLCKgBnasP86FRgQ p preview堆积条形图我添加了一个平分线和一个自定义x invert函数 以便您可以读取每个月的值 问题是 当我添加此自定义函数时 团队
  • OpenXML SDK 2.0 与 Aspose 在 .NET 中生成服务器端 Word 2007 文档

    我将在 Net 中启动一个服务器端办公自动化项目 以下是计划的主要活动 创建一个word文档 使用现有的包含封面 页眉 页脚 目录的 Word 文档模板 保存存档 嵌入文件并调整大小 HTML 图像 Word Excel TOC 生成和格式
  • 我无法从数据库 PostgreSQL 生成 Hibernate 映射文件和 POJO?

    已经在数据库 PostgreSQL 中创建了表和关系 但是当我想生成 Hibernate 映射文件和 POJO 时 它们没有生成 我应用了所有适当的步骤hibernate cfg xml一代和hibernate reveng xml 我认为
  • 如何在没有数据库的情况下配置 Ruby on Rails?

    对于当前不需要数据库的小型网站项目来说 使用 Ruby on Rails 会很方便 我知道我可以在 MySQL 中创建一个空数据库并从那里开始 但是有人知道在没有数据库的情况下运行 Rails 的更好方法吗 Thanks For Rails
  • 对于矩阵向量乘法,行优先排序是否更有效?

    If M是一个 n x m 矩阵并且v and u是向量 那么就索引而言 矩阵向量乘法看起来像u i sum M i j v j 1 lt j lt m Since v是一个向量 对于面向数值计算的语言 其元素可能存储在连续的内存位置中 如
  • python 在pdf文件中搜索

    这是pdf结构的一部分 5 0 obj lt lt Length 56 gt gt stream BT F1 12 Tf 100 700 Td 15 TL JavaScript example Tj ET endstream endobj
  • matplotlib中如何限制y轴高度?

    如何限制matplotlib图中y轴的高度 我正在尝试显示 x 轴 并降低该一维图的图形高度 我尝试过设置刻度 图形大小 tight layout 边距等 但没有成功 另外 无论我选择什么限制 更改 ylimit 都会跨越整个图形高度 im
  • Maven 故障安全插件不运行并行测试

    我有一个 Maven POM 文件 当我提供并行执行选项时 我在日志中没有看到任何并行执行的迹象 XML 调试让我抓狂 有什么想法这里出了什么问题吗
  • 查找两条曲线之间的重叠面积

    我一直在努力寻找解决方案来找到两条曲线之间的重叠区域 我处理的不是具有已知参数的概率密度函数 而是通过平滑经验数据点获得的曲线 我发现的唯一提示是计算不重叠的区域 如这段代码 来自here https www researchgate ne
  • 使用 python 求解 7000x7000 线性系统时的最佳性能方法

    我需要一种有效的方法来反转 python 中的 7000x7000 空气动力学影响系数 密集 矩阵 在使用 FORTRAN 例程之前 我已经开始使用 LAPACK 中的 LU 分解例程来处理问题 我已经看到它在其他相关应用程序中的使用非常有