The attached example exposes a bug in some implementations of MatGetSubMatrices():<div><div>The extracted submatrices have row such as these: </div><div><div>row 0: (0, 4) (3, -1) (1, -1) </div><div>row 1: (0, -1) (1, 4) (4, -1) (2, -1) </div>
</div><div>Because of there are inversions among the column indices, binary search will fail for XXXAIJ matrices.</div><div>This is the result of using unsorted column ISs in MatGetSubMatrices().</div><div>According to our documentation, it's legal. Still, some (but not all) implementations (e.g., MPIDENSE, SEQAIJ) </div>
<div>will reject unsorted column indices. </div><div><br></div><div><br></div><div>Note that this subroutine is used by PCASM, where the user's subdomain ISs are </div><div>used to extract the subdomain submatrices. It looks like the problem never got exposed because</div>
<div>by default the subdomain ISs are sorted by ASM internally, unless the user explicitly sets a flag </div><div>(and, incidentally, it doesn't look like the flag can be set on the command line).</div>
<div><br></div><div>The problem seems to exist for the following matrix types:</div><div>MPIAIJ, MPIBAIJ, MPISBAIJ</div><div><br></div><div><div>Unless someone objects I'm going to push a fix to petsc-dev that will check that </div>
<div>all column ISs are sorted in the three implementations above and will throw an error otherwise. </div></div><div><br></div><div>Dmitry.</div><div><br></div><div>PS For those that want to run this test:</div><div><div>
make runsubmat_test NP=2 ARGS="-mat_type mpiaij -sort_rows -show_inversions"</div><div><br></div><div>NP must be 1 if using a sequential -mat_type.</div><div><br></div><div>-show_inversions requires MatGetRow(), which doesn't work with symmetric types.</div>
</div><div><br></div><div>-sort_rows will sort the row ISs, but not the column ISs (that's how one gets past the first line of defense </div><div>for matrix types like MPIDENSE or SEQAIJ, amazingly, which bail as soon as they see unsorted row ISs); </div>
<div>symmetric matrix type will not like this, since they need identical row and column ISs</div><div><br></div>
<div>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):</div><div><div>Original matrix</div><div>Matrix Object: 1 MPI processes</div>
<div> type: seqaij</div><div>row 0: (0, 4) (1, -1) (4, -1) </div><div>row 1: (0, -1) (1, 4) (2, -1) (5, -1) </div><div>row 2: (1, -1) (2, 4) (3, -1) (6, -1) </div><div>row 3: (2, -1) (3, 4) (7, -1) </div><div>
row 4: (0, -1) (4, 4) (5, -1) (8, -1) </div>
<div>row 5: (1, -1) (4, -1) (5, 4) (6, -1) (9, -1) </div><div>row 6: (2, -1) (5, -1) (6, 4) (7, -1) (10, -1) </div><div>row 7: (3, -1) (6, -1) (7, 4) (11, -1) </div><div>row 8: (4, -1) (8, 4) (9, -1) (12, -1) </div>
<div>row 9: (5, -1) (8, -1) (9, 4) (10, -1) (13, -1) </div><div>row 10: (6, -1) (9, -1) (10, 4) (11, -1) (14, -1) </div><div>row 11: (7, -1) (10, -1) (11, 4) (15, -1) </div><div>row 12: (8, -1) (12, 4) (13, -1) </div>
<div>row 13: (9, -1) (12, -1) (13, 4) (14, -1) </div><div>row 14: (10, -1) (13, -1) (14, 4) (15, -1) </div><div>row 15: (11, -1) (14, -1) (15, 4) </div><div>[0:1]: Number of subdomains: 2:</div><div>[0:1]: Subdomain row IS 0:</div>
<div>Number of indices in set 12</div><div>0 0</div><div>1 1</div><div>2 2</div><div>3 3</div><div>4 4</div><div>5 5</div><div>6 6</div><div>7 7</div><div>8 8</div><div>9 9</div><div>10 10</div><div>11 11</div><div>[0:1]: Subdomain col IS 0:</div>
<div>Number of indices in set 12</div><div>0 0</div><div>1 4</div><div>2 8</div><div>3 1</div><div>4 5</div><div>5 9</div><div>6 2</div><div>7 6</div><div>8 10</div><div>9 3</div><div>10 7</div><div>11 11</div><div>[0:1]: Submatrix 0:</div>
<div>Matrix Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 4) (3, -1) (1, -1) </div><div>row 1: (0, -1) (3, 4) (6, -1) (4, -1) </div><div>row 2: (3, -1) (6, 4) (9, -1) (7, -1) </div><div>row 3: (6, -1) (9, 4) (10, -1) </div>
<div>row 4: (0, -1) (1, 4) (4, -1) (2, -1) </div><div>row 5: (3, -1) (1, -1) (4, 4) (7, -1) (5, -1) </div><div>row 6: (6, -1) (4, -1) (7, 4) (10, -1) (8, -1) </div><div>row 7: (9, -1) (7, -1) (10, 4) (11, -1) </div>
<div>row 8: (1, -1) (2, 4) (5, -1) </div><div>row 9: (4, -1) (2, -1) (5, 4) (8, -1) </div><div>row 10: (7, -1) (5, -1) (8, 4) (11, -1) </div><div>row 11: (10, -1) (8, -1) (11, 4) </div><div>***Inversion in row 0: col[2] = 1 < 3 = col[1]</div>
<div>***Inversion in row 1: col[3] = 4 < 6 = col[2]</div><div>***Inversion in row 2: col[3] = 7 < 9 = col[2]</div><div>***Inversion in row 4: col[3] = 2 < 4 = col[2]</div><div>***Inversion in row 5: col[1] = 1 < 3 = col[0]</div>
<div>***Inversion in row 5: col[4] = 5 < 7 = col[3]</div><div>***Inversion in row 6: col[1] = 4 < 6 = col[0]</div><div>***Inversion in row 6: col[4] = 8 < 10 = col[3]</div><div>***Inversion in row 7: col[1] = 7 < 9 = col[0]</div>
<div>***Inversion in row 9: col[1] = 2 < 4 = col[0]</div><div>***Inversion in row 10: col[1] = 5 < 7 = col[0]</div><div>***Inversion in row 11: col[1] = 8 < 10 = col[0]</div><div>[0:1]: Subdomain row IS 1:</div>
<div>
Number of indices in set 12</div><div>0 4</div><div>1 5</div><div>2 6</div><div>3 7</div><div>4 8</div><div>5 9</div><div>6 10</div><div>7 11</div><div>8 12</div><div>9 13</div><div>10 14</div><div>11 15</div><div>[0:1]: Subdomain col IS 1:</div>
<div>Number of indices in set 12</div><div>0 4</div><div>1 8</div><div>2 12</div><div>3 5</div><div>4 9</div><div>5 13</div><div>6 6</div><div>7 10</div><div>8 14</div><div>9 7</div><div>10 11</div><div>11 15</div><div>[0:1]: Submatrix 1:</div>
<div>Matrix Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 4) (3, -1) (1, -1) </div><div>row 1: (0, -1) (3, 4) (6, -1) (4, -1) </div><div>row 2: (3, -1) (6, 4) (9, -1) (7, -1) </div><div>row 3: (6, -1) (9, 4) (10, -1) </div>
<div>row 4: (0, -1) (1, 4) (4, -1) (2, -1) </div><div>row 5: (3, -1) (1, -1) (4, 4) (7, -1) (5, -1) </div><div>row 6: (6, -1) (4, -1) (7, 4) (10, -1) (8, -1) </div><div>row 7: (9, -1) (7, -1) (10, 4) (11, -1) </div>
<div>row 8: (1, -1) (2, 4) (5, -1) </div><div>row 9: (4, -1) (2, -1) (5, 4) (8, -1) </div><div>row 10: (7, -1) (5, -1) (8, 4) (11, -1) </div><div>row 11: (10, -1) (8, -1) (11, 4) </div><div>***Inversion in row 0: col[2] = 1 < 3 = col[1]</div>
<div>***Inversion in row 1: col[3] = 4 < 6 = col[2]</div><div>***Inversion in row 2: col[3] = 7 < 9 = col[2]</div><div>***Inversion in row 4: col[3] = 2 < 4 = col[2]</div><div>***Inversion in row 5: col[1] = 1 < 3 = col[0]</div>
<div>***Inversion in row 5: col[4] = 5 < 7 = col[3]</div><div>***Inversion in row 6: col[1] = 4 < 6 = col[0]</div><div>***Inversion in row 6: col[4] = 8 < 10 = col[3]</div><div>***Inversion in row 7: col[1] = 7 < 9 = col[0]</div>
<div>***Inversion in row 9: col[1] = 2 < 4 = col[0]</div><div>***Inversion in row 10: col[1] = 5 < 7 = col[0]</div><div>***Inversion in row 11: col[1] = 8 < 10 = col[0]</div><div>[0]Total space allocated 16 bytes</div>
<div>[ 0]16 bytes main() line 99 in submat_test.c</div></div><div><br></div>
</div>