Block matrices and Schur complement preconditioning

[Jed Brown]

For context, I'm thinking of LNKS optimization or coupled indefinite
systems, perhaps with many fields, where the Jacobian is applied matrix-free
(using a MatShell in my usage).  In these circumstances, it isn't useful to
actually form all the blocks in the Jacobian when constructing the
preconditioning matrix, just those which are the operator or preconditioning
matrix in a KSP inside the Schur complement/reduced space preconditioner.
In order to make the preconditioner slightly more generic, it would be nice
to have a matrix type which is really just a wrapper for the blocks which
are actually needed by the preconditioner.  For a simple concrete example,
consider a mixed discretization of the Stokes problem  J = [A B'; B 0]  where
J is applied matrix-free.  For preconditioning, we'll need an approximate
Schur complement  S = -B \tilde{A^{-1}} B'  where  \tilde{A^{-1}}  may be
V-cycle of AMG applied to  \hat{A}  (an approximation to  A) which is the
only matrix which needs to be actually formed.  Normally the pressure mass
matrix  M  would also be formed to precondition  S.  Now we could define the
preconditioning matrix  P = [\hat{A} 0; 0 M], but I don't like it.

Of course  \hat{A}  and  M  need not have the same matrix type, but it seems
logical to assemble them in the Jacobian assembly function.  What I
currently do is to put them in the PCShell context and assemble them in
PCSetUp(), but this means that the PC needs access to the problem
description.  We could put them in P = [\hat{A} 0; 0 M] (of type MatShell)
and extract the pieces with MatGetSubMatrix(), but it seems to me that
MatGetSubMatrix() ought to be able to succeed with any valid pair of IS.

So perhaps it would be useful to have another function MatGetSubBlock()
which just takes a pair of integers and returns the block if it is
available.  That is, we could have MatGetSubBlock(P,0,0,&Ahat) give the
explicit viscosity matrix and MatGetSubBlock(J,1,0,&B) give a MatShell which
implements B, but MatGetSubBlock(P,1,0,&Bhat) return NULL since an explicit
form of that matrix is not available (or it could return the MatShell B).
On the other hand, what should MatGetSubBlock(P,1,1,&C) give?  Logically, it
should be the zero matrix or NULL, but the preconditioner needs to get
access to M somehow.  So it looks like we are back to the original situation
where the assembly code needs to know details about the preconditioner or
the preconditioner needs to know how to assemble the matrices it needs.  We
can get slightly more separation by using P as a container for those
matrices that the preconditioner would need, but now its looking like P
isn't a matrix at all (it wouldn't really make sense to implement
MatMult()).

Any ideas for a better way?

Jed
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 197 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20080811/2ca2a2e0/attachment.sig>