[petsc-users] Allocating the diagonal for MatMPIAIJSetPreallocation

Fri Apr 1 12:44:29 CDT 2022

On Fri, Apr 1, 2022 at 1:08 PM Samuel Estes <samuelestes91 at gmail.com> wrote:

> Thank you for all the help. I think I'm mostly clear now.
>
> I think my only remaining question relates to point 1:
> So let's call the global matrix across all processors A and let's say it's
> a square matrix of size n. Then each processor has a subset of the rows of
> A so that there's in principle a contiguous rectangular chunk of A on each
> proc. Now we want to use this preallocator matrix to determine the non-zero
> structure of each processors chunk of A. Let's call this matrix P. You say
> that P should also be square. But I would think that, at least locally, P
> should be rectangular on each proc and basically have the same size as the
> local chunks of A. Do you mean that P is square when considered globally
> across all processors? I guess my question is a bit subtle, but basically,
> is P essentially a parallel Mat type which exists nominally across all
> procs, or is it a serial type object which is local to each proc? I think I
> was viewing it in the second way since it doesn't obviously require
> communication but I guess if you consider it in the first way it makes
> sense that it would be nxn.
>

The idea here is that the internal structure of P does not matter. it has
the same interface as the matrix A, so fom your point of view they are
identical.

  Thanks,

     Matt

> On Fri, Apr 1, 2022 at 12:00 PM Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Fri, Apr 1, 2022 at 12:57 PM Samuel Estes <samuelestes91 at gmail.com>
>> wrote:
>>
>>> Ah, so let me see if I understand this now. So basically you preallocate
>>> for this "dummy matrix" of type MatPreallocator. But rather than simply
>>> calling MatXXXAIJSetPreallocation() to allocate space, you actually tell
>>> this "dummy matrix" exactly where the non-zeros will be. You then call this
>>> MatPreallocatorPreallocate() routine to pass the non-zero structure of the
>>> "dummy matrix" to the actual matrix that you're using for the linear
>>> system. You then call MatSetValues() (or a variant thereof) to set the
>>> values of this preallocated matrix.
>>>
>>> A couple more clarifying questions:
>>> 1. So this MatPreallocator matrix should have a rectangular structure of
>>> (rEnd-rStart) rows and n columns (where n is the total size of the system?
>>> I'm assuming here that the global matrix across all processors is square.
>>>
>>
>> This "dummy matrix" will look exactly like your system matrix, meaning it
>> is square, and you tell it the nonzeros by just calling MatSetValues(),
>> exactly as you would for your system matrix. That is why you can reuse
>> exactly the same code.
>>
>>
>>> 2. You preallocate for the dummy matrix by calling matsetvalues rather
>>> than the usual way of preallocating?
>>>
>>
>> Yes.
>>
>>
>>> 3. That's what you mean by looping twice? You loop over once using
>>> matsetvalues to preallocate, then you have to loop using matsetvalues again
>>> a second time to actually set the values of the parallel Mat you actually
>>> use to solve the system?
>>>
>>
>> Yes.
>>
>>   Thanks,
>>
>>     Matt
>>
>>
>>> On Fri, Apr 1, 2022 at 11:50 AM Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Fri, Apr 1, 2022 at 12:45 PM Samuel Estes <samuelestes91 at gmail.com>
>>>> wrote:
>>>>
>>>>> Thanks! This seems like it might be what I need. I'm still a little
>>>>> unclear on how it works though? My problem is basically that for any given
>>>>> row, I know the total non-zeros but not how many occur in the diagonal vs
>>>>> off-diagonal. Without knowledge of the underlying grid, I'm not sure how
>>>>> there could be a black box utility to figure this out? Am I
>>>>> misunderstanding how this is used?
>>>>>
>>>>
>>>> So each process gets a stack of rows, [rStart, rEnd). A nonzero (r, c)
>>>> is in the diagonal block if rStart <= c < rEnd. So if you know (r, c) for
>>>> each nonzero, you know whether it is in the diagonal block.
>>>>
>>>>   Thanks,
>>>>
>>>>      Matt
>>>>
>>>>
>>>>> On Fri, Apr 1, 2022 at 11:34 AM Matthew Knepley <knepley at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> On Fri, Apr 1, 2022 at 12:27 PM Samuel Estes <samuelestes91 at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> I have a problem in which I know (roughly) the number of non-zero
>>>>>>> entries for each row of a matrix but I don't have a convenient way of
>>>>>>> determining whether they belong to the diagonal or off-diagonal part of the
>>>>>>> parallel matrix. Is there some way I can just allocate the total number of
>>>>>>> non-zeros in a row regardless of which part they belong to? I'm assuming
>>>>>>> that this is not possible but I just wanted to check. It seems like it
>>>>>>> should be possible in principle since the matrix is only split by the rows
>>>>>>> but the columns of a row are all on the same processor (at least as I
>>>>>>> understand it). Thanks!
>>>>>>>
>>>>>>
>>>>>> In serial, the matrix is stored by rows. In parallel, it is split
>>>>>> into a diagonal and off-diagonal block, so that we can overlap
>>>>>> communication and computation in the matvec.
>>>>>>
>>>>>> However, we have a convenience structure for figuring this out,
>>>>>> called MatPreallocator,
>>>>>> https://petsc.org/main/docs/manualpages/Mat/MATPREALLOCATOR.html
>>>>>> In my code, I wrote a loop around the code that filled up my matrix,
>>>>>> which executed twice. On the first pass, I fed in the MatPreallocator
>>>>>> matrix. When this finished
>>>>>> you can call
>>>>>> https://petsc.org/main/docs/manualpages/Mat/MatPreallocatorPreallocate.html#MatPreallocatorPreallocate
>>>>>> on your system amtrix, then on the second
>>>>>> pass use the system matrix. This was only a few extra lines of code
>>>>>> for me. If you want to optimize further, you can have a flag that only
>>>>>> computes the values on the
>>>>>> second pass.
>>>>>>
>>>>>>   Thanks,
>>>>>>
>>>>>>       Matt
>>>>>>
>>>>>> 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/>
>>>>
>>>
>>
>> --
>> 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/20220401/27fb88be/attachment-0001.html>