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 <[email protected]> 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