[petsc-dev] Bug in MatGetSubMatrices_MPIXAIJ()

[Dmitry Karpeev]

The attached example exposes a bug in some implementations of
MatGetSubMatrices():
The  extracted submatrices have row such as these:
row 0: (0, 4)  (3, -1)  (1, -1)
row 1: (0, -1)  (1, 4)  (4, -1)  (2, -1)
Because of there are inversions among the column indices, binary search
will fail for XXXAIJ matrices.
This is the result of using unsorted column ISs in MatGetSubMatrices().
According to our documentation, it's legal.  Still, some (but not all)
implementations (e.g., MPIDENSE, SEQAIJ)
will reject unsorted column indices.


Note that this subroutine is used by PCASM, where the user's subdomain ISs
are
used to extract the subdomain submatrices.  It looks like the problem never
got exposed because
by default the subdomain ISs are sorted by ASM internally, unless the user
explicitly sets a flag
(and, incidentally, it doesn't look like the flag can be set on the command
line).

The problem seems to exist for the following matrix types:
MPIAIJ, MPIBAIJ, MPISBAIJ

Unless someone objects I'm going to push a fix to petsc-dev that will check
that
all column ISs are sorted in the three implementations above and will throw
an error otherwise.

Dmitry.

PS For those that want to run this test:
make runsubmat_test NP=2 ARGS="-mat_type mpiaij -sort_rows -show_inversions"

NP must be 1 if using a sequential -mat_type.

-show_inversions requires MatGetRow(), which doesn't work with symmetric
types.

-sort_rows will sort the row ISs, but not the column ISs (that's how one
gets past the first line of defense
for matrix types like MPIDENSE or SEQAIJ, amazingly, which bail as soon as
they see unsorted row ISs);
symmetric matrix type will not like this, since they need identical row and
column ISs

Here's the output of make runsubmat_test NP=2 ARGS="-mat_type mpiaij
-sort_rows -show_inversions" (I can't figure out what's causing the memory
leak):
Original matrix
Matrix Object: 1 MPI processes
  type: seqaij
row 0: (0, 4)  (1, -1)  (4, -1)
row 1: (0, -1)  (1, 4)  (2, -1)  (5, -1)
row 2: (1, -1)  (2, 4)  (3, -1)  (6, -1)
row 3: (2, -1)  (3, 4)  (7, -1)
row 4: (0, -1)  (4, 4)  (5, -1)  (8, -1)
row 5: (1, -1)  (4, -1)  (5, 4)  (6, -1)  (9, -1)
row 6: (2, -1)  (5, -1)  (6, 4)  (7, -1)  (10, -1)
row 7: (3, -1)  (6, -1)  (7, 4)  (11, -1)
row 8: (4, -1)  (8, 4)  (9, -1)  (12, -1)
row 9: (5, -1)  (8, -1)  (9, 4)  (10, -1)  (13, -1)
row 10: (6, -1)  (9, -1)  (10, 4)  (11, -1)  (14, -1)
row 11: (7, -1)  (10, -1)  (11, 4)  (15, -1)
row 12: (8, -1)  (12, 4)  (13, -1)
row 13: (9, -1)  (12, -1)  (13, 4)  (14, -1)
row 14: (10, -1)  (13, -1)  (14, 4)  (15, -1)
row 15: (11, -1)  (14, -1)  (15, 4)
[0:1]: Number of subdomains: 2:
[0:1]: Subdomain row IS 0:
Number of indices in set 12
0 0
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
[0:1]: Subdomain col IS 0:
Number of indices in set 12
0 0
1 4
2 8
3 1
4 5
5 9
6 2
7 6
8 10
9 3
10 7
11 11
[0:1]: Submatrix 0:
Matrix Object: 1 MPI processes
  type: seqaij
row 0: (0, 4)  (3, -1)  (1, -1)
row 1: (0, -1)  (3, 4)  (6, -1)  (4, -1)
row 2: (3, -1)  (6, 4)  (9, -1)  (7, -1)
row 3: (6, -1)  (9, 4)  (10, -1)
row 4: (0, -1)  (1, 4)  (4, -1)  (2, -1)
row 5: (3, -1)  (1, -1)  (4, 4)  (7, -1)  (5, -1)
row 6: (6, -1)  (4, -1)  (7, 4)  (10, -1)  (8, -1)
row 7: (9, -1)  (7, -1)  (10, 4)  (11, -1)
row 8: (1, -1)  (2, 4)  (5, -1)
row 9: (4, -1)  (2, -1)  (5, 4)  (8, -1)
row 10: (7, -1)  (5, -1)  (8, 4)  (11, -1)
row 11: (10, -1)  (8, -1)  (11, 4)
***Inversion in row 0: col[2] = 1 < 3 = col[1]
***Inversion in row 1: col[3] = 4 < 6 = col[2]
***Inversion in row 2: col[3] = 7 < 9 = col[2]
***Inversion in row 4: col[3] = 2 < 4 = col[2]
***Inversion in row 5: col[1] = 1 < 3 = col[0]
***Inversion in row 5: col[4] = 5 < 7 = col[3]
***Inversion in row 6: col[1] = 4 < 6 = col[0]
***Inversion in row 6: col[4] = 8 < 10 = col[3]
***Inversion in row 7: col[1] = 7 < 9 = col[0]
***Inversion in row 9: col[1] = 2 < 4 = col[0]
***Inversion in row 10: col[1] = 5 < 7 = col[0]
***Inversion in row 11: col[1] = 8 < 10 = col[0]
[0:1]: Subdomain row IS 1:
Number of indices in set 12
0 4
1 5
2 6
3 7
4 8
5 9
6 10
7 11
8 12
9 13
10 14
11 15
[0:1]: Subdomain col IS 1:
Number of indices in set 12
0 4
1 8
2 12
3 5
4 9
5 13
6 6
7 10
8 14
9 7
10 11
11 15
[0:1]: Submatrix 1:
Matrix Object: 1 MPI processes
  type: seqaij
row 0: (0, 4)  (3, -1)  (1, -1)
row 1: (0, -1)  (3, 4)  (6, -1)  (4, -1)
row 2: (3, -1)  (6, 4)  (9, -1)  (7, -1)
row 3: (6, -1)  (9, 4)  (10, -1)
row 4: (0, -1)  (1, 4)  (4, -1)  (2, -1)
row 5: (3, -1)  (1, -1)  (4, 4)  (7, -1)  (5, -1)
row 6: (6, -1)  (4, -1)  (7, 4)  (10, -1)  (8, -1)
row 7: (9, -1)  (7, -1)  (10, 4)  (11, -1)
row 8: (1, -1)  (2, 4)  (5, -1)
row 9: (4, -1)  (2, -1)  (5, 4)  (8, -1)
row 10: (7, -1)  (5, -1)  (8, 4)  (11, -1)
row 11: (10, -1)  (8, -1)  (11, 4)
***Inversion in row 0: col[2] = 1 < 3 = col[1]
***Inversion in row 1: col[3] = 4 < 6 = col[2]
***Inversion in row 2: col[3] = 7 < 9 = col[2]
***Inversion in row 4: col[3] = 2 < 4 = col[2]
***Inversion in row 5: col[1] = 1 < 3 = col[0]
***Inversion in row 5: col[4] = 5 < 7 = col[3]
***Inversion in row 6: col[1] = 4 < 6 = col[0]
***Inversion in row 6: col[4] = 8 < 10 = col[3]
***Inversion in row 7: col[1] = 7 < 9 = col[0]
***Inversion in row 9: col[1] = 2 < 4 = col[0]
***Inversion in row 10: col[1] = 5 < 7 = col[0]
***Inversion in row 11: col[1] = 8 < 10 = col[0]
[0]Total space allocated 16 bytes
[ 0]16 bytes main() line 99 in submat_test.c
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120523/3d339bfe/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: submat_test.c
Type: text/x-csrc
Size: 5143 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120523/3d339bfe/attachment.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: makefile
Type: application/octet-stream
Size: 462 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120523/3d339bfe/attachment.obj>