[petsc-users] Matrix preallocation

[Mark Adams]

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?
>

Yes. The integer is ignored if you provide the array. Just make sure that
the integer is an upper bound on the number of non-zeros in any row. If
not, MatSetValues can be very slow. Otherwise, the performance hit is not
bad.


> 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.
>

That is what the "X" is for. It is basically syntactic sugar to call the
right version of preallocate (seq, mpi, blocked versions, hypre, etc.)


> 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.
>

Sorry,


>
> A related followup question: Is it good practice to use this function or
> should I just use the other routines like MatSeqAIJSetPreallocation() and
> MatMPIAIJSetPreallocation()?
>

Yes this is the way to go.
TL;DR
Before we had a sparse matrix type explosion with GPUs, you could call both
the MPI and Seq versions and be fine. You could call the blocked versions I
suppose to be safe.
But now there are several GPU/device enabled matrix types and "X" is
convenient.


>
> 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?
>

You could do that but no need. As I said, the old way was to call both.


> 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?
>

There is probably a default. It should work w/o any of these calls, but
performance would suffer if the default (10 maybe) is too small for you.

Mark


> I hope that my questions were clear. Let me know if they need
> clarification and thanks in advance for the help!
>
> Sam
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220205/cc5f01bc/attachment.html>