[petsc-users] Question about using MatCreateAIJ
Barry Smith
bsmith at petsc.dev
Sun Jun 18 10:22:05 CDT 2023
I am concerned this is not good advice being provided.
Let's back up and look more closely at your use case.
* What is the ratio of new nonzero locations added compared to the initial number of nonzeros for your code, for each of your iterations?
* Is it possible for many iterations, no or very few new nonzeros are being added?
* Are many previous nonzero values becoming zero (or unneeded) later? Again as a ratio compared to the initial number of nonzeros?
* Can you quantify the difference in time between initially filling the matrix and then refilling it using the reset preallocation and not using the reset preallocation?
The effect you report that resetting the preallocation results in slower code is possible if relatively few additional nonzero locations are being created.
After a matrix is assembled with a given nonzero structure (regardless of how it was filled it, using preallocation or not), setting nonzero values into new locations will be slow due to needing to do possibly many mallocs() (as much as one for each new nonzero introduced). Resetting the initially provided preallocation removes the need for all the new mallocs(), but at the expense of needing additional bookkeeping while setting the values. If you did not preallocate originally, then there is no way to prevent the additional mallocs() the second time through, so if you never preallocate but need to add many new nonzero locations adding the new nonzeros will be time consuming; hence in that situation providing the initial preallocation (taking into account all future nonzeros appearing) will pay off in possibly a big way.
I will look at the bug you report for when MatResetPreallocation()is called before the first matrix assembly and see if it can be fixed.
Barry
> On Jun 18, 2023, at 9:01 AM, Junchao Zhang <junchao.zhang at gmail.com> wrote:
>
> yeah, see a C example at https://gitlab.com/petsc/petsc/-/blob/main/src/mat/tests/ex259.c
>
> I guess you can code in this outline with petsc-3.19
>
> MatCreate
> MatSetSizes
> MatSetFromOptions
> iteration-loop start
>
> MatResetPreallocation(A,...)
>
> Fill Matrix A with MatSetValues
>
> MatAssemblyXX(A_WAVE,MAT_FINAL_ASSEMBLY,IERR)
>
> iteration-loop end
>
>
> --Junchao Zhang
>
>
> On Sun, Jun 18, 2023 at 7:37 AM Matthew Knepley <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
>> On Sun, Jun 18, 2023 at 4:26 AM Karsten Lettmann <karsten.lettmann at uni-oldenburg.de <mailto: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 )
>>>
>>>
>>>
>>> <fgcfkbigjgelpmgo.png>
>>>
>>>
>>>
>>>
>>>
>>> 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 <mailto: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 <mailto: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 <mailto: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 <mailto: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/6c9e2ec2/attachment-0001.html>
More information about the petsc-users
mailing list