On Mon, Sep 6, 2021 at 12:22 PM Matteo Semplice < [email protected]> wrote:
> > Il 31/08/21 17:32, Jed Brown ha scritto: > > Matteo Semplice <[email protected]> writes: > > > >> Hi. > >> > >> We are writing a code for a FD scheme on an irregular domain and thus > >> the local stencil is quite variable: we have inner nodes, boundary nodes > >> and inactive nodes, each with their own stencil type and offset with > >> respect to the grid node. We currently create a matrix with > >> DMCreateMatrix on a DMDA and for now have set the option > >> MAT_NEW_NONZERO_LOCATIONS to PETSC_TRUE, but its time to render the code > >> memory-efficient. The layout created automatically is correct for inner > >> nodes, wrong for boundary ones (off-centered stencils) and redundant for > >> outer nodes. > >> > >> After the preprocessing stage (including stencil creation) we'd be in > >> position to set the nonzero pattern properly. > >> > >> Do we need to start from a Mat created by CreateMatrix? Or is it ok to > >> call DMCreateMatrix (so that the splitting among CPUs and the block size > >> are set by PETSc) and then call a MatSetPreallocation routine? > > You can call MatXAIJSetPreallocation after. It'll handle all matrix > types so you don't have to shepherd data for all the specific > preallocations. > > Hi. > > Actually I am still struggling with this... Let me explain. > > My code relies on a SNES and the matrix I need to preallocate is the > Jacobian. So I do: > > in the main file > ierr = DMCreateMatrix(ctx.daAll,&ctx.J);CHKERRQ(ierr); > ierr = setJacobianPattern(ctx,ctx.J);CHKERRQ(ierr); //calling > MatXAIJSetPreallocation on the second argument > I do not understand this. DMCreateMatrix() has already preallocated _and_ filled the matrix with zeros. Additional preallocation statements will not do anything (I think). > ierr = MatSetOption(ctx.J,MAT_NEW_NONZERO_LOCATIONS,*******); > CHKERRQ(ierr);//allow new nonzeros? > ierr = SNESSetFunction(snes,ctx.F ,FormFunction,(void *) &ctx); > CHKERRQ(ierr); > ierr = SNESSetJacobian(snes,ctx.J,ctx.J,FormJacobian,(void *) &ctx); > CHKERRQ(ierr); > > ierr = FormSulfationRHS(ctx, ctx.U0, ctx.RHS);CHKERRQ(ierr); > ierr = SNESSolve(snes,ctx.RHS,ctx.U); CHKERRQ(ierr); > > and > > PetscErrorCode FormJacobian(SNES snes,Vec U,Mat J, Mat P,void *_ctx) > > does (this is a 2 dof finite difference problem, so logically 2x2 blocks > in J) > > ierr = setJac00(....,P) //calls to MatSetValues in the 00 block > > ierr = setJac01(....,P) //calls to MatSetValues in the 01 block > > ierr = setJac1X(....,P) //calls to MatSetValues in the 10 ans 11 block > > ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > > If I run with MAT_NEW_NONZERO_LOCATIONS=TRUE, all runs fine and using > the -info option I see that no mallocs are performed during Assembly. > > Computing F > 0 SNES Function norm 7.672682917097e+02 > Computing J > [0] MatAssemblyEnd_SeqAIJ(): Matrix size: 71874 X 71874; storage space: > 17661 unneeded,191714 used > [0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0 > [0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 27 > [0] MatCheckCompressedRow(): Found the ratio (num_zerorows > 0)/(num_localrows 71874) < 0.6. Do not use CompressedRow routines. > > If I omit the call to setJacobianPattern, info reports a nonsero number > of mallocs, so somehow the setJacobianPattern routine should be doing > its job correctly. > Hmm, this might mean that the second preallocation call is wiping out the info in the first. Okay, I will go back and look at the code. > However, if I run with MAT_NEW_NONZERO_LOCATIONS=FALSE, the Jacobian is > entirely zero and no error messages appear until the KSP tries to do its > job: > This sounds like your setJacobianPattern() is not filling the matrix with zeros, so that the insertions make new nonzeros. It is hard to make sense of this string of events without the code. > Computing F > 0 SNES Function norm 7.672682917097e+02 > Computing J > [0] MatAssemblyEnd_SeqAIJ(): Matrix size: 71874 X 71874; storage space: > 209375 unneeded,0 used > [0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0 > [0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 0 > [0] MatCheckCompressedRow(): Found the ratio (num_zerorows > 71874)/(num_localrows 71874) > 0.6. Use CompressedRow routines. > ... and then KSP complains! > > I have tried adding MAT_FLUSH_ASSEMBLY calls inside the subroutines, but > nothing changes. > > So I have 2 questions: > > 1. If, as a temporary solution, I leave MAT_NEW_NONZERO_LOCATIONS=TRUE, > am I going to incur in performance penalties even if no new nonzeros are > created by my assembly routine? > If you are worried about performance, I think the option you want is MAT_NEW_NONZERO_ALLOCATION_ERR since allocations, not new nonzeros, are what causes performance problems. Thanks, Matt > 2. Can you guess what is causing the problem? > > Thanks > > Matteo > > > > > > > -- 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 https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
