Interfaces for Schur complements and FieldSplit

[Jed Brown]

Barry Smith wrote:
> 
> On May 20, 2009, at 3:47 PM, Jed Brown wrote:
> 
>> Do you have a preferred interface for specifying a preconditioning
>> matrix for A in Mat_SchurComplement
> 
>    I would add a new second argument to MatCreateSchurComplement() Ap
> that would be the pmat passed into the KSPSetOperators(). To get the
> current behavior one would use MatCreateSchurComplement(A,A,....)
> 
>> and for the Schur complement in PC_FieldSplit?
> 
> PCFieldSplitSchurPrecondition() could be modified to take an
> additional argument that is used as the preconditioner matrix.

This is what I had in mind so I'll do it, but maybe not right away.

>> It looks like I must implement MatGetSubMatrix for my matrix-free
>> operators.  This is of course only possible for very special IS.  Is
>> there any chance of a higher-level description of the submatrix?  (I'm
>> not sure ther is a better way, but as is, both PC_FieldSplit and my
>> Mat_Shell already need to have duplicate representations of the
>> submatrices but we'll be talking by complementing each IS.)
>>
>   Not likely. The interface is to provide an IS. It will always be an
> IS, but of course IS is an abstract base class and perhaps additional
> ISs could be introduced that are "higher level" but in the end they
> still have to be an abstract representation of a bunch of integers :-(

That's fine, and it's only done once so even the naive representation
isn't a performance issue.

>> Do you have an interface in mind for specifying reduced forms for block
>> preconditioners?  For example, Galerkin discretization of incompressible
>> flow gives
>>
>> J = [A  B']
>>    [B  0 ]
>>
>> The reduced forms
>>
>> [A    ]   or    [A  B']
>> [B   S]         [   S ]
>>
>> are appropriate for left- and right-preconditioned GMRES respectively
>> (the preconditioned operator has minimal degree 2 and all unit
>> eigenvalues).  Also,
>>
>> [A    ]
>> [    S]
>>
>> is popular for symmetric problems.
>>
>>
>    Not sure what you mean "specifying". See the manual page for
> Block_Preconditioners (bottom of fieldsplit.c) where it gives the
> cases. Sadly it does not tell you directly how to run each different
> case (It should) but I think you just need to use the proper
> -pc_fieldsplit_type additive,multiplicitive,symmetric_multiplicative

I've looked at that.  This is different than what is described there
because all of my examples are variants of a Schur-complement
preconditioner where as those examples are block relaxation of some
sort.  As I see it, there are two sources of preconditioners for these
systems, relaxation and factorization.  So additive, multiplicative, and
symmetric_multiplicative are all block relaxation where as Schur is the
canonical example of factorization.  What you have set up is nice and
flexible for relaxation, but I'm having trouble seeing how to extend the
factorization choices.

Some more thoughts on this:

So factorization always produces Schur complements and it is rare to get
a scalable Schur solver without finding a way to interpret the Schur
complement physically, or at least use some approximate commutators, to
build a preconditioner.  Since factorization of larger systems produces
very complex Schur complements they will be very hard to understand, so
I think it is useful to do both relaxation and factorization.  There is
a good chance that the right way to do this is to nest a factorization
preconditioner as one block inside a larger split based on relaxation.

Unfortunately, factorization preconditioners don't necessarily require
preconditioners for diagonal blocks.  For instance, consider the
Jacobian in LNKS.  This is a 3x3 (state, adjoint, design) system which
is indefinite about the adjoint.  The factorization you end up with
needs preconditioners for the forward and adjoint models which sit at
the (adjoint,state) and (state,adjoint) blocks.  This case is not
reducible to a 2x2 form and the factorization requires pivoting.  The
reduced Hessian (the block that ends up on the diagonal in the design
space) is of course a nasty thing and you might not have a
preconditioner available for it.

Within the factorization variants, it is very common that you want to
drop either the upper or lower factor.  That is, suppose there is no
pivoting and you have factors J = L D U.  If the system is nonsymmetric,
and sometimes if it is, you usually want to solve with (LD) or (DU), but
not with the whole (LDU).  This is because the effect of the dropping
the factor is well-controlled by GMRES (usually the preconditioned
operator has a low-degree minimal polynomial) and you need one more
inner preconditioner application if you use the full factorization.

I think it would be very useful to come up with an interface for this so
that the user only needs to specify a factorization and provide
preconditioners for the resulting Schur complements.  Maybe such
problems never really get bigger than 3x3 or 4x4 so they can be
special-cased, but I think it ought to be possible to support any size
and arbitrary factorizations.


A related issue is telling the application (UserJacobian) what auxiliary
systems are going to be needed by the preconditioner.  For example, if
you switch from factorization to relaxation, the application no longer
needs to assemble an auxiliary system used to precondition the Schur
complement.  I can think of conventions like attaching a Mat with a
special name to the (global, packed) preconditioning matrix (can be
wrapped in an API call), but that would mean that the PC needs to have
seen the matrix *before* SNESComputeFunction is called, and it just
seems dirty and backwards.  Probably the right thing is to have a
PCFieldSplit API, so UserJacobian would unwrap the PC from SNES and
query for which pieces are needed.  (My motivation here is that we
definitely don't want UserJacobian snooping -pc_fieldsplit_* to find out
what will be needed.)  For more general factorization, just identifying
the blocks of the factorization seems tricky (at least coupled to the
way you specify the factorization).


Thanks for the discussion.

Jed

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 260 bytes
Desc: OpenPGP digital signature
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090521/9264e44d/attachment.sig>