<html><head>
<meta http-equiv="Content-Type" content="text/html; charset=Windows-1252">
</head>
<body>
<p>Thanks to both of you! <br>
</p>
@Matthew: <br>
<p> Indeed my setJacobianPattern() makes the calls to
MatXAIJSetPreallocation, but does not insert zeros in the matrix.
</p>
<p> I had missed that actual insertions of zeros was needed
before calling SNESSolve.<br>
</p>
<p>@Barry: <br>
</p>
<p> Good idea: I'll study your DMCreateMatrix routines.<br>
</p>
<p>Thanks</p>
<p> Matteo<br>
</p>
<div class="moz-cite-prefix">Il 07/09/21 03:48, Barry Smith ha
scritto:<br>
</div>
<blockquote type="cite" cite="mid:9EE66D0B-604F-4683-A42A-D487223DEE0B@petsc.dev">
<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="" moz-do-not-send="true">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="" moz-do-not-send="true">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="" moz-do-not-send="true">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 class="" clear="all">
<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="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" originalsrc="http://www.cse.buffalo.edu/~knepley/" shash="cvscZ81RV4MSK5H6vFiM83zydoSYtY6cCmo04oh2y3qSR4XsuwWVCN9F2VCToslZEHBQ7NGoJFKGxbZw7yr3BDqf8oarM5hb6tipAu5TGhA4ps48bt0JL0wRdi9olMkOgNlLryybL2vGC/LCQDy9ewWBAx7ytdvXeiL9nOGIncA=" target="_blank" class="" moz-do-not-send="true">https://www.cse.buffalo.edu/~knepley/</a><br class="">
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</blockquote>
<pre class="moz-signature" cols="72">--
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</pre>
</body>
</html>