1 static char help[] = "Extract submatrices using unsorted indices. For SEQSBAIJ either sort both rows and columns, or sort none.\n\n";
2 /*
3 Take a 4x4 grid and form a 5-point stencil graph Laplacian over it.
4 Partition the grid into two subdomains by splitting into two in the j-direction (slowest varying).
5 Impose an overlap of 1 and order the subdomains with the j-direction varying fastest.
6 Extract the subdomain submatrices, one per rank.
7 */
8 /* Results:
9 Sequential:
10 - seqaij: will error out, if rows or columns are unsorted
11 - seqbaij: will extract submatrices correctly even for unsorted row or column indices
12 - seqsbaij: will extract submatrices correctly even for unsorted row and column indices (both must be sorted or not);
13 CANNOT automatically report inversions, because MatGetRow is not available.
14 MPI:
15 - mpiaij: will error out, if columns are unsorted
16 - mpibaij: will error out, if columns are unsorted.
17 - mpisbaij: will error out, if columns are unsorted; even with unsorted rows will produce correct submatrices;
18 CANNOT automatically report inversions, because MatGetRow is not available.
19 */
20
21 #include <petscmat.h>
22 #include <petscis.h>
23
main(int argc,char ** args)24 int main(int argc, char **args)
25 {
26 Mat A, *S;
27 IS rowis[2], colis[2];
28 PetscInt n, N, i, j, k, l, nsub, Jlow[2] = {0, 1}, *jlow, Jhigh[2] = {3, 4}, *jhigh, row, col, *subindices, ncols;
29 const PetscInt *cols;
30 PetscScalar v;
31 PetscMPIInt rank, size, p, inversions, total_inversions;
32 PetscBool sort_rows, sort_cols, show_inversions;
33
34 PetscFunctionBeginUser;
35 PetscCall(PetscInitialize(&argc, &args, NULL, help));
36 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
37 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
38 PetscCheck(size < 3, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "A uniprocessor or two-processor example only.");
39
40 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
41 if (size > 1) {
42 n = 8;
43 N = 16;
44 } else {
45 n = 16;
46 N = 16;
47 }
48 PetscCall(MatSetSizes(A, n, n, N, N));
49 PetscCall(MatSetFromOptions(A));
50 PetscCall(MatSetUp(A));
51
52 /* Don't care if the entries are set multiple times by different procs. */
53 for (i = 0; i < 4; ++i) {
54 for (j = 0; j < 4; ++j) {
55 row = j * 4 + i;
56 v = -1.0;
57 if (i > 0) {
58 col = row - 1;
59 PetscCall(MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES));
60 }
61 if (i < 3) {
62 col = row + 1;
63 PetscCall(MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES));
64 }
65 if (j > 0) {
66 col = row - 4;
67 PetscCall(MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES));
68 }
69 if (j < 3) {
70 col = row + 4;
71 PetscCall(MatSetValues(A, 1, &row, 1, &col, &v, INSERT_VALUES));
72 }
73 v = 4.0;
74 PetscCall(MatSetValues(A, 1, &row, 1, &row, &v, INSERT_VALUES));
75 }
76 }
77 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
78 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
79 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original matrix\n"));
80 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
81
82 if (size > 1) {
83 nsub = 1; /* one subdomain per rank */
84 } else {
85 nsub = 2; /* both subdomains on rank 0 */
86 }
87 if (rank) {
88 jlow = Jlow + 1;
89 jhigh = Jhigh + 1;
90 } else {
91 jlow = Jlow;
92 jhigh = Jhigh;
93 }
94 sort_rows = PETSC_FALSE;
95 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sort_rows", &sort_rows, NULL));
96 sort_cols = PETSC_FALSE;
97 PetscCall(PetscOptionsGetBool(NULL, NULL, "-sort_cols", &sort_cols, NULL));
98 for (l = 0; l < nsub; ++l) {
99 PetscCall(PetscMalloc1(12, &subindices));
100 k = 0;
101 for (i = 0; i < 4; ++i) {
102 for (j = jlow[l]; j < jhigh[l]; ++j) {
103 subindices[k] = j * 4 + i;
104 k++;
105 }
106 }
107 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 12, subindices, PETSC_OWN_POINTER, rowis + l));
108 if ((sort_rows && !sort_cols) || (!sort_rows && sort_cols)) {
109 PetscCall(ISDuplicate(rowis[l], colis + l));
110 } else {
111 PetscCall(PetscObjectReference((PetscObject)rowis[l]));
112 colis[l] = rowis[l];
113 }
114 if (sort_rows) PetscCall(ISSort(rowis[l]));
115 if (sort_cols) PetscCall(ISSort(colis[l]));
116 }
117
118 PetscCall(MatCreateSubMatrices(A, nsub, rowis, colis, MAT_INITIAL_MATRIX, &S));
119
120 show_inversions = PETSC_FALSE;
121
122 PetscCall(PetscOptionsGetBool(NULL, NULL, "-show_inversions", &show_inversions, NULL));
123
124 inversions = 0;
125 for (p = 0; p < size; ++p) {
126 if (p == rank) {
127 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Number of subdomains: %" PetscInt_FMT ":\n", rank, size, nsub));
128 for (l = 0; l < nsub; ++l) {
129 PetscInt i0, i1;
130 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain row IS %" PetscInt_FMT ":\n", rank, size, l));
131 PetscCall(ISView(rowis[l], PETSC_VIEWER_STDOUT_SELF));
132 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Subdomain col IS %" PetscInt_FMT ":\n", rank, size, l));
133 PetscCall(ISView(colis[l], PETSC_VIEWER_STDOUT_SELF));
134 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%" PetscInt_FMT ":%" PetscInt_FMT "]: Submatrix %" PetscInt_FMT ":\n", rank, size, l));
135 PetscCall(MatView(S[l], PETSC_VIEWER_STDOUT_SELF));
136 if (show_inversions) {
137 PetscCall(MatGetOwnershipRange(S[l], &i0, &i1));
138 for (i = i0; i < i1; ++i) {
139 PetscCall(MatGetRow(S[l], i, &ncols, &cols, NULL));
140 for (j = 1; j < ncols; ++j) {
141 if (cols[j] < cols[j - 1]) {
142 PetscCall(PetscPrintf(PETSC_COMM_SELF, "***Inversion in row %" PetscInt_FMT ": col[%" PetscInt_FMT "] = %" PetscInt_FMT " < %" PetscInt_FMT " = col[%" PetscInt_FMT "]\n", i, j, cols[j], cols[j - 1], j - 1));
143 inversions++;
144 }
145 }
146 PetscCall(MatRestoreRow(S[l], i, &ncols, &cols, NULL));
147 }
148 }
149 }
150 }
151 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
152 }
153 if (show_inversions) {
154 PetscCallMPI(MPI_Reduce(&inversions, &total_inversions, 1, MPIU_INT, MPI_SUM, 0, PETSC_COMM_WORLD));
155 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "*Total inversions: %" PetscInt_FMT "\n", total_inversions));
156 }
157 PetscCall(MatDestroy(&A));
158
159 for (l = 0; l < nsub; ++l) {
160 PetscCall(ISDestroy(&rowis[l]));
161 PetscCall(ISDestroy(&colis[l]));
162 }
163 PetscCall(MatDestroySubMatrices(nsub, &S));
164 PetscCall(PetscFinalize());
165 return 0;
166 }
167