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

[Matthew Knepley]

On Thu, Apr 9, 2015 at 3:17 AM, Patrick Farrell <
patrick.farrell at maths.ox.ac.uk> wrote:

> 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.
>

So your code is more like Fieldsplit than the SNESMultiblock?


> 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.
>

Just so I understand, you are talking about a method that projects the
problem onto
a subspace, so

  F(u) = 0     ==>     P^T F (u) = 0

which I assume has Jacobian

  P^T F' P

and you would do something like

  P^T \hat J P

so you want to know P when calculating a MF PC. Is this right? I can see two
related ways to do this. First, the information about P should be available
from
the solver, which you get passed in the PC construction. Second, we could
have
a mode for these projected solvers where the Jacobian PC is also projected.

  Thanks,

     Matt


> Thanks for your advice!
>
> Patrick
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20150409/1f2289c7/attachment.html>