> DMPlex does it. I was just thinking about how you managed this in DMPlex. I will take a look at this implementation.
Vijay On Wed, Sep 18, 2013 at 12:10 PM, Matthew Knepley <[email protected]> wrote: > On Wed, Sep 18, 2013 at 10:06 AM, Vijay S. Mahadevan <[email protected]> > 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 <[email protected]> wrote: >> > >> > On Sep 18, 2013, at 10:52 AM, "Vijay S. Mahadevan" <[email protected]> >> > 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 <[email protected]> >> >> wrote: >> >>> >> >>> The code did not "crash", it ended with an error message >> >>> >> >>> >> >>> On Sep 18, 2013, at 10:32 AM, "Vijay S. Mahadevan" <[email protected]> >> >>> 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
