PC_FieldSplitLink

[Jed Brown]

What is the use case where the linked list, as opposed to array storage is
important?

I want to extend the factorization part to support an arbitrary number
of fields with pivoting and possibly dropping coupling.  Suppose we have
a matrix with block form:

  A_00  A_01  A_02  A_03
  A_10  A_11  A_12  A_13
  A_20  A_21  A_22  A_23
  A_30  A_31  A_32  A_33

Introduce the notation

  S^{i,j}_{k,l} = A_{k,l} - A_{k,j} inv(A_{i,j}) A_{i,l}

which is the Schur complement obtained by eliminating the ij block in kl.
The matrix A factors as LU where

L = 

  1
  A_10 inv(A_00)    1
  A_20 inv(A_00)    S^00_21 inv(S^00_11)     1
  A_30 inv(A_00)    S^00_31 inv(S^00_11)     S^{01,01}_32 inv(S^{01,01}_22)    1

U =

  A_00        A_01          A_02             A_03
              S^00_11       S^00_12          S^00_13
                            S^{01,01}_22     S^{01,01}_23
                                             S^{012,012}_33

Evidently we need to be able to apply several Schur complements and
solve with those on the diagonal of U.  In many cases, the user will
want to drop some or all of L and possibly off-diagonal elements in U.
As an implementation detail, consider

  S^{01,01}_22 = [A_20  A_21] inv([A_00  A_01]) [A_02]
                                 ([A_10  A_11]) [A_12]

The whole point of using factorization is that we don't have a good
preconditioner for the entire block

  [A_00  A_01]
  [A_10  A_11]

so this would be applied using factorization.  But, the user will likely
have used some physical argument to provide a preconditioner for
S^{01,01}_22 so we'll get that using MatGetSchurComplement().

Anyway, I'm having some trouble making this fit nicely into
PC_FieldSplitLink when pivoting is used.  I'm considering using a CSR format
to hold the factors that the user wants to have.  Note that with support for
factoring 3x3 systems with pivoting, we can make all the LNKS
preconditioning variants available as runtime options.


We could unify relaxation and factorization by providing the option to just
use the diagonal block instead of the Schur complement.  That is, choosing
submatrices instead of Schur complements and using a diagonal factorization
matrix is Jacobi.  The upper-triangular piece is Gauss-Seidel (equivalent to
the lower-triangular version after a permutation).  Suppose the user
specified the nonzero pattern of the factored matrix with the convention
that a 0 represents a dropped block, 1 is just a submatrix, and 2 uses the
Schur complement.  Then

  1 0 0 0
  0 1 0 0
  0 0 1 0
  0 0 0 1

is Jacobi (pure relaxation),

  1 1 1 1
  0 1 1 1
  0 0 1 1
  0 0 0 1

is Gauss-Seidel (also pure relaxation),

  2 2 2 2
  2 2 2 2
  2 2 2 2
  2 2 2 2

is the full factorization,

  2 2 2 2
  0 2 2 2
  0 0 2 2
  0 0 0 2

drops the lower factor (which is often desired for right-preconditioned
GMRES),

  2 2 1 1
  0 2 1 1
  0 0 1 1
  0 0 0 1

relaxes fields 2 and 3 multiplicatively, then solves for fields 0 and 1
using factorization.  We would expect this hybrid to work well when stepping
over a stiff wave coupling fields 0 and 1, but fields 2 and 3 interact on a
longer time scale.  To use this final variant efficiently, the user would
only need to provide a preconditioner for S^00_11.

Does this unification sound sensible?  It would be an overhaul of
PC_FieldSplit, although I think I can make it almost transparent to the user
and it should simplify some aspects.

A complication is that when the Jacobian is not assembled, it is easier to
apply an entire block-row of the Jacobian (Afield) than a single block.  The
current implementation of multiplicative relaxation is designed with this in
mind, but that means that the diagonal terms are needlessly applied.  It
seems to me that any matrix format that can provide the rows could also
provide the off-diagonal block (even if only by scattering the subvector to
a full vector and applying full rows).

As a related issue, do you think an implementation of MatGetSubMatrix_MFFD
would be a good thing?  It would simplify PC_FieldSplit, but might not be
worthwhile unless it would be useful elsewhere.

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/20090710/fc68b839/attachment.sig>