[petsc-dev] Promote MatShellSetOperation() and MatShellSetContext() to Mat

[Jed Brown]

On Mon, Jun 6, 2011 at 21:21, Barry Smith <bsmith at mcs.anl.gov> wrote:

> Yes, that is the whole point (in the same way I want DM's propagated down
> inside mulgrid and fieldsplit created subksps and pcs) So for example it can
> be set in the main program into the SNES (or TS) and is accessible deep
> inside some nested preconditioner that is a problem specific PC beasty
> thing.
>

So with PCFieldSplit, a DMComposite should be taken apart and have just the
pieces forwarded on. The global DMComposite wouldn't be much use for a
nested PC.

The same occurs with the application context. I imagine that each branch of
the PCFieldSplit is almost an independent code that implements the right
interface. It should not need to take apart the global app context to find
itself. So here's some top-level code. The vht object is the top level app
context.

      PCFieldSplitSetIS(pc,"s",vht->all.sblock);
      PCFieldSplitSetIS(pc,"e",vht->all.eblock);
      PCFieldSplitGetSubKSP(pc,&nsub,&subksp);
      stk_ksp = subksp[0];
      KSPGetPC(stk_ksp,&stk_pc);
      PCSetType(stk_pc,PCFIELDSPLIT);
      PCFieldSplitSetIS(stk_pc,"u",vht->stokes->ublock);
      PCFieldSplitSetIS(stk_pc,"p",vht->stokes->pblock);
      KSPSetApplicationContext(stk_ksp,vht);

It would be better to replace the last five lines with

      StokesKSPConfigure(vht->stokes,stk_ksp);

Which would end up using KSPSetApplicationContext() to set the stokes
context on the nested ksp. Note that the Stokes context no longer needs to
find itself within a larger scope, and it no longer even has access to that
higher level.


> > What happens when a component is accessed from multiple places?
>
>    Not sure what you mean by multiple places.
>

For example, a PC containing expensive setup (e.g. direct solve) might be
used as part of a stand-alone linear solve and as part of a fieldsplit
preconditioner for some larger nonlinear problem. If the user had an
application context on that thing, it probably shouldn't be overwritten by
SNESSetApplicationContext().

The argument here is just that an inner object might not have a unique
parent. In such cases, it would seem easier to manage if the first parent
(or a direct call by the user) took precedence instead of the last parent.


> >
> > I wonder if this forwarding should be made to not overwrite what is
> already there. That is, if the user sets it explicitly, then it will not be
> overwritten by a parent.
>
>    Hmm, yes this gets tricky. I general I think it should overwrite and the
> user needs to set it differently if they want after they set it on the
> outside.  So they do TSSetApplicationContext() and then TSGetSNES,
> SNESGetKSP, KSPGetPC etc etc to change it down there.
>

This is fine, but adds a little order dependence at configuration time. It
does not seem too onerous if the user does something like
StokesKSPConfigure() above.


> > How should MatShellSetOperation behave? Should it error or do nothing
> when the matrix type is not MATSHELL?
>
>    Error
>
>
Okay, this is not the current behavior and this pre-historic snippet
indicates the current behavior is intentional.

  if (op == MATOP_DESTROY) {
    ierr = PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
    if (flg) {

> If the current behavior is to be preserved, what is your recommended way
> to get a context into the function overloading a method on a non-MATSHELL
> Mat?
>
>     PetscObjectCompose()
>

This means that the user's context has to be a PetscObject. I'm always torn
whether to make my application context a PetscObject. I usually end up doing
so, but it seems to me that it should only be necessary for extraordinay
things. Maybe this qualifies.


>
>    If there is one particular pattern of a non-MATSHELL matrix needing a
> user provided function overloading on a regular basis we need to examine
> that pattern and see what is most appropriate.
>

I use it to provide matrix-free implementations of the global operation with
MatNest. Thus I'm not using MatMult_Nest(), but I am taking advantage of
MatGetSubMatrix_Nest() and similar. The reason for the global operation is
that it makes one traversal of the mesh instead of one traversal per
"block". I don't want to depend on an assembled format because assembling
the true Jacobian is very expensive and the assembled matrix is slow to work
with.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20110606/b8aa3c4f/attachment.html>