Block matrices and Schur complement preconditioning

[Matthew Knepley]

Just so I can remove this from my todo mail pile, it appears Barry is
doing this.
The actually Schur stuff is a new KSP thing so we do not link Mat to KSP, but
it should work as you want.

   Matt

On Mon, Aug 11, 2008 at 4:16 PM, Jed Brown <jed at 59a2.org> wrote:
> 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
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which
their experiments lead.
-- Norbert Wiener