"Diagonal" Columns of MPIAIJ Matrix?
Andrew Colombi
acolombi at gmail.com
Wed Jul 2 16:33:41 CDT 2008
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
>>>>>
>>>>>
>>>>
>>>
>>
>
>
More information about the petsc-users
mailing list