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:<br>
<br>(please excuse the embedded HTML, I need it to specify a fixed-width font)<br> <br><span style="font-family: courier new,monospace;"> </span><font face="courier new,monospace">c0 c1 c2 c3 c4 c5 c6 c7</font><br>
<font face="courier new,monospace">p0 r0 X X X X<br>p0 r1 X X X X<br>p1 r2 X X X X<br>p1 r3 X X X X<br></font><br><br><font face="arial,sans-serif">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.<br>
</font><br>Mat A;<br>int rows = 4;<br>int cols = 8;<br><br>if (mpi_rank() == 0)<br>{<br> MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, rows, cols,<br> 2, PETSC_NULL, 1, PETSC_NULL, &A);<br>
MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);<br><br> static const int R = 2;<br> int set_rows[R] = {0, 1};<br> static const int C = 2;<br> int set_cols[C] = {0, 1};<br> double vals[C*R] = {1., 1., 1., 1.};<br>
MatSetValues(A, R, set_rows, C, set_cols, vals, ADD_VALUES);<br>}<br>else<br>{<br> MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,rows, cols,<br> 2, PETSC_NULL, 1, PETSC_NULL, &A);<br>
MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR);<br><br> static const int R = 2;<br> int set_rows[R] = {2,3};<br> static const int C = 2;<br> int set_cols[C] = {2,3};<br> double vals[C*R] = {1., 1., 1., 1.};<br>
MatSetValues(A, R, set_rows, C, set_cols, vals, ADD_VALUES);<br>}<br><br>MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);<br>MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);<br><br>On Wed, Jul 2, 2008 at 6:35 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
> By definition, the diagonal portion is square. Thus you retrieve the<br>> local rows [start, end), and then the diagonal block is<br>><br>> [start, end) x [start, end)<br>><br>> which you can do with<br>
><br>> <a href="http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatGetOwnershipRange.html">http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatGetOwnershipRange.html</a><br>
><br>> Matt<br>><br>> On Wed, Jul 2, 2008 at 2:52 AM, Andrew Colombi <<a href="mailto:acolombi@gmail.com">acolombi@gmail.com</a>> wrote:<br>>> I'm letting PETSC_DECIDE how to distribute the rows and columns of my<br>
>> non-square matrix, and I want to retrieve the columns it's decided to<br>>> make the "Diagonal" for a particular process. How can I do that? I'm<br>>> probably missing something obvious here, but I checked the docs.<br>
>><br>>> In case some context helps (or inspires a work-around):<br>>><br>>> I want this information to count the number of d_nnz and o_nnz to<br>>> preallocate perfectly. This shouldn't be hard to do, but I need to<br>
>> know the boundary between off-diagonal and diagonal to do it.<br>>><br>>> Thanks,<br>>> -Andrew<br>>><br>>><br>><br>><br>><br>> --<br>> What most experimenters take for granted before they begin their<br>
> experiments is infinitely more interesting than any results to which<br>> their experiments lead.<br>> -- Norbert Wiener<br>><br>><br><br>