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