1 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel MatSBAIJ format.\n";
2 /* Example of usage:
3 mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat
4 */
5 #include <petscmat.h>
6
main(int argc,char ** args)7 int main(int argc, char **args)
8 {
9 Mat A, Atrans, sA, *submatA, *submatsA;
10 PetscMPIInt size, rank;
11 PetscInt bs = 1, mbs = 10, ov = 1, i, j, k, *rows, *cols, nd = 2, *idx, rstart, rend, sz, M, N, Mbs;
12 PetscScalar *vals, rval, one = 1.0;
13 IS *is1, *is2;
14 PetscRandom rand;
15 PetscBool flg, TestOverlap, TestSubMat, TestAllcols, test_sorted = PETSC_FALSE;
16 PetscInt vid = -1;
17 PetscLogStage stages[2];
18
19 PetscFunctionBeginUser;
20 PetscCall(PetscInitialize(&argc, &args, NULL, help));
21 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23
24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_mbs", &mbs, NULL));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-view_id", &vid, NULL));
29 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_overlap", &TestOverlap));
30 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_submat", &TestSubMat));
31 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_allcols", &TestAllcols));
32 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sorted", &test_sorted, NULL));
33
34 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
35 PetscCall(MatSetSizes(A, mbs * bs, mbs * bs, PETSC_DECIDE, PETSC_DECIDE));
36 PetscCall(MatSetType(A, MATBAIJ));
37 PetscCall(MatSeqBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL));
38 PetscCall(MatMPIBAIJSetPreallocation(A, bs, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
39
40 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
41 PetscCall(PetscRandomSetFromOptions(rand));
42
43 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
44 PetscCall(MatGetSize(A, &M, &N));
45 Mbs = M / bs;
46
47 PetscCall(PetscMalloc1(bs, &rows));
48 PetscCall(PetscMalloc1(bs, &cols));
49 PetscCall(PetscMalloc1(bs * bs, &vals));
50 PetscCall(PetscMalloc1(M, &idx));
51
52 /* Now set blocks of values */
53 for (j = 0; j < bs * bs; j++) vals[j] = 0.0;
54 for (i = 0; i < Mbs; i++) {
55 cols[0] = i * bs;
56 rows[0] = i * bs;
57 for (j = 1; j < bs; j++) {
58 rows[j] = rows[j - 1] + 1;
59 cols[j] = cols[j - 1] + 1;
60 }
61 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
62 }
63 /* second, add random blocks */
64 for (i = 0; i < 20 * bs; i++) {
65 PetscCall(PetscRandomGetValue(rand, &rval));
66 cols[0] = bs * (PetscInt)(PetscRealPart(rval) * Mbs);
67 PetscCall(PetscRandomGetValue(rand, &rval));
68 rows[0] = rstart + bs * (PetscInt)(PetscRealPart(rval) * mbs);
69 for (j = 1; j < bs; j++) {
70 rows[j] = rows[j - 1] + 1;
71 cols[j] = cols[j - 1] + 1;
72 }
73
74 for (j = 0; j < bs * bs; j++) {
75 PetscCall(PetscRandomGetValue(rand, &rval));
76 vals[j] = rval;
77 }
78 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
79 }
80
81 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
82 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
83
84 /* make A a symmetric matrix: A <- A^T + A */
85 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
86 PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN));
87 PetscCall(MatDestroy(&Atrans));
88 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
89 PetscCall(MatEqual(A, Atrans, &flg));
90 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric");
91 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
92 PetscCall(MatDestroy(&Atrans));
93
94 /* create a SeqSBAIJ matrix sA (= A) */
95 PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
96 if (vid >= 0 && vid < size) {
97 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "A:\n"));
98 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
99 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "sA:\n"));
100 PetscCall(MatView(sA, PETSC_VIEWER_STDOUT_WORLD));
101 }
102
103 /* Test sA==A through MatMult() */
104 PetscCall(MatMultEqual(A, sA, 10, &flg));
105 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Error in MatConvert(): A != sA");
106
107 /* Test MatIncreaseOverlap() */
108 PetscCall(PetscMalloc1(nd, &is1));
109 PetscCall(PetscMalloc1(nd, &is2));
110
111 for (i = 0; i < nd; i++) {
112 if (!TestAllcols) {
113 PetscCall(PetscRandomGetValue(rand, &rval));
114 sz = (PetscInt)((0.5 + 0.2 * PetscRealPart(rval)) * mbs); /* 0.5*mbs < sz < 0.7*mbs */
115
116 for (j = 0; j < sz; j++) {
117 PetscCall(PetscRandomGetValue(rand, &rval));
118 idx[j * bs] = bs * (PetscInt)(PetscRealPart(rval) * Mbs);
119 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
120 }
121 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
122 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
123 if (rank == vid) {
124 PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n", rank, i, sz));
125 PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
126 }
127 } else { /* Test all rows and columns */
128 sz = M;
129 PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is1 + i));
130 PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is2 + i));
131
132 if (rank == vid) {
133 PetscBool colflag;
134 PetscCall(ISIdentity(is2[i], &colflag));
135 PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] is2[%" PetscInt_FMT "], colflag %d\n", rank, i, colflag));
136 PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
137 }
138 }
139 }
140
141 PetscCall(PetscLogStageRegister("MatOv_SBAIJ", &stages[0]));
142 PetscCall(PetscLogStageRegister("MatOv_BAIJ", &stages[1]));
143
144 /* Test MatIncreaseOverlap */
145 if (TestOverlap) {
146 PetscCall(PetscLogStagePush(stages[0]));
147 PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));
148 PetscCall(PetscLogStagePop());
149
150 PetscCall(PetscLogStagePush(stages[1]));
151 PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
152 PetscCall(PetscLogStagePop());
153
154 if (rank == vid) {
155 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from BAIJ:\n", rank));
156 PetscCall(ISView(is1[0], PETSC_VIEWER_STDOUT_SELF));
157 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from SBAIJ:\n", rank));
158 PetscCall(ISView(is2[0], PETSC_VIEWER_STDOUT_SELF));
159 }
160
161 for (i = 0; i < nd; ++i) {
162 PetscCall(ISEqual(is1[i], is2[i], &flg));
163 if (!flg) {
164 if (rank == 0) {
165 PetscCall(ISSort(is1[i]));
166 PetscCall(ISSort(is2[i]));
167 }
168 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i);
169 }
170 }
171 }
172
173 /* Test MatCreateSubmatrices */
174 if (TestSubMat) {
175 if (test_sorted) {
176 for (i = 0; i < nd; ++i) PetscCall(ISSort(is1[i]));
177 }
178 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
179 PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_INITIAL_MATRIX, &submatsA));
180
181 PetscCall(MatMultEqual(A, sA, 10, &flg));
182 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A != sA");
183
184 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
185 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
186 PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_REUSE_MATRIX, &submatsA));
187 PetscCall(MatMultEqual(A, sA, 10, &flg));
188 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatCreateSubmatrices(): A != sA");
189
190 PetscCall(MatDestroySubMatrices(nd, &submatA));
191 PetscCall(MatDestroySubMatrices(nd, &submatsA));
192 }
193
194 /* Free allocated memory */
195 for (i = 0; i < nd; ++i) {
196 PetscCall(ISDestroy(&is1[i]));
197 PetscCall(ISDestroy(&is2[i]));
198 }
199 PetscCall(PetscFree(is1));
200 PetscCall(PetscFree(is2));
201 PetscCall(PetscFree(idx));
202 PetscCall(PetscFree(rows));
203 PetscCall(PetscFree(cols));
204 PetscCall(PetscFree(vals));
205 PetscCall(MatDestroy(&A));
206 PetscCall(MatDestroy(&sA));
207 PetscCall(PetscRandomDestroy(&rand));
208 PetscCall(PetscFinalize());
209 return 0;
210 }
211
212 /*TEST
213
214 test:
215 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
216 output_file: output/empty.out
217
218 test:
219 suffix: 2
220 nsize: {{3 4}}
221 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
222 output_file: output/empty.out
223
224 test:
225 suffix: 3
226 nsize: {{3 4}}
227 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols
228 output_file: output/empty.out
229
230 test:
231 suffix: 3_sorted
232 nsize: {{3 4}}
233 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted
234 output_file: output/empty.out
235
236 test:
237 suffix: 4
238 nsize: {{3 4}}
239 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols
240 output_file: output/empty.out
241
242 TEST*/
243