1 static char help[] = "Tests MatIncreaseOverlap() and MatCreateSubmatrices() for the parallel case.\n\
2 This example is similar to ex40.c; here the index sets used are random.\n\
3 Input arguments are:\n\
4 -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\
5 -nd <size> : > 0 no of domains per processor \n\
6 -ov <overlap> : >=0 amount of overlap between domains\n\n";
7
8 #include <petscmat.h>
9
main(int argc,char ** args)10 int main(int argc, char **args)
11 {
12 PetscInt nd = 2, ov = 1, i, j, lsize, m, n, *idx, bs;
13 PetscMPIInt rank, size;
14 PetscBool flg;
15 Mat A, B, *submatA, *submatB;
16 char file[PETSC_MAX_PATH_LEN];
17 PetscViewer fd;
18 IS *is1, *is2;
19 PetscRandom r;
20 PetscBool test_unsorted = PETSC_FALSE;
21 PetscScalar rand;
22
23 PetscFunctionBeginUser;
24 PetscCall(PetscInitialize(&argc, &args, NULL, help));
25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
27 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
28 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_unsorted", &test_unsorted, NULL));
31
32 /* Read matrix A and RHS */
33 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
34 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
35 PetscCall(MatSetType(A, MATAIJ));
36 PetscCall(MatSetFromOptions(A));
37 PetscCall(MatLoad(A, fd));
38 PetscCall(PetscViewerDestroy(&fd));
39
40 /* Read the same matrix as a seq matrix B */
41 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, file, FILE_MODE_READ, &fd));
42 PetscCall(MatCreate(PETSC_COMM_SELF, &B));
43 PetscCall(MatSetType(B, MATSEQAIJ));
44 PetscCall(MatSetFromOptions(B));
45 PetscCall(MatLoad(B, fd));
46 PetscCall(PetscViewerDestroy(&fd));
47
48 PetscCall(MatGetBlockSize(A, &bs));
49
50 /* Create the Random no generator */
51 PetscCall(MatGetSize(A, &m, &n));
52 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
53 PetscCall(PetscRandomSetFromOptions(r));
54
55 /* Create the IS corresponding to subdomains */
56 PetscCall(PetscMalloc1(nd, &is1));
57 PetscCall(PetscMalloc1(nd, &is2));
58 PetscCall(PetscMalloc1(m, &idx));
59 for (i = 0; i < m; i++) idx[i] = i;
60
61 /* Create the random Index Sets */
62 for (i = 0; i < nd; i++) {
63 /* Skip a few,so that the IS on different procs are different*/
64 for (j = 0; j < rank; j++) PetscCall(PetscRandomGetValue(r, &rand));
65 PetscCall(PetscRandomGetValue(r, &rand));
66 lsize = (PetscInt)(rand * (m / bs));
67 /* shuffle */
68 for (j = 0; j < lsize; j++) {
69 PetscInt k, swap, l;
70
71 PetscCall(PetscRandomGetValue(r, &rand));
72 k = j + (PetscInt)(rand * ((m / bs) - j));
73 for (l = 0; l < bs; l++) {
74 swap = idx[bs * j + l];
75 idx[bs * j + l] = idx[bs * k + l];
76 idx[bs * k + l] = swap;
77 }
78 }
79 if (!test_unsorted) PetscCall(PetscSortInt(lsize * bs, idx));
80 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i));
81 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i));
82 PetscCall(ISSetBlockSize(is1[i], bs));
83 PetscCall(ISSetBlockSize(is2[i], bs));
84 }
85
86 if (!test_unsorted) {
87 PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
88 PetscCall(MatIncreaseOverlap(B, nd, is2, ov));
89
90 for (i = 0; i < nd; ++i) {
91 PetscCall(ISSort(is1[i]));
92 PetscCall(ISSort(is2[i]));
93 }
94 }
95
96 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
97 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
98
99 /* Now see if the serial and parallel case have the same answers */
100 for (i = 0; i < nd; ++i) {
101 PetscCall(MatEqual(submatA[i], submatB[i], &flg));
102 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%" PetscInt_FMT "-th parallel submatA != seq submatB", i);
103 }
104
105 /* Free Allocated Memory */
106 for (i = 0; i < nd; ++i) {
107 PetscCall(ISDestroy(&is1[i]));
108 PetscCall(ISDestroy(&is2[i]));
109 }
110 PetscCall(MatDestroySubMatrices(nd, &submatA));
111 PetscCall(MatDestroySubMatrices(nd, &submatB));
112
113 PetscCall(PetscRandomDestroy(&r));
114 PetscCall(PetscFree(is1));
115 PetscCall(PetscFree(is2));
116 PetscCall(MatDestroy(&A));
117 PetscCall(MatDestroy(&B));
118 PetscCall(PetscFree(idx));
119 PetscCall(PetscFinalize());
120 return 0;
121 }
122
123 /*TEST
124
125 build:
126 requires: !complex
127
128 test:
129 nsize: 3
130 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
131 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 5 -ov 2
132 output_file: output/empty.out
133
134 test:
135 suffix: 2
136 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -ov 2
137 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
138 output_file: output/empty.out
139
140 test:
141 suffix: unsorted_baij_mpi
142 nsize: 3
143 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
144 args: -f ${DATAFILESPATH}/matrices/cfd.1.10 -nd 8 -mat_type baij -test_unsorted
145 output_file: output/empty.out
146
147 test:
148 suffix: unsorted_baij_seq
149 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
150 args: -f ${DATAFILESPATH}/matrices/cfd.1.10 -nd 8 -mat_type baij -test_unsorted
151 output_file: output/empty.out
152
153 test:
154 suffix: unsorted_mpi
155 nsize: 3
156 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
157 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -test_unsorted
158 output_file: output/empty.out
159
160 test:
161 suffix: unsorted_seq
162 requires: datafilespath double !defined(PETSC_USE_64BIT_INDICES) !complex
163 args: -f ${DATAFILESPATH}/matrices/arco1 -nd 8 -test_unsorted
164 output_file: output/empty.out
165
166 TEST*/
167