我正在研究块预处理，似乎为他们编写定制的 Krylov 求解器是很常见的。在每个求解器中，带有预处理器的单个块线性系统通过几次 AMG V 循环或 Gauss Seidel 扫描来近似。
然而，这些方法在 scipy 中用作平滑器，甚至 FEniC 支持 PETSc。我想让它们作为本机求解器，就像在 MATLAB 中可以做的那样。例如，只是为了独立运行一个 AMG 循环。但是，我在 python 环境（FEniCs + Scipy + Numpy）中工作，需要相同的指导。有没有我可以用来执行这些的库？
FEniCS 倾向于隐藏有关它构建的实际矩阵的细节，并防止对其进行轻松操作。据我所知，这是一个设计决策，因为他们正试图创建一个包罗万象的封闭管道来离散化和求解 PDES。矩阵存储和矩阵元素的特定细节被认为是只有开发人员才应该关注的实现细节。可用的线性代数后端格式因版本而异，破坏了旧代码。没有内置函数可以将稀疏矩阵导出为 scipy 格式（.toarray() 转换为密集矩阵！）；您必须指定后端类型，然后专门为它编写自己的导出函数。
对于只想指定 PDE 并获得解决方案的人来说，这可能很好，但对于像你我这样想要尝试构建不同的预条件子、求解器、研究矩阵的光谱特性等的人来说，这令人沮丧。所以，通常我做的是：
无论如何，我经常这样做，以至于我只花了几分钟就编写了以下执行此过程的简单代码，在一个简单的椭圆 PDE 上测试了一个简单的 Gauss-Seidel 求解器和来自 pyAMG 的多重网格求解器：
from dolfin import *
import numpy as np
import scipy.sparse as sps
import scipy.sparse.linalg as spla
from pyamg import *
dolfin.parameters.linear_algebra_backend = 'Eigen'
def fenics_eigen_to_scipy_csr(assembled_fenics_form):
row,col,val = as_backend_type(assembled_fenics_form).data()
M = sps.csr_matrix((val,col,row))
M.eliminate_zeros()
return M
def do_gauss_seidel(b, A, tol=1e-5, maxiter=1000):
dd = A.diagonal()
D = sps.dia_matrix(A.shape)
D.setdiag(dd)
L = sps.tril(A, -1)
U = sps.triu(A, 1)
solve_DL = spla.factorized(D + L)
normb = np.linalg.norm(b)
x = np.zeros(b.shape)
for ii in range(maxiter):
normr = np.linalg.norm(b - A*x)
if normr/normb < tol:
break
x = solve_DL(b - U*x)
return x
mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 1.5), 100, 200)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(nabla_grad(u), nabla_grad(v))*dx + 10*u*v*dx
A = fenics_eigen_to_scipy_csr(assemble(a))
b = np.random.rand(V.dim())
x_true = spla.spsolve(A, b)
x_gs = do_gauss_seidel(b, A, maxiter=1000)
print 'relative error in gauss-seidel solve:', np.linalg.norm(x_true - x_gs)/np.linalg.norm(x_true)
ml = ruge_stuben_solver(A)
print ml
x_amg = ml.solve(b, tol=1e-10, maxiter=2)
print "relative error in amg solve:", np.linalg.norm(x_true - x_amg)/np.linalg.norm(x_true)
# Get the vector back into a FEniCS Function and plot it:
xtrue_fct = Function(V)
xtrue_fct.vector()[:] = np.copy(x_true)
plot(xtrue_fct, interactive=True)
希望这可以帮助。
pyAMG是一个很好的代数多重网格求解器库。我不在 Python 中使用它，但我使用了Julia 包装的版本，并且发现它足够好，我很快就会在我的DifferentialEquations.jl中使用它。
对于其他人，您可能想尝试petsc4py。这些是 PETSc 的包装器。