[petsc-dev] Preconditioning systems with a submatrix of a larger matrix

[Patrick Farrell]

Hi Matt, Jason, and everyone else,

I have a big system (for two variables, u and alpha) and I'd like to solve this
by "alternate minimisation": block Gauss-Seidel-Newton, first solving for u and
then for alpha ("smaller systems"). I'd like to write the code in a generic way,
like a nonlinear fieldsplit, and in an ideal world contribute this back
upstream. It would be very useful for lots of things (like PDE-constrained
optimisation, where you cycle through forward - adjoint - gradient a few times
before starting Newton).

I have the first draft (in Python) written and working; it takes in a list of
index sets, creates subsneses for each subproblem, and wires them up with
residual and Jacobian functions that take the relevant subsets.  I've only
implemented the multiplicative version so far, although implementing the
additive variant would be straightforward.

I'm stuck on one point: I'd like to do this matrix-free, where the Jacobian of
the big problem is a shell matrix (for deflation, but also for time-dependent
PDE-constrained optimisation). It's straightforward (or will be when Jason's
commits land) to set up the sub-Jacobians with MatCreateSubMatrix. But I'm not
sure how to precondition the small Jacobian arising in applying Newton to each
subproblem. Given a preconditioner object for the big coupled system, is there a
way to "notify" the PC to say: here is your operator, which is the big Jacobian
with an index set attached? In my case with deflation, if the PC object could
know the right index set, I can construct an effective preconditioner for the
problem; I just don't know how to communicate that information from the solver
to the PC.

The same question arises in matrix-free preconditioning of variational
inequalities with active set methods, so I thought you might have thought about
this before.

Thanks for your advice!

Patrick