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?


2 Answers 2


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))
    return M

def do_gauss_seidel(b, A, tol=1e-5, maxiter=1000):
    dd = A.diagonal()
    D = sps.dia_matrix(A.shape)
    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:
        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.

  • $\begingroup$ Great insight into FEniCs Nick! I took a long time to realise this myself about it. So it seems that some of the methods like Gauss Seidel I have to write down own code ( thanks for your version ;) ). For AMG it is indeed pyAMG. I need to wade through its documentation a little to see which routines should help. Your pointer was definitely helpful. Thanks! $\endgroup$
    – Aseem Dua
    Commented Jul 26, 2016 at 10:18
  • 1
    $\begingroup$ I like pyAMG's root_node_solver; it is based on some very advanced modern techniques, and I've gotten good results from it on difficult 4th order problems that other AMG methods struggle on. Also, a lot of good documentation for pyAMG exists within the code itself - the comments at the beginnings of functions are very good and explain a lot about how to use the function, etc. $\endgroup$
    – Nick Alger
    Commented Jul 26, 2016 at 12:58

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.

  • $\begingroup$ Doesn't PETSc also just provide for amg and gauss seidel as preconditioners instead of standalone solvers? $\endgroup$
    – Aseem Dua
    Commented Jul 25, 2016 at 21:45
  • $\begingroup$ Yes, and you can use petsc4py for those. But adding that huge library just for the AMG can be overkill. I gave you an extra little choice in case you needed minimal dependencies. However, PETSc (via petsc4py) is definitely the way to go if you need all the power. $\endgroup$ Commented Jul 25, 2016 at 21:48
  • $\begingroup$ I am actually looking for amg and gauss seidel type methods as standalone solvers. AMG as a preconditioner is there in native scipy itself. $\endgroup$
    – Aseem Dua
    Commented Jul 26, 2016 at 6:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.