[petsc-users] Mat preallocation for SNES jacobian [WAS Re: Mat preallocation in case of variable stencils]
Matteo Semplice
matteo.semplice at uninsubria.it
Tue Sep 7 03:01:39 CDT 2021
Thanks to both of you!
@Matthew:
Indeed my setJacobianPattern() makes the calls to
MatXAIJSetPreallocation, but does not insert zeros in the matrix.
I had missed that actual insertions of zeros was needed before
calling SNESSolve.
@Barry:
Good idea: I'll study your DMCreateMatrix routines.
Thanks
Matteo
Il 07/09/21 03:48, Barry Smith ha scritto:
> 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
>> <mailto: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/
>> <https://eur01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.cse.buffalo.edu%2F~knepley%2F&data=04%7C01%7Cmatteo.semplice%40uninsubria.it%7C181ba7fff01f4de6308608d971a1966e%7C9252ed8bdffc401c86ca6237da9991fa%7C0%7C0%7C637665761093056131%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=NMMtdL%2Bt78kPAV3%2Fce1RCpzD4uA%2FBIElk5qJRWRaYjs%3D&reserved=0>
>
--
Prof. Matteo Semplice
Università degli Studi dell’Insubria
Dipartimento di Scienza e Alta Tecnologia – DiSAT
Professore Associato
Via Valleggio, 11 – 22100 Como (CO) – Italia
tel.: +39 031 2386316
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210907/87449ab6/attachment-0001.html>
More information about the petsc-users
mailing list