[petsc-users] Error when using TSSetRHSJacobian+DMSetMatrixPreallocateOnly

[Matthew Knepley]

On Wed, Sep 18, 2013 at 10:06 AM, Vijay S. Mahadevan <vijay.m at gmail.com>wrote:

> >    Until you can provide the preallocation  you need to simply call
> MatSetUp() and then MatAssemblyBegin/End() on your matrix inside your
> DMCreateMatrix_Moab(). The expectation is that the DM matrix provided is
> "ready to go", we didn't expect to handle matrices that were "incomplete"
> in other parts of the code.
>
> Interesting. Calling MatAssemblyBegin/End() without any MatSetValues*
> wipes out the preallocation information. With -info turned on, I get
> the following for some sample run.
>
> [0] MatAssemblyEnd_SeqAIJ(): Matrix size: 5192 X 5192; storage space:
> 6906304 unneeded,0 used
> [0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0
> [0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 0
>
> I do understand the expectation that matrix from DMCreateMatrix needs
> to be complete but just that traversing through the mesh and setting
> the right NNZ is difficult (I dont think this should exist in the
> library code since this is physics/problem dependent). Somewhat
> manageable for DMDA since it relies on a predefined stencil but I
> couldn't wrap my head around doing this for a generic enough
> implementation. I will think about this some more.
>

DMPlex does it.

   Matt


> > The "user matrix" must have been assembled otherwise it would have
> stopped above.
>
> It did.
>
> Vijay
>
> On Wed, Sep 18, 2013 at 11:17 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > On Sep 18, 2013, at 10:52 AM, "Vijay S. Mahadevan" <vijay.m at gmail.com>
> wrote:
> >
> >>>    Why are you setting this flag?  If you don't set it then the code
> will go through.
> >>
> >> Yes it passes through with DMDA. But with my custom DM implementation
> >> (MOAB based), this fails unless I explicitly set the NNZ pattern in
> >> DMCreateMatrix. I currently don't do that and the only way to
> >> replicate this behavior with DMDA was setting the flag explicitly.
> >
> >    Until you can provide the preallocation  you need to simply call
> MatSetUp() and then MatAssemblyBegin/End() on your matrix inside your
> DMCreateMatrix_Moab(). The expectation is that the DM matrix provided is
> "ready to go", we didn't expect to handle matrices that were "incomplete"
> in other parts of the code.
> >
> >>
> >>>   There is no reason to use a user created matrix since the DM can
> provide the correct one. One of the "main jobs" of DM is to create that
> matrix so we really don't intend people to use a DM but then create the
> matrix themselves.
> >>
> >> Agreed, and that's what I expected when I passed in NULL args instead
> >> of user matrix. But the DM created matrix errors out during Duplicate
> >> while user provided matrix passes through, even though both of them
> >> are unassembled.
> >
> > PetscErrorCode  MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M)
> > {
> >   PetscErrorCode ierr;
> >   Mat            B;
> >   PetscInt       i;
> >
> >   PetscFunctionBegin;
> >   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
> >   PetscValidType(mat,1);
> >   PetscValidPointer(M,3);
> >   if (!mat->assembled)
> SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for
> unassembled matrix");
> >
> > The "user matrix" must have been assembled otherwise it would have
> stopped above.
> >
> >
> >> And this is what is confusing.
> >>
> >> Vijay
> >>
> >> On Wed, Sep 18, 2013 at 10:46 AM, Barry Smith <bsmith at mcs.anl.gov>
> wrote:
> >>>
> >>>   The code did not "crash", it ended with an error message
> >>>
> >>>
> >>> On Sep 18, 2013, at 10:32 AM, "Vijay S. Mahadevan" <vijay.m at gmail.com>
> wrote:
> >>>
> >>>> All,
> >>>>
> >>>> When I call DMSetMatrixPreallocateOnly
> >>>
> >>>    Why are you setting this flag?  If you don't set it then the code
> will go through.
> >>>
> >>>> and then TSSetRHSJacobian
> >>>> without a matrix input argument (NULL instead), then the code crashes
> >>>> during TSSolve and more specifically during MatDuplicate. This is
> >>>> replicable with src/ts/examples/tutorials/ex2.c. Here's the error with
> >>>
> >>>   Jed is reworking the logic here so the code will change.
> >>>
> >>>
> >>>> stack trace.
> >>>>
> >>>> [0]PETSC ERROR: --------------------- Error Message
> >>>>
> >>>>
> >>>> I hope this is not supposed to happen. The fix seems to be to call
> >>>> TSSetRHSJacobian with a user created matrix instead of letting DM
> >>>> create one,
> >>>
> >>>   There is no reason to use a user created matrix since the DM can
> provide the correct one. One of the "main jobs" of DM is to create that
> matrix so we really don't intend people to use a DM but then create the
> matrix themselves.
> >>>
> >>>   Barry
> >>>
> >>>
> >>>> which rectifies this issue. This behavior feels
> >>>> inconsistent and is there a way to fix this.
> >>>>
> >>>> Vijay
> >>>>
> >>>> PS: Just replace the following lines in ex2.c to replicate the error.
> >>>>
> >>>> /*  ierr =
> TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);CHKERRQ(ierr); */
> >>>> ierr =
> TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
> >>>> ierr = DMSetMatrixPreallocateOnly(appctx.da,PETSC_TRUE);CHKERRQ(ierr);
> >>>
> >
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130918/615d1652/attachment.html>