[petsc-dev] LSC question

[Jed Brown]

On Fri, Mar 9, 2012 at 09:29, Barry Smith <bsmith at mcs.anl.gov> wrote:

>      As you know, currently when using PCFIELDSPLIT  we have the
> presentation:
>
>      ilink = jac->head;
>      ierr  = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
>      ierr  =
> MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr);
>      ierr  = ISDestroy(&ccis);CHKERRQ(ierr);
>      ilink = ilink->next;
>      ierr  = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr);
>      ierr  =
> MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
>      ierr  = ISDestroy(&ccis);CHKERRQ(ierr);
>      /* Use mat[0] (diagonal block of the real matrix) preconditioned by
> pmat[0] */
>      ierr  =
> MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr);
>
>      I think it makes sense for the L to be attached by the user to the
> user provided pc->mat and pc->pmat , (for example the ComputeFunction() the
> user provides could attach the L directly to the Jacobian they produce)
> with our current design this means PCFIELDSPLIT here needs to pull the L
> out of the pc->mat and put it into the jac->schur. Kind of ugly but
> workable.
>
>      An alternative is to refactor MatCreateSchurComplement() to take the
> original matrix and the IS as arguments
>

That routine is named MatGetSchurComplement() and the user can implement
it. That is the right approach. It could have a convention for forwarding L
and Lp on, but the user could just call MatGetSchurComplement() and compose
"L"/"Lp" with the Schur complement. This is better because there are
multiple ways to take a Schur complement, thus "L" cannot be defined
without knowing the index sets.

Additionally, the original matrix does not currently retain ownership, so
the user can't actually attach anything to the Schur complement without
implementing MatGetSchurComplement() themselves.

Worse yet, MatGetSchurComplement() is probably going to call
MatGetSubMatrix() with exactly the same arguments as the outer
PCFieldSplit, because both of them will need the submatrices.


> and basically suck in all the code above into the
> MatCreateSchurComplement() this routine would automatically pull the L out
> of the input matrix and stick it into the created Schur matrix,
>

This does not work.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120309/c9e09fe6/attachment.html>