[petsc-users] Mat preallocation for SNES jacobian [WAS Re: Mat preallocation in case of variable stencils]

[Matthew Knepley]

On Mon, Sep 6, 2021 at 12:22 PM Matteo Semplice <
matteo.semplice at uninsubria.it> wrote:

>
> Il 31/08/21 17:32, Jed Brown ha scritto:
> > Matteo Semplice <matteo.semplice at uninsubria.it> 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210906/196a816c/attachment.html>