<div dir="ltr">yeah, see a C example at <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex259.c">https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex259.c</a><div><br></div><div>I guess you can code in this outline with petsc-3.19</div><div><br></div><div><code>MatCreate</code><br></div><div><code><code>MatSetSizes</code><br></code></div><div><code><code><code>MatSetFromOptions</code><br></code></code></div><div><p><font face="monospace">iteration-loop start</font></p><p><font face="monospace"> MatResetPreallocation(A,...)</font></p><p><font face="monospace"> Fill Matrix A with MatSetValues</font></p><p><font face="monospace"> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)</font></p><p><font face="monospace"> iteration-loop end</font><br></p><div><br clear="all"><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr">--Junchao Zhang</div></div></div><br></div></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sun, Jun 18, 2023 at 7:37 AM Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Sun, Jun 18, 2023 at 4:26 AM Karsten Lettmann <<a href="mailto:karsten.lettmann@uni-oldenburg.de" target="_blank">karsten.lettmann@uni-oldenburg.de</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">
<div>
<p>Dear Matthew,</p>
<p><br>
</p>
<p>thanks for you help.</p>
<p><br>
</p>
<p>1) I tested your suggestion to pass NULL as the arguments for the
MatXAIJSetPreallocation.</p>
<p>So the old version was:</p>
<p>CALL
MatCreateMPIAIJ(MPI_GROUP,N_local,N_local,N_global,N_global,0,DNNZ,0,ONNZ,A,IERR)<br>
</p>
<p>And after you suggestion it is now:</p>
<p> CALL MatCreate(MPI_GROUP,A,IERR)<br>
CALL MatSetType(A,MATAIJ,IERR)<br>
CALL MatSetSizes(A,N_local,N_local,N_global,N_global,IERR)<br>
CALL
MatXAIJSetPreallocation(A,1,DNNZ,ONNZ,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,IERR)<br>
</p>
<p><br>
</p>
<p> Setting block-size = 1.<br>
</p>
<p><br>
</p>
<p>2) Concerning the error with MatResetPreallocation:</p>
<p>We have an iterative loop, in which the matrix A is filled very
often with different non-zero structure.<br>
Further, I read in the manual pages that due to performance
issues, one should preallocate enough space, as operations as
matsetvalues might be time consuming due to additional <br>
on-demand allocations.</p>
<p><br>
</p>
<p>So I did the following coding in principle:</p>
<p><br>
</p>
<p> Set matrix A the first time with preallocation</p>
<p><br>
</p>
<p> iteration-loop start</p>
<p> MatResetPreallocation(A,...)</p>
<p> MatZeroEntries (A)<br>
</p>
<p> Fill Matrix A</p>
<p> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)</p>
<p><br>
</p>
<p> iteration-loop end<br>
</p>
<p><br>
</p>
<p>With these settings, the code run with 2 CPU.<br>
But with 1 CPU I got an error, which was in MatResetPreallocation.<br>
I could not understand, why the above code works with 2 CPU but
not with 1 CPU.</p>
<p>At the moment, I believe the reason for this error seems to be a
pre-check, that is done in SeqAIJ but not in MPIAIJ fo a valid and
present matrix A.</p>
<p>(Now, an image is included showing the codings of :<br>
<a href="https://petsc.org/release/src/mat/impls/aij/seq/aij.c.html#MatResetPreallocation_SeqAIJ" target="_blank">https://petsc.org/release/src/mat/impls/aij/seq/aij.c.html#MatResetPreallocation_SeqAIJ</a><br>
<a href="https://petsc.org/release/src/mat/impls/aij/mpi/mpiaij.c.html#MatResetPreallocation_MPIAIJ" target="_blank">https://petsc.org/release/src/mat/impls/aij/mpi/mpiaij.c.html#MatResetPreallocation_MPIAIJ</a>
) <br>
</p>
<p><br>
</p>
<p><img src="cid:188ce81459a2eb3091d1" alt=""></p>
<p><br>
</p>
<p><br>
</p>
<p>So, it seems for me at the moment, that the first
MatResetPreallocation (when the iteration loop is entered the
first time) is done on an not-assembled matrix A.<br>
So for one CPU I got an error, while 2 CPUs seem to have been more
tolerant.<br>
(I'm not sure, if this interpretation is correct.)<br>
</p>
<p><br>
</p>
<p>So, I changed the coding in that way, that I check the assembly
status before the preallocation.</p>
<p><br>
</p>
<p>Using the coding:</p>
<p> CALL MatAssembled(A,A_assembled,ierr)<br>
IF (A_assembled .eqv. PETSC_TRUE) then<br>
CALL MatResetPreallocation(A,ierr)<br>
ENDIF <br>
</p>
<p>then worked for 1 and 2 CPU.<br>
</p>
<p><br>
</p>
<p><br>
</p>
<p>3) There was another finding, which hopefully is correct.<br>
</p>
<p><br>
</p>
<p>Actually, I did this MatResetPreallocation to have a better
performance when filling the matrix A later, as was suggested on
the manual pages.<br>
<br>
</p>
<p>However, I found (if I did nothing wrong) that this
MatResetPreallocation was much more time consuming than the
additional (and unwanted) allocations done during filling the
matrix.<br>
<br>
</p>
<p>So, in the end, my code seems to be faster, when <b>not</b>
doing this in the iteration loop:<br>
</p>
<p> CALL MatAssembled(A,A_assembled,ierr)<br>
IF (A_assembled .eqv. PETSC_TRUE) then<br>
CALL MatResetPreallocation(A,ierr)<br>
ENDIF </p>
<p><br>
</p>
<p>As I told you, I'm a beginner to PETSC and I do not know, if I
have done it correctly???</p></div></blockquote><div>I think this may be correct now. We have rewritten Mat so that inserting values is much more efficient, and</div><div>can be done online, so preallocation is not really needed anymore. It is possible that this default mechanism</div><div>is faster than the old preallocation.</div><div><br></div><div>I would try the code without preallocation, using the latest release, and see how it performs.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>
<p>Best, Karsten</p>
<br>
<blockquote type="cite">
<div>
<div dir="ltr">
<div class="gmail_quote">
<div>The arguments are a combination of the AIJ and SBAIJ
arguments. You can just pass NULL for the SBAIJ args.<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<p>Then I ran into issues with Resetpreallocation, that
I do not understand.</p>
</div>
</blockquote>
<div>Can you send the error?</div>
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt <br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<p>I want this, because we have an iterative procedure,
where the matrix A_wave and its non-zero structure are
changing very often.<br>
</p>
<p>I try to find the reason for my problem.</p>
<p><br>
</p>
<p><br>
</p>
<p>I really thank you for your answer, that helped me to
understand things a bit.<br>
</p>
<p><br>
</p>
<p>I wish you all the best, Karsten<br>
</p>
<p><br>
</p>
<p><br>
</p>
<div>Am 15.06.23 um 16:51 schrieb Matthew Knepley:<br>
</div>
<blockquote type="cite">
<div style="border:1pt solid rgb(0,51,51);padding:2pt">
<p style="line-height:10pt;background:rgb(213,234,255)"><span>ACHTUNG!</span><span> Diese E-Mail kommt von Extern!
<span style="color:rgb(255,0,0)">WARNING!</span>
This email originated off-campus.<br>
</span></p>
</div>
<div>
<div dir="ltr">
<div dir="ltr">On Thu, Jun 15, 2023 at 8:32 AM
Karsten Lettmann <<a href="mailto:karsten.lettmann@uni-oldenburg.de" target="_blank">karsten.lettmann@uni-oldenburg.de</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">
Dear all,<br>
<br>
<br>
I'm quite new to PETSC. So I hope the
following questions are not too <br>
stupid.<br>
<br>
<br>
1) We have a (Fortran) code, that we want to
update from an older PETSC <br>
version (petsc.2.3.3-p16) to a newer version.<br>
<br>
Inside the old code, for creating matrices A,
there are function calls <br>
of the from:<br>
MatCreateMPIAIJ<br>
<br>
In the reference page for this old version it
says:<br>
When calling this routine with a single
process communicator, a matrix <br>
of type SEQAIJ is returned.<br>
<br>
So I assume the following behavior of this old
routine:<br>
- for N_proc == 1:<br>
a matrix of type SEQAIJ is returned.<br>
<br>
- for N_proc > 1:<br>
a matrix of type MPIAIJ is returned<br>
<br>
<br>
<br>
2a) So, in the new code, we want to have a
similar behavior.<br>
<br>
I found that this function is not present any
more in the newer PETSC <br>
versions.<br>
<br>
Instead, one might use: MatCreateAIJ(….)<br>
( <a href="https://petsc.org/release/manualpages/Mat/MatCreateAIJ/" rel="noreferrer" target="_blank">
https://petsc.org/release/manualpages/Mat/MatCreateAIJ/</a> )<br>
<br>
If I understand the reference page of the
function correctly, then, <br>
actually, a similar behavior should be
expected:<br>
<br>
- for N_proc == 1:<br>
a matrix of type SEQAIJ is returned.<br>
<br>
- for N_proc > 1:<br>
a matrix of type MPIAIJ is returned<br>
<br>
<br>
2b) However, on the reference page, there is
the note:<br>
<br>
It is recommended that one use the
MatCreate(), MatSetType() and/or <br>
MatSetFromOptions(), MatXXXXSetPreallocation()
paradigm instead of this <br>
routine directly.<br>
<br>
So, if I want the behavior above, it is
recommended to code it like <br>
this, isn't it:<br>
<br>
If (N_Proc == 1)<br>
<br>
MatCreate(.. ,A ,...)<br>
MatSetType(…,A, MATSEQAIJ,..)<br>
MatSetSizes(…,A, ..)<br>
MatSeqAIJSetPreallocation(,...A,...)<br>
<br>
else<br>
<br>
MatCreate(.. ,A ,...)<br>
MatSetType(…,A, MATMPIAIJ,..)<br>
MatSetSizes(…,A, ..)<br>
MatMPIAIJSetPreallocation(,...A,...)<br>
</blockquote>
<div><br>
</div>
<div>You can use</div>
<div><br>
</div>
<div> MatCreate(comm, &A);</div>
<div> MatSetType(A, MATAIJ);</div>
<div> MatSetSizes(A, ...);</div>
<div> MatXAIJSetPreallocation(A, ...);</div>
<div><br>
</div>
<div>We recommend this because we would like to
get rid of the convenience functions that</div>
<div>wrap up exactly this code.</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">
end<br>
<br>
<br>
<br>
3) So my questions are:<br>
<br>
- Is my present understanding correct?<br>
<br>
If yes:<br>
<br>
- Why might using MatCreateAIJ(….) for my case
not be helpful?<br>
<br>
- So, why is it recommended to use the way 2b)
instead of this <br>
MatCreateAIJ(….) ?<br>
<br>
<br>
Best, Karsten<br>
<br>
<br>
<br>
<br>
-- <br>
ICBM<br>
Section: Physical Oceanography<br>
Universitaet Oldenburg<br>
Postfach 5634<br>
D-26046 Oldenburg<br>
Germany<br>
<br>
Tel: +49 (0)441 798 4061<br>
email: <a href="mailto:karsten.lettmann@uni-oldenburg.de" target="_blank">karsten.lettmann@uni-oldenburg.de</a><br>
<br>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
<span class="gmail_signature_prefix">-- </span><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>
</div>
</blockquote>
<pre cols="72">--
ICBM
Section: Physical Oceanography
Universitaet Oldenburg
Postfach 5634
D-26046 Oldenburg
Germany
Tel: +49 (0)441 798 4061
email: <a href="mailto:karsten.lettmann@uni-oldenburg.de" target="_blank">karsten.lettmann@uni-oldenburg.de</a></pre>
</div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
<span class="gmail_signature_prefix">-- </span><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>
</div>
</blockquote>
<pre cols="72">--
ICBM
Section: Physical Oceanography
Universitaet Oldenburg
Postfach 5634
D-26046 Oldenburg
Germany
Tel: +49 (0)441 798 4061
email: <a href="mailto:karsten.lettmann@uni-oldenburg.de" target="_blank">karsten.lettmann@uni-oldenburg.de</a></pre>
</div>
</blockquote></div><br clear="all"><div><br></div><span class="gmail_signature_prefix">-- </span><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>
</blockquote></div>