On Mar 9, 2012, at 11:10 AM, Jed Brown wrote:
> 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()
Sorry I missed that.
> 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.
Please modify MatGetSchurComplement() to forward the L and Lp. I think is
should be automatic.
> 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.
Yes, this is why we cannot use MatGetSchurComplement() inside the PCFIELDSPLIT
So we are stuck with having PCFIELDSPLIT forward the and L and Lp from the
outer matrix. Will you add this?
Barry
>
> 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.
