[petsc-users] Error when using TSSetRHSJacobian+DMSetMatrixPreallocateOnly

[Vijay S. Mahadevan]

>    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.

> 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);
>>>
>