Preconditioning for saddle point problems

Lisandro Dalcin dalcinl at
Tue Apr 29 08:44:24 CDT 2008

Well, I've worked hard on similar methods, but for incompressible NS
equations (pressure-convection preconditioners, Elman et al.). I
abandoned temporarily this research, but I was not able to get decent
results. However, for Stokes flow it seens to work endeed, but never
studied this seriously.

I'll comment you the degree of abstraction I could achieve. In my base
FEM code, I have a global [F, G; D C] matrix (I use stabilized
methods) built from standard linear elements and partitioned across
processors in a way inherited by the mesh partitioner (metis). So the
F, G, D, C entries are all 'interleaved' at each proc.

In order to extract the blocks as parallel matrices from the goblal
saddle-point parallel matrix, I used MatGetSubmatrix, for this I
needed to build two index set for momentum eqs and continuity eqs
local at each proc but in global numbering. Those index set are the
only input required (apart from the global matrix) to build the

On 4/29/08, Jed Brown <jed at> wrote:
> This issue arose while discussing improved preconditioners for a colleague's
>  libmesh-based ice flow model.  Since libmesh is one of several finite element
>  libraries that use PETSc as the backend, perhaps the abstraction should be
>  implemented in PETSc.
>  As a model, consider the Stokes problem with mixed discretization in
>  matrix form
>  (1)  [A, B1'; B2, 0] [u; p] = [f; 0]
>  where A is the viscous velocity system, B1' is a gradient and B2 is divergence.
>  The system (1) is indefinite so most preconditioners and many direct methods
>  will fail, however A is a uniformly elliptic operator so we can use exotic
>  preconditioners.  An effective preconditioner is to solve the Schur complement
>  problem
>  (2)  S p = -B2 A^{-1} f
>  where S = -B2 A^{-1} B1' is the Schur complement.  Once p is determined, we can
>  find u by solving
>  (3)  A u = f - B1' p
>  If we solve (2,3) exactly, this is the ultimate preconditioner and we converge
>  in one iteration.  It is important to note that applying S in (2) requires
>  solving a system with A.  This can be implemented in PETSc using a MatShell
>  object.  The classical Uzawa algorithm is a Richardson iteration on this
>  problem.  A better method (in my experience) is to precondition a Krylov
>  iteration (FGMRES) on (1) with an approximate solve of (2,3).  For instance I
>  have had success with 3 iterations of GMRES on (2) where the A^{-1} in S is
>  replaced by one V-cycle of algebraic multigrid and the A^{-1} appearing
>  explicitly in (2,3) are replaced by 4 iterations of GMRES preconditioned with
>  AMG.
>  It is desirable to have all the details of this iteration controllable from the
>  command line.  For an example implementation (a spectral method) of this
>  problem, check out StokesMatMultXXX() and StokesPCApply() at
>  Note that typically, S is preconditioned with the mass matrix (or a
>  ``pressure-laplacian'' at higher Reynolds number) but since this code uses a
>  collocation method and is only intended for slow flow, there is no
>  preconditioner for S.
>  My question for the PETSc developers: how much of this can/should be abstracted
>  out?  While it's not difficult to code, perhaps there should at least be an
>  example.  If it would be useful, I can strip my experiment above down to make it
>  suitable.
>  PS: An excellent review paper:
>  @article{benzi2005nss,
>   title={{Numerical solution of saddle point problems}},
>   author={Benzi, M. and Golub, G.H. and Liesen, J.},
>   journal={Acta Numerica},
>   volume={14},
>   pages={1--137},
>   year={2005},
>   publisher={Cambridge Univ Press}
>  }
>  Jed

Lisandro Dalcín
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594

More information about the petsc-dev mailing list