Interfaces for Schur complements and FieldSplit

[Jed Brown]

Jed Brown wrote:

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

I thought about this a bit more and have a proposed interface for field
split preconditioners based on factorization.  Let's give it some time
and see if it has any fatal flaws.  I'll be visiting Argonne on June 30
and July 1 so perhaps we can discuss it more then.

My proposal is that the user registers "symbolic factorization" which
would include pivoting, fill for the incomplete factors, and perhaps
names for the resulting terms for which a preconditioner would be
needed.  These symbolic factorization strategies would be indexed by
name (FList), and some common cases could be provided.  So the runtime
interface would look like

 -pc_fieldsplit_factor_strategy lnks_p2

Note that the location of the complement is not a complete description
without the pivoting strategy, and would not be robust to changing the
number of fields in the split or changing the fill, which is my
motivation for using names.

PC_FieldSplit would perform "numeric factorization" using the strategy
specified (e.g. named on the command line).  UserJacobian would call

  PCFieldSplitFactorGetSubKSP(PC,const char *name,KSP*);

a possible alternative is

  PCFieldSplitFactorGetSubKSP(PC,const char *strategy,PetscInt n,KSP*);

where the numbering convention needs to be explained in the
documentation for each strategy.  If this returns non-null, it means
that the factorization strategy in use needs a preconditioner for the
named operator.  PC_FieldSplit_Factor would be providing the operator
(via Mat_SchurComplement, or just from MatGetSubMatrix) so UserJacobian
would not need to change that, but would be providing a preconditioning
matrix based on the PC (since the preconditioning matrix is usually
MatShell or something exotic).  For example,

  S = -B Au^{-1} B^T

might be preconditioned with

  P = C Ap^{-1} Mp

where C is a rediscretization of the elliptic operator (-B Mu^{-1} B^T),
Ap is a translation of Au to the pressure space (i.e. discretizing a
scalar analogue of the vector equations discretized by Au), and Mu,Mp
are (possibly lumped) velocity and pressure mass matrices.  This
preconditioner would have a name (here PCKLW for Kay-Loghin-Wathen) and
requires a special matrix type (basically Mat_Composite).  Of course

  P^{-1} = Mp^{-1} Ap C^{-1}

is what we care about.  In my case, UserJacobian would be providing
matrix-free implementations of Mp, Ap, and C (because these require zero
setup, typically only Ap would be used with the others silent because of
preonly), assembling a lumped mass matrix Mp, and assembling a
preconditioning matrix for C.  We would then wrap this all up in a
suitable matrix type (understood by PCKLW) as the preconditioning
matrix, or call something like PCKLWSetOperators().

PCBFBT (or the equivalent but better-scaled least squares commutator,
call it PCLSC) could be put under this framework.  The only unique thing
about it is that it uses a different representation for the Schur
complement, namely

  P = C (B Mu^{-1} Au Mu^{-1} B^T)^{-1} C

where C has the same meaning as above, a discretization of the elliptic
operator

  -B Mu^{-1} B^T

(perhaps computed algebraically).


In light of these examples, do you have a preference between using
matrix types that can carry this information (the pieces used by PCKLW
and PCLSC) and calling the PC interface directly, as in
PC{KLW,LSC}SetOperators?




So I'm thinking of a new type with definition something like

struct _p_PCFieldSplitFactorStrategy {
  PETSCHEADER(int);
  PetscInt splitsize;            /* number of fields in split       */
  PetscInt *leftperm,*rightperm; /* P and Q  in  P L D U Q          */
  PetscInt *fill;                /* mask nonzero locations in L+D+U */
  /* perhaps also names for terms needing preconditioning */
};

interface:
  PCFieldSplitFactorStrategyRegister(const char[],const char[],const char[],PetscErrorCode(*)(PCFieldSplitFactorStrategy));

either a convention is used for numbering the factors:
  PCFieldSplitFactorGetSubKSP(PC,const char*,PetscInt,KSP*);
or they are explicitly named:
  PCFieldSplitFactorGetSubKSP(PC,const char*,KSP*); /* 


This gives runtime control in the form

 -fieldsplit_factor_lnks_p2_8_{ksp,pc}_*

If instead the resulting factors had names, we could have

 -fieldsplit_factor_reduced_hessian_{ksp,pc}_*

or the more verbose

 -fieldsplit_factor_lnks_p2_reduced_hessian_{ksp,pc}_*



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/288939b2/attachment.sig>