On Sep 6, 2021, at 7:34 PM, Matthew Knepley <[email protected]
<mailto:[email protected]>> wrote:
On Mon, Sep 6, 2021 at 12:22 PM Matteo Semplice
<[email protected]
<mailto:[email protected]>> wrote:
Il 31/08/21 17:32, Jed Brown ha scritto:
> Matteo Semplice <[email protected]
<mailto:[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/
<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>
