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?
Python environments for AMG and Gauss Seidel as solvers instead of preconditioners
linear algebra;python;preconditioning;multigrid
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:Get the matrix out of FEniCS and into scipyDo the linear algebra in scipy, using my own code or packages like pyamgBring 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 npimport scipy.sparse as spsimport scipy.sparse.linalg as splafrom 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 Mdef 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 xmesh = 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*dxA = 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 mlx_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.