正确使用scipy的sparse.linalg.spilu

计算科学 稀疏矩阵 scipy 预处理
2021-12-17 16:22:51

我正在尝试使用 scipy 的spilu例程作为预处理器,但我发现我的应用程序性能不佳(解决由时间相关 ADR PDE 的 DG 离散化产生的全局线性系统)。

在我编写测试以调整选项并计算预处理器对我的线性系统的特征值的影响之前,我想确保我以scipy预期的方式应用预处理矩阵,因为文档spilu非常简陋。 ..

elif solver == 'spILUPCG': P = sp.sparse.linalg.spilu(A.tocsc()) P = P.L * P.U x = sp.sparse.linalg.cg(A, b, M=P)

2个回答

如果给出的近似分解,你不会想使用作为预处理器(这大约是),你会想使用(大约是)。即使那样,您也不希望将它们显式地乘以,它可能很密集。而是将视为级联两个求解步骤的操作/运算符(首先由求解,然后由求解)。LUAP=LUAP=U1L1A1PPLU

我希望有某种方法可以将该操作作为闭包/回调注入 precondition cg()可悲的是,我不是 scipy 专家,所以我不知道确切的方法。但是在某个地方 scipy 应该有一个抽象来应用线性运算符的操作,我希望它们的矩阵对象和它们的预处理器对象都能实现这个抽象。我想我会尝试只传递你的对象本身。P

这个问题有一个例子,说明如何使用形状M的 scipy 稀疏矩阵创建预处理器ANxN

from scipy.sparse.linalg import LinearOperator, spilu
ilu = spilu(A)
Mx = lambda x: ilu.solve(x)
M = LinearOperator((N, N), Mx)