Python environments for AMG and Gauss Seidel as solvers instead of preconditioners

Question

I am working on block preconditioning and seemingly it is common to write customised Krylov solvers for them. Within each solver, the individual block linear system with preconditioners are approximated by a few runs of AMG V-cycle or a Gauss Seidel sweep.

These methods are however used as smoothers in scipy, or even FEniCs supported PETSc. I would like to have them as native solvers, like it is possible to do in MATLAB. Just to run a single AMG cycle standalone, for example. However, I am working in a python environment (FEniCs + Scipy + Numpy) and need guidance for the same. Are there any libraries that I can make use of to carry out these?

Answer 1

FEniCS tends to hide the details about the actual matrices it builds, and prevent easy manipulation of them. As far as I can tell this is a design decision, as they are trying to create an all-encompassing closed pipeline for discretizing and solving PDES. The particular details of matrix storage and matrix elements are considered an implementation detail that only developers should be concerned with. The available linear algebra backend formats change from version to version, breaking old code. There is no built in function to export a sparse matrix to the scipy format (.toarray() converts to a dense matrix!); you have to specify a backend type and then write your own exporting function specifically for it.

For someone who just wants to specify a PDE and get a solution, maybe this is fine, but for someone like you or me who wants to experiment with building different preconditioners, solvers, investigate spectral properties of the matrices, etc, this is frustrating. So, typically what I do is:

1. Get the matrix out of FEniCS and into scipy
2. Do the linear algebra in scipy, using my own code or packages like pyamg
3. Bring the results back into FEniCS for post-processing.

Anyways, I do this so often that it only took me a few minutes to whip up the following simple code that does this process, testing both a simple Gauss-Seidel solver and a multigrid solver from pyAMG on a simple elliptic PDE:

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)

Hope this helps.

Answer 2

pyAMG is a good library for algebraic multigrid solvers. I don't use it in Python, but I have used the Julia wrapped version and have found it good enough that I will be using it in my DifferentialEquations.jl pretty soon.

For others, you may want to try petsc4py. These are wrappers for PETSc.