[petsc-users] Matrix preallocation

[Matthew Knepley]

On Sun, Feb 6, 2022 at 2:39 PM Samuel Estes <samuelestes91 at gmail.com> wrote:

> First of all, thank you so much for the detailed answers!
> That clears up most of my confusion. Just to clarify let me make sure I
> understand everything:
> 1. So it seems that your advice is just to call the
> MatXAIJSetPreallocation() routine (rather than make separate calls to other
> preallocation routines such as MatSeqAIJPreallcation() and
> MatMPIAIJPreallcation())? This will preallocate for Seq and MPI matrices
> (among other types) and the decision will be made at runtime?
>

Yes.


> 2. So the old way was just to call MatSeqAIJPreallcation() and
> MatMPIAIJPreallcation()? I had thought that if you called the sequential
> version for a parallel matrix or vice versa that the program would crash.
> It seems that instead nothing happens, so by calling both you were covered
> in either case. Is this correct? Essentially, one would always execute and
> the other would do nothing?
>

Yes, PETSc works like Objective-C, in that if you call a routine for a
subclass, and the object is of a different subclass, it is just ignored.


> 3. If the dnnzu and onnzu arrays for the upper triangular parts seem
> unnecessary for me I can just call the MatXAIJSetPreallocation() routine
> with those values set to 'NULL'? I'm really just looking for a blend of the
> MatSeqAIJPreallcation() and MatMPIAIJPreallcation() routines so it seems
> that I only need dnnz and onnz arrays.
>

Yes.


> 4. If we ignore the dnnzu and onnzu upper triangular parameters, then the
> MatXAIJSetPreallocation() routine has the same parameters as the
> MatMPIAIJSetPreallocation() routine without the option to give constant
> integer values for diagonal and off-diagonal parts of the matrix so it
> seems clear to me how this routine would work for a parallel matrix when
> the upper triangular parameters are NULL. If the matrix is sequential then
> does MPIXAIJSetPreallocation() just ignore the onnz parameter and just use
> the dnnz parameter as the sole argument for determining preallocation?
>

Yes.

  Thanks,

    Matt


> Thanks again for all the help. I think I understand well enough to use
> this routine now. The questions above are mostly just to clarify and
> double-check that I understood your previous responses correctly.
>
> Sam
>
> On Sat, Feb 5, 2022 at 10:12 AM Mark Adams <mfadams at lbl.gov> wrote:
>
>>
>>
>> On Sat, Feb 5, 2022 at 10:36 AM Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Fri, Feb 4, 2022 at 11:47 PM Samuel Estes <samuelestes91 at gmail.com>
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>> I have a very basic question about matrix preallocation. I am trying to
>>>> use the MatCreate(), MatSetFromOptions(), MatXXXXSetPreallocation()
>>>> paradigm. I thought that I should use the MatXAIJSetPreallocation() routine
>>>> since the code may be run with a SeqAIJ or MPIAIJ matrix but I do not
>>>> understand all of the inputs required for the
>>>> MatXAIJSetPreallocation routine. In particular, the dnnzu and
>>>> onnzu variables don't quite make sense to me. Can these be NULL? I was
>>>> basically just hoping for a routine that would preallocate for either a
>>>> sequential or parallel matrix depending on what was given at runtime. This
>>>> routine seems to be what I want but I don't understand it very well and the
>>>> documentation hasn't helped me to figure it out.
>>>>
>>>
>>> The example for this is here
>>>
>>>
>>> https://petsc.org/main/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html#MatMPIAIJSetPreallocation
>>>
>>> Maybe we should copy it to the XAIJ page as well. Does this help explain
>>> the arguments?
>>>
>>>
>>>> A related followup question: Is it good practice to use this function
>>>> or should I just use the other routines like MatSeqAIJSetPreallocation()
>>>> and MatMPIAIJSetPreallocation()?
>>>>
>>>> And finally my last question: if I were to use the
>>>> MatSeqAIJSetPreallocation()/MatMPIAIJSetPreallocation() routines for
>>>> preallocating memory, is it common to just call MatGetType() then call the
>>>> appropriate routine depending on whether or not the matrix is parallel or
>>>> not? I ask because when I have tested these routines out, it seems
>>>> that MatSeqAIJSetPreallocation() works even for parallel matrices which is
>>>> a bit confusing. I'm assuming that it just sets the diagonal part of the
>>>> matrix?
>>>>
>>>
>>> No, it definitely will not preallocate in parallel, so something else is
>>> happening.
>>>
>>
>> With only MatSeqAIJSetPreallocation in parallel it would not do any
>> preallocation, in which case it would fall back to dynamic allocation. I'm
>> guessing that is what is happening and it will be very slow.
>>
>>
>>>
>>>   Thanks,
>>>
>>>      Matt
>>>
>>>
>>>> I hope that my questions were clear. Let me know if they need
>>>> clarification and thanks in advance for the help!
>>>>
>>>> Sam
>>>>
>>>
>>>
>>> --
>>> 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/>
>>>
>>

-- 
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/20220206/e3e75e22/attachment-0001.html>