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