PC_FieldSplitLink

[Barry Smith]

   I just used a link list to support flexibly adding more fields  
without needing to allocate a new array and copying over the old values.
If array storage makes the implementation easier that is fine with  
me.  Since there are not many blocks why bother with CSR to store the
(say 5 by 5 = 25 values) I would just use a dense array.

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

    Yes, I've been wanting that for a while.

    Barry

On Jul 10, 2009, at 6:32 AM, Jed Brown wrote:

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