"Diagonal" Columns of MPIAIJ Matrix?

[Barry Smith]

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