"Diagonal" Columns of MPIAIJ Matrix?

[Barry Smith]

On Jul 2, 2008, at 3:25 PM, Andrew Colombi wrote:

> Is MatGetOwnershipRangeColumn a new function?

    Looks like it.

> because I don't think it
> I have it.  I'm using petsc-2.3.3-p13.  I searched petsc/include for
> RangeColumn and didn't find anything.
>
> I could switch over to MatPreallocateInitialize, but I'm not entirely
> sure how it works.  The docs reference a chapter in the user manual
> about performance,
>
> "See the chapter in the users manual on performance for more details"
>
> But I can't find anything in the PDF manual that mentions
> PreallocateInitialize or PreallocateSet (I searched the text for that
> string).

    The manual pages for these macros are pretty extensive; they are not
discussed in the manual.
>
>
> My first question is, how does MatPreallocateInitialize allocate dnz
> and onz? It says we should not malloc or free them, but it takes
> PetscInt *dnz not PetscInt **dnz?  My understanding of how these
> utilities work is that you do your normal matrix assembly, but use
> MatPreallocateSet instead of MatSetValues.  Then it just does some
> counting and when you're done dnz and onz are correct.  Am I correct?

   Yes. See src/dm/da/src/utils/fdda.c for many examples of its use.

    Barry

>
> Thanks,
> -Andrew
>
> On Wed, Jul 2, 2008 at 1:52 PM, Barry Smith <bsmith at mcs.anl.gov>  
> wrote:
>>
>> MatGetOwnershipRangeColumn() will give you the needed information.  
>> If you are coding from C/C++
>> you might consider using the macro MatPreallocateInitialize() and  
>> friends. These utilities make it
>> somewhat easier, in some circumstances, to count the entries in the  
>> diagonal and off-diagonal blocks.
>>
>>  Barry
>>
>> On Jul 2, 2008, at 11:58 AM, Andrew Colombi wrote:
>>
>>> Hmm... this has not been my experience.  When A is not square I'm  
>>> finding the that the diagonal portion is also not square.  Instead  
>>> it assigns columns to the diagonal portion by evenly distributing  
>>> the columns among processors.  For example, a 4 by 8 would look  
>>> like:
>>>
>>> (please excuse the embedded HTML, I need it to specify a fixed- 
>>> width font)
>>>
>>>     c0  c1  c2  c3  c4  c5  c6  c7
>>> p0 r0  X   X   X   X
>>> p0 r1  X   X   X   X
>>> p1 r2                  X   X   X   X
>>> p1 r3                  X   X   X   X
>>>
>>>
>>> X marks the diagonal portions for each processor.  I have  
>>> attempted to verify this using the program below.  If you run this  
>>> code it will complain that processor 1 must allocate new memory  
>>> for "(0, 3)" (which I am interpreting to mean the local (0, 3),  
>>> otherwise known as (2, 3)).  Further more you can stop the  
>>> complaint by incrementing the off-diagonal preallocation, and  
>>> decrementing the diagonal preallocation.  Note, that specifying 0  
>>> is some special value that isn't 0. If you specify 0 you won't get  
>>> any malloc errors.
>>>
>>> Mat A;
>>> int rows = 4;
>>> int cols = 8;
>>>
>>> if (mpi_rank() == 0)
>>> {
>>>   MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,  
>>> rows, cols,
>>>                   2, PETSC_NULL, 1, PETSC_NULL, &A);
>>>   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);
>>>
>>>   static const int R = 2;
>>>   int set_rows[R] = {0, 1};
>>>   static const int C = 2;
>>>   int set_cols[C] = {0, 1};
>>>   double vals[C*R] = {1., 1., 1., 1.};
>>>   MatSetValues(A, R, set_rows, C, set_cols, vals, ADD_VALUES);
>>> }
>>> else
>>> {
>>>   MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE,  
>>> PETSC_DECIDE,rows, cols,
>>>                   2, PETSC_NULL, 1, PETSC_NULL, &A);
>>>   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);
>>>
>>>   static const int R = 2;
>>>   int set_rows[R] = {2,3};
>>>   static const int C = 2;
>>>   int set_cols[C] = {2,3};
>>>   double vals[C*R] = {1., 1., 1., 1.};
>>>   MatSetValues(A, R, set_rows, C, set_cols, vals, ADD_VALUES);
>>> }
>>>
>>> MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
>>> MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
>>>
>>> On Wed, Jul 2, 2008 at 6:35 AM, Matthew Knepley  
>>> <knepley at gmail.com> wrote:
>>>> By definition, the diagonal portion is square. Thus you retrieve  
>>>> the
>>>> local rows [start, end), and then the diagonal block is
>>>>
>>>> [start, end) x [start, end)
>>>>
>>>> which you can do with
>>>>
>>>> http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatGetOwnershipRange.html
>>>>
>>>>  Matt
>>>>
>>>> On Wed, Jul 2, 2008 at 2:52 AM, Andrew Colombi  
>>>> <acolombi at gmail.com> wrote:
>>>>> I'm letting PETSC_DECIDE how to distribute the rows and columns  
>>>>> of my
>>>>> non-square matrix, and I want to retrieve the columns it's  
>>>>> decided to
>>>>> make the "Diagonal" for a particular process.  How can I do  
>>>>> that?  I'm
>>>>> probably missing something obvious here, but I checked the docs.
>>>>>
>>>>> In case some context helps (or inspires a work-around):
>>>>>
>>>>> I want this information to count the number of d_nnz and o_nnz to
>>>>> preallocate perfectly.  This shouldn't be hard to do, but I need  
>>>>> to
>>>>> know the boundary between off-diagonal and diagonal to do it.
>>>>>
>>>>> Thanks,
>>>>> -Andrew
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>>
>>>
>>
>