如何在scipy中实现ILU预处理器？

Question

对于scipy.sparse.linalg中的迭代解算器，例如bicg、gmres等，有一个选项可以为矩阵A添加预处理器。然而，documentation对于我应该给出什么作为预条件并不十分清楚。如果我使用ilu = sp.sparse.linalg.spilu(A)，ilu不是任何矩阵，而是包含许多事物的对象。你知道吗

有人问了一个类似的问题here关于python2.7，但我不适合我（python3.7，scipyversion1.1.0）

所以我的问题是如何将不完全LU预条件加入到这些迭代算法中？你知道吗

Answer 1

作为预条件，bicg或gmres接受

• 稀疏矩阵
• 密集矩阵
• 线性算子

在您的例子中，预条件来自一个因式分解，因此它必须作为一个线性运算符传递。你知道吗

因此，您可能希望从通过spilu获得的ILU分解中显式定义一个线性运算符。一些类似的东西：

sA_iLU = sparse.linalg.spilu(sA)
M = sparse.linalg.LinearOperator((nrows,ncols), sA_iLU.solve)

这里，sA是CSC格式的稀疏矩阵，M现在将是提供给迭代解算器的预条件线性运算符。你知道吗

---

一个完整的例子based on the question you mentioned：

import numpy as np
from scipy import sparse
from scipy.sparse import linalg

A = np.array([[ 0.4445,  0.4444, -0.2222],
              [ 0.4444,  0.4445, -0.2222],
              [-0.2222, -0.2222,  0.1112]])
sA = sparse.csc_matrix(A)

b = np.array([[ 0.6667], 
              [ 0.6667], 
              [-0.3332]])

sA_iLU = sparse.linalg.spilu(sA)
M = sparse.linalg.LinearOperator((3,3), sA_iLU.solve)

x = sparse.linalg.gmres(A,b,M=M)

print(x)

注意事项：

• 我实际上是以密集矩阵为例，而在你的例子中，从一个有代表性的稀疏矩阵开始会更有意义。你知道吗
• 线性运算符M的大小是硬编码的。你知道吗
• ILU尚未配置，但使用默认值。你知道吗
• 这几乎就是人们所建议的in the comments to that aforementioned answer，然而，我不得不做一些小的调整以使它与Python3兼容。你知道吗