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