Actual source code: ex92.c
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>
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];
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));
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));
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));
40: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
41: PetscCall(PetscRandomSetFromOptions(rand));
43: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
44: PetscCall(MatGetSize(A, &M, &N));
45: Mbs = M / bs;
47: PetscCall(PetscMalloc1(bs, &rows));
48: PetscCall(PetscMalloc1(bs, &cols));
49: PetscCall(PetscMalloc1(bs * bs, &vals));
50: PetscCall(PetscMalloc1(M, &idx));
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: }
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: }
81: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
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: if (flg) {
91: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
92: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric");
93: PetscCall(MatDestroy(&Atrans));
95: /* create a SeqSBAIJ matrix sA (= A) */
96: PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
97: if (vid >= 0 && vid < size) {
98: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "A:\n"));
99: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
100: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "sA:\n"));
101: PetscCall(MatView(sA, PETSC_VIEWER_STDOUT_WORLD));
102: }
104: /* Test sA==A through MatMult() */
105: PetscCall(MatMultEqual(A, sA, 10, &flg));
106: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Error in MatConvert(): A != sA");
108: /* Test MatIncreaseOverlap() */
109: PetscCall(PetscMalloc1(nd, &is1));
110: PetscCall(PetscMalloc1(nd, &is2));
112: for (i = 0; i < nd; i++) {
113: if (!TestAllcols) {
114: PetscCall(PetscRandomGetValue(rand, &rval));
115: sz = (PetscInt)((0.5 + 0.2 * PetscRealPart(rval)) * mbs); /* 0.5*mbs < sz < 0.7*mbs */
117: for (j = 0; j < sz; j++) {
118: PetscCall(PetscRandomGetValue(rand, &rval));
119: idx[j * bs] = bs * (PetscInt)(PetscRealPart(rval) * Mbs);
120: for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
121: }
122: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is1 + i));
123: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sz * bs, idx, PETSC_COPY_VALUES, is2 + i));
124: if (rank == vid) {
125: PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n", rank, i, sz));
126: PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
127: }
128: } else { /* Test all rows and columns */
129: sz = M;
130: PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is1 + i));
131: PetscCall(ISCreateStride(PETSC_COMM_SELF, sz, 0, 1, is2 + i));
133: if (rank == vid) {
134: PetscBool colflag;
135: PetscCall(ISIdentity(is2[i], &colflag));
136: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] is2[%" PetscInt_FMT "], colflag %d\n", rank, i, colflag));
137: PetscCall(ISView(is2[i], PETSC_VIEWER_STDOUT_SELF));
138: }
139: }
140: }
142: PetscCall(PetscLogStageRegister("MatOv_SBAIJ", &stages[0]));
143: PetscCall(PetscLogStageRegister("MatOv_BAIJ", &stages[1]));
145: /* Test MatIncreaseOverlap */
146: if (TestOverlap) {
147: PetscCall(PetscLogStagePush(stages[0]));
148: PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));
149: PetscCall(PetscLogStagePop());
151: PetscCall(PetscLogStagePush(stages[1]));
152: PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
153: PetscCall(PetscLogStagePop());
155: if (rank == vid) {
156: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from BAIJ:\n", rank));
157: PetscCall(ISView(is1[0], PETSC_VIEWER_STDOUT_SELF));
158: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] IS from SBAIJ:\n", rank));
159: PetscCall(ISView(is2[0], PETSC_VIEWER_STDOUT_SELF));
160: }
162: for (i = 0; i < nd; ++i) {
163: PetscCall(ISEqual(is1[i], is2[i], &flg));
164: if (!flg) {
165: if (rank == 0) {
166: PetscCall(ISSort(is1[i]));
167: PetscCall(ISSort(is2[i]));
168: }
169: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i);
170: }
171: }
172: }
174: /* Test MatCreateSubmatrices */
175: if (TestSubMat) {
176: if (test_sorted) {
177: for (i = 0; i < nd; ++i) PetscCall(ISSort(is1[i]));
178: }
179: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
180: PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_INITIAL_MATRIX, &submatsA));
182: PetscCall(MatMultEqual(A, sA, 10, &flg));
183: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A != sA");
185: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
186: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
187: PetscCall(MatCreateSubMatrices(sA, nd, is1, is1, MAT_REUSE_MATRIX, &submatsA));
188: PetscCall(MatMultEqual(A, sA, 10, &flg));
189: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "MatCreateSubmatrices(): A != sA");
191: PetscCall(MatDestroySubMatrices(nd, &submatA));
192: PetscCall(MatDestroySubMatrices(nd, &submatsA));
193: }
195: /* Free allocated memory */
196: for (i = 0; i < nd; ++i) {
197: PetscCall(ISDestroy(&is1[i]));
198: PetscCall(ISDestroy(&is2[i]));
199: }
200: PetscCall(PetscFree(is1));
201: PetscCall(PetscFree(is2));
202: PetscCall(PetscFree(idx));
203: PetscCall(PetscFree(rows));
204: PetscCall(PetscFree(cols));
205: PetscCall(PetscFree(vals));
206: PetscCall(MatDestroy(&A));
207: PetscCall(MatDestroy(&sA));
208: PetscCall(PetscRandomDestroy(&rand));
209: PetscCall(PetscFinalize());
210: return 0;
211: }
213: /*TEST
215: test:
216: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
217: output_file: output/ex92_1.out
219: test:
220: suffix: 2
221: nsize: {{3 4}}
222: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat
223: output_file: output/ex92_1.out
225: test:
226: suffix: 3
227: nsize: {{3 4}}
228: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols
229: output_file: output/ex92_1.out
231: test:
232: suffix: 3_sorted
233: nsize: {{3 4}}
234: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted
235: output_file: output/ex92_1.out
237: test:
238: suffix: 4
239: nsize: {{3 4}}
240: args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols
241: output_file: output/ex92_1.out
243: TEST*/