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
