[petsc-users] Question about using MatCreateAIJ
Matthew Knepley
knepley at gmail.com
Sun Jun 18 07:37:27 CDT 2023
On Sun, Jun 18, 2023 at 4:26 AM Karsten Lettmann <
karsten.lettmann at uni-oldenburg.de> wrote:
> Dear Matthew,
>
>
> thanks for you help.
>
>
> 1) I tested your suggestion to pass NULL as the arguments for the
> MatXAIJSetPreallocation.
>
> So the old version was:
>
> CALL
> MatCreateMPIAIJ(MPI_GROUP,N_local,N_local,N_global,N_global,0,DNNZ,0,ONNZ,A,IERR)
>
> And after you suggestion it is now:
>
> CALL MatCreate(MPI_GROUP,A,IERR)
> CALL MatSetType(A,MATAIJ,IERR)
> CALL MatSetSizes(A,N_local,N_local,N_global,N_global,IERR)
> CALL
> MatXAIJSetPreallocation(A,1,DNNZ,ONNZ,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,IERR)
>
>
> Setting block-size = 1.
>
>
> 2) Concerning the error with MatResetPreallocation:
>
> We have an iterative loop, in which the matrix A is filled very often with
> different non-zero structure.
> 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
> on-demand allocations.
>
>
> So I did the following coding in principle:
>
>
> Set matrix A the first time with preallocation
>
>
> iteration-loop start
>
> MatResetPreallocation(A,...)
>
> MatZeroEntries (A)
>
> Fill Matrix A
>
> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)
>
>
> iteration-loop end
>
>
> With these settings, the code run with 2 CPU.
> But with 1 CPU I got an error, which was in MatResetPreallocation.
> I could not understand, why the above code works with 2 CPU but not with 1
> CPU.
>
> 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.
>
> (Now, an image is included showing the codings of :
>
> https://petsc.org/release/src/mat/impls/aij/seq/aij.c.html#MatResetPreallocation_SeqAIJ
>
> https://petsc.org/release/src/mat/impls/aij/mpi/mpiaij.c.html#MatResetPreallocation_MPIAIJ
> )
>
>
>
>
> 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.
> So for one CPU I got an error, while 2 CPUs seem to have been more
> tolerant.
> (I'm not sure, if this interpretation is correct.)
>
>
> So, I changed the coding in that way, that I check the assembly status
> before the preallocation.
>
>
> Using the coding:
>
> CALL MatAssembled(A,A_assembled,ierr)
> IF (A_assembled .eqv. PETSC_TRUE) then
> CALL MatResetPreallocation(A,ierr)
> ENDIF
>
> then worked for 1 and 2 CPU.
>
>
>
> 3) There was another finding, which hopefully is correct.
>
>
> Actually, I did this MatResetPreallocation to have a better performance
> when filling the matrix A later, as was suggested on the manual pages.
>
> 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.
>
> So, in the end, my code seems to be faster, when *not* doing this in the
> iteration loop:
>
> CALL MatAssembled(A,A_assembled,ierr)
> IF (A_assembled .eqv. PETSC_TRUE) then
> CALL MatResetPreallocation(A,ierr)
> ENDIF
>
>
> As I told you, I'm a beginner to PETSC and I do not know, if I have done
> it correctly???
>
I think this may be correct now. We have rewritten Mat so that inserting
values is much more efficient, and
can be done online, so preallocation is not really needed anymore. It is
possible that this default mechanism
is faster than the old preallocation.
I would try the code without preallocation, using the latest release, and
see how it performs.
Thanks,
Matt
> Best, Karsten
>
> The arguments are a combination of the AIJ and SBAIJ arguments. You can
> just pass NULL for the SBAIJ args.
>
>> Then I ran into issues with Resetpreallocation, that I do not understand.
>>
> Can you send the error?
>
> Thanks,
>
> Matt
>
>> I want this, because we have an iterative procedure, where the matrix
>> A_wave and its non-zero structure are changing very often.
>>
>> I try to find the reason for my problem.
>>
>>
>>
>> I really thank you for your answer, that helped me to understand things a
>> bit.
>>
>>
>> I wish you all the best, Karsten
>>
>>
>>
>> Am 15.06.23 um 16:51 schrieb Matthew Knepley:
>>
>> ACHTUNG! Diese E-Mail kommt von Extern! WARNING! This email originated
>> off-campus.
>> On Thu, Jun 15, 2023 at 8:32 AM Karsten Lettmann <
>> karsten.lettmann at uni-oldenburg.de> wrote:
>>
>>> Dear all,
>>>
>>>
>>> I'm quite new to PETSC. So I hope the following questions are not too
>>> stupid.
>>>
>>>
>>> 1) We have a (Fortran) code, that we want to update from an older PETSC
>>> version (petsc.2.3.3-p16) to a newer version.
>>>
>>> Inside the old code, for creating matrices A, there are function calls
>>> of the from:
>>> MatCreateMPIAIJ
>>>
>>> In the reference page for this old version it says:
>>> When calling this routine with a single process communicator, a matrix
>>> of type SEQAIJ is returned.
>>>
>>> So I assume the following behavior of this old routine:
>>> - for N_proc == 1:
>>> a matrix of type SEQAIJ is returned.
>>>
>>> - for N_proc > 1:
>>> a matrix of type MPIAIJ is returned
>>>
>>>
>>>
>>> 2a) So, in the new code, we want to have a similar behavior.
>>>
>>> I found that this function is not present any more in the newer PETSC
>>> versions.
>>>
>>> Instead, one might use: MatCreateAIJ(….)
>>> ( https://petsc.org/release/manualpages/Mat/MatCreateAIJ/ )
>>>
>>> If I understand the reference page of the function correctly, then,
>>> actually, a similar behavior should be expected:
>>>
>>> - for N_proc == 1:
>>> a matrix of type SEQAIJ is returned.
>>>
>>> - for N_proc > 1:
>>> a matrix of type MPIAIJ is returned
>>>
>>>
>>> 2b) However, on the reference page, there is the note:
>>>
>>> It is recommended that one use the MatCreate(), MatSetType() and/or
>>> MatSetFromOptions(), MatXXXXSetPreallocation() paradigm instead of this
>>> routine directly.
>>>
>>> So, if I want the behavior above, it is recommended to code it like
>>> this, isn't it:
>>>
>>> If (N_Proc == 1)
>>>
>>> MatCreate(.. ,A ,...)
>>> MatSetType(…,A, MATSEQAIJ,..)
>>> MatSetSizes(…,A, ..)
>>> MatSeqAIJSetPreallocation(,...A,...)
>>>
>>> else
>>>
>>> MatCreate(.. ,A ,...)
>>> MatSetType(…,A, MATMPIAIJ,..)
>>> MatSetSizes(…,A, ..)
>>> MatMPIAIJSetPreallocation(,...A,...)
>>>
>>
>> You can use
>>
>> MatCreate(comm, &A);
>> MatSetType(A, MATAIJ);
>> MatSetSizes(A, ...);
>> MatXAIJSetPreallocation(A, ...);
>>
>> We recommend this because we would like to get rid of the convenience
>> functions that
>> wrap up exactly this code.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> end
>>>
>>>
>>>
>>> 3) So my questions are:
>>>
>>> - Is my present understanding correct?
>>>
>>> If yes:
>>>
>>> - Why might using MatCreateAIJ(….) for my case not be helpful?
>>>
>>> - So, why is it recommended to use the way 2b) instead of this
>>> MatCreateAIJ(….) ?
>>>
>>>
>>> Best, Karsten
>>>
>>>
>>>
>>>
>>> --
>>> ICBM
>>> Section: Physical Oceanography
>>> Universitaet Oldenburg
>>> Postfach 5634
>>> D-26046 Oldenburg
>>> Germany
>>>
>>> Tel: +49 (0)441 798 4061
>>> email: karsten.lettmann at uni-oldenburg.de
>>>
>>>
>>
>> --
>> 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/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>> --
>> ICBM
>> Section: Physical Oceanography
>> Universitaet Oldenburg
>> Postfach 5634
>> D-26046 Oldenburg
>> Germany
>>
>> Tel: +49 (0)441 798 4061
>> email: karsten.lettmann at uni-oldenburg.de
>>
>>
>
> --
> 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/
> <http://www.cse.buffalo.edu/~knepley/>
>
> --
> ICBM
> Section: Physical Oceanography
> Universitaet Oldenburg
> Postfach 5634
> D-26046 Oldenburg
> Germany
>
> Tel: +49 (0)441 798 4061
> email: karsten.lettmann at uni-oldenburg.de
>
>
--
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/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230618/17c3b232/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fgcfkbigjgelpmgo.png
Type: image/png
Size: 147416 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230618/17c3b232/attachment-0001.png>
More information about the petsc-users
mailing list