[petsc-dev] LSC question

[Barry Smith]

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.