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

[Barry Smith]

  Matteo,

   I think it might be best if you simply took our "standard" DMCreateMatrix for DMDA routine and modified exactly for your needs. You can find the source code in src/dm/impls/da/fdda.c There are a variety of routines such as DMCreateMatrix_DA_2d_MPIAIJ_Fill() DMCreateMatrix_DA_2d_MPIAIJ() etc. 

    Pick the one that matches your needs, copy it and modify to do the exact preallocation and then filling with zeros for interior points, boundary points etc. I don't think "fixing" the incorrect default behavior after the fact will work.

  Barry


> On Sep 6, 2021, at 7:34 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Mon, Sep 6, 2021 at 12:22 PM Matteo Semplice <matteo.semplice at uninsubria.it <mailto:matteo.semplice at uninsubria.it>> wrote:
> 
> Il 31/08/21 17:32, Jed Brown ha scritto:
> > Matteo Semplice <matteo.semplice at uninsubria.it <mailto: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/f669789d/attachment.html>