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
