"Diagonal" Columns of MPIAIJ Matrix?

[Barry Smith]

   The code used to divide it up by default is PetscMapSetUp() in src/ 
vec/impls/mpi/pmap.c


    Barry

On Jul 2, 2008, at 4:33 PM, Andrew Colombi wrote:

> Is it fair to say that MatPreallocateInitialize and MatPreallocateSet
> need to be in the same scope for this to work?  If so, that will never
> work with my code. (Unless I can call repeatedly without having it
> reinitialize things.)
>
> I already had everything setup to work with my own bookkeeping (STL
> vector and binary_search, in my tests STL set took too much memory).
> But I need get the diagonal columns!
>
> So be it.  I'm willing to hack.  May I assume PETSc distributes
> columns with the following algorithm:
>
> int diag_cols = total_cols / processors;
> if (mpi_rank() < total_cols % processors)
>    diag_cols += 1;
>
> Then assume the offsets follow the MPI rank? (i.e. higher rank has
> later columns)
>
> Thanks,
> -Andrew
>
> On Wed, Jul 2, 2008 at 3:34 PM, Barry Smith <bsmith at mcs.anl.gov>  
> wrote:
>>
>> 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
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>>
>