Actual source code: ex183.c
1: static char help[] = "Example of extracting an array of MPI submatrices from a given MPI matrix.\n"
2: "This test can only be run in parallel.\n"
3: "\n";
5: #include <petscmat.h>
7: PetscErrorCode MyISView(IS *rowis, IS *colis, PetscInt gs, PetscInt ss, PetscViewer viewer)
8: {
9: PetscViewer subviewer = NULL;
11: PetscFunctionBeginUser;
12: PetscCheck(ss <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ss must be less than or equal to 1");
13: PetscCall(PetscViewerASCIIPrintf(viewer, "Row IS %" PetscInt_FMT "\n", gs));
14: PetscCall(PetscViewerGetSubViewer(viewer, ss > -1 ? PetscObjectComm((PetscObject)rowis[ss]) : PETSC_COMM_SELF, &subviewer));
15: if (ss > -1) PetscCall(ISView(rowis[ss], subviewer));
16: PetscCall(PetscViewerRestoreSubViewer(viewer, ss > -1 ? PetscObjectComm((PetscObject)rowis[ss]) : PETSC_COMM_SELF, &subviewer));
17: PetscCall(PetscViewerASCIIPrintf(viewer, "Col IS %" PetscInt_FMT "\n", gs));
18: PetscCall(PetscViewerGetSubViewer(viewer, ss > -1 ? PetscObjectComm((PetscObject)rowis[ss]) : PETSC_COMM_SELF, &subviewer));
19: if (ss > -1) PetscCall(ISView(colis[ss], subviewer));
20: PetscCall(PetscViewerRestoreSubViewer(viewer, ss > -1 ? PetscObjectComm((PetscObject)rowis[ss]) : PETSC_COMM_SELF, &subviewer));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: int main(int argc, char **args)
25: {
26: Mat A, *submats;
27: MPI_Comm subcomm;
28: PetscMPIInt rank, size, subrank, subsize, color;
29: PetscInt m, n, N, bs, rstart, rend, i, j, k, total_subdomains, hash, nsubdomains = 1;
30: PetscInt nis, *cols, gnsubdomains, gsubdomainnums[1], gsubdomainperm[1], s, gs;
31: PetscInt *rowindices, *colindices, idx, rep;
32: PetscScalar *vals;
33: IS rowis[1], colis[1];
34: PetscViewer viewer;
35: PetscBool permute_indices, flg;
36: VecType vtype, subtype;
38: PetscFunctionBeginUser;
39: PetscCall(PetscInitialize(&argc, &args, NULL, help));
40: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
41: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
43: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex183", "Mat");
44: m = 5;
45: PetscCall(PetscOptionsInt("-m", "Local matrix size", "MatSetSizes", m, &m, &flg));
46: total_subdomains = size - 1;
47: PetscCall(PetscOptionsRangeInt("-total_subdomains", "Number of submatrices where 0 < n < comm size", "MatCreateSubMatricesMPI", total_subdomains, &total_subdomains, &flg, 1, size));
48: permute_indices = PETSC_FALSE;
49: PetscCall(PetscOptionsBool("-permute_indices", "Whether to permute indices before breaking them into subdomains", "ISCreateGeneral", permute_indices, &permute_indices, &flg));
50: hash = 7;
51: PetscCall(PetscOptionsInt("-hash", "Permutation factor, which has to be relatively prime to M = size*m (total matrix size)", "ISCreateGeneral", hash, &hash, &flg));
52: rep = 2;
53: PetscCall(PetscOptionsRangeInt("-rep", "Number of times to carry out submatrix extractions; currently only 1 & 2 are supported", NULL, rep, &rep, &flg, 1, 2));
54: PetscOptionsEnd();
56: viewer = PETSC_VIEWER_STDOUT_WORLD;
57: /* Create logically sparse, but effectively dense matrix for easy verification of submatrix extraction correctness. */
58: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
59: PetscCall(MatSetSizes(A, m, m, PETSC_DECIDE, PETSC_DECIDE));
60: PetscCall(MatSetFromOptions(A));
61: PetscCall(MatSetUp(A));
62: PetscCall(MatGetSize(A, NULL, &N));
63: PetscCall(MatGetLocalSize(A, NULL, &n));
64: PetscCall(MatGetBlockSize(A, &bs));
65: PetscCall(MatSeqAIJSetPreallocation(A, n, NULL));
66: PetscCall(MatMPIAIJSetPreallocation(A, n, NULL, N - n, NULL));
67: PetscCall(MatSeqBAIJSetPreallocation(A, bs, n / bs, NULL));
68: PetscCall(MatMPIBAIJSetPreallocation(A, bs, n / bs, NULL, (N - n) / bs, NULL));
69: PetscCall(MatSeqSBAIJSetPreallocation(A, bs, n / bs, NULL));
70: PetscCall(MatMPISBAIJSetPreallocation(A, bs, n / bs, NULL, (N - n) / bs, NULL));
72: PetscCall(PetscMalloc2(N, &cols, N, &vals));
73: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
74: for (j = 0; j < N; ++j) cols[j] = j;
75: for (i = rstart; i < rend; i++) {
76: for (j = 0; j < N; ++j) vals[j] = i * 10000 + j;
77: PetscCall(MatSetValues(A, 1, &i, N, cols, vals, INSERT_VALUES));
78: }
79: PetscCall(PetscFree2(cols, vals));
80: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
81: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
83: PetscCall(PetscViewerASCIIPrintf(viewer, "Initial matrix:\n"));
84: PetscCall(MatView(A, viewer));
86: /*
87: Create subcomms and ISs so that each rank participates in one IS.
88: The IS either coalesces adjacent rank indices (contiguous),
89: or selects indices by scrambling them using a hash.
90: */
91: k = size / total_subdomains + (size % total_subdomains > 0); /* There are up to k ranks to a color */
92: color = rank / k;
93: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, rank, &subcomm));
94: PetscCallMPI(MPI_Comm_size(subcomm, &subsize));
95: PetscCallMPI(MPI_Comm_rank(subcomm, &subrank));
96: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
97: nis = 1;
98: PetscCall(PetscMalloc2(rend - rstart, &rowindices, rend - rstart, &colindices));
100: for (j = rstart; j < rend; ++j) {
101: if (permute_indices) {
102: idx = (j / bs * hash) * bs + j % bs;
103: } else {
104: idx = j;
105: }
106: rowindices[j - rstart] = idx % N;
107: colindices[j - rstart] = (idx + m) % N;
108: }
109: PetscCall(ISCreateGeneral(subcomm, rend - rstart, rowindices, PETSC_COPY_VALUES, &rowis[0]));
110: PetscCall(ISCreateGeneral(subcomm, rend - rstart, colindices, PETSC_COPY_VALUES, &colis[0]));
111: PetscCall(ISSetBlockSize(rowis[0], bs));
112: PetscCall(ISSetBlockSize(colis[0], bs));
113: PetscCall(ISSort(rowis[0]));
114: PetscCall(ISSort(colis[0]));
115: PetscCall(PetscFree2(rowindices, colindices));
116: /*
117: Now view the ISs. To avoid deadlock when viewing a list of objects on different subcomms,
118: we need to obtain the global numbers of our local objects and wait for the corresponding global
119: number to be viewed.
120: */
121: PetscCall(PetscViewerASCIIPrintf(viewer, "Subdomains"));
122: if (permute_indices) PetscCall(PetscViewerASCIIPrintf(viewer, " (hash=%" PetscInt_FMT ")", hash));
123: PetscCall(PetscViewerASCIIPrintf(viewer, ":\n"));
124: PetscCall(PetscViewerFlush(viewer));
126: nsubdomains = 1;
127: for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
128: PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD, 1, (PetscObject *)rowis, &gnsubdomains, gsubdomainnums));
129: PetscCall(PetscSortIntWithPermutation(nsubdomains, gsubdomainnums, gsubdomainperm));
130: for (gs = 0, s = 0; gs < gnsubdomains; ++gs) {
131: PetscInt ss;
132: if (s < nsubdomains) {
133: ss = gsubdomainperm[s];
134: if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
135: ++s;
136: } else ss = -1;
137: } else ss = -1;
138: PetscCall(MyISView(rowis, colis, gs, ss, viewer));
139: }
140: PetscCall(PetscViewerFlush(viewer));
141: PetscCall(ISSort(rowis[0]));
142: PetscCall(ISSort(colis[0]));
143: nsubdomains = 1;
144: PetscCall(MatCreateSubMatricesMPI(A, nsubdomains, rowis, colis, MAT_INITIAL_MATRIX, &submats));
145: PetscCall(MatGetVecType(A, &vtype));
146: PetscCall(MatGetVecType(submats[0], &subtype));
147: PetscCall(PetscStrcmp(vtype, subtype, &flg));
148: PetscCheck(flg, subcomm, PETSC_ERR_PLIB, "Submatrix vector type %s differs from the parent vector type %s", subtype, vtype);
149: PetscCall(MatGetBlockSizes(submats[0], &i, &j));
150: PetscCheck(i == bs && j == bs, subcomm, PETSC_ERR_PLIB, "Submatrix block sizes %" PetscInt_FMT " and %" PetscInt_FMT " differ from the requested block size %" PetscInt_FMT, i, j, bs);
151: /*
152: Now view the matrices. To avoid deadlock when viewing a list of objects on different subcomms,
153: we need to obtain the global numbers of our local objects and wait for the corresponding global
154: number to be viewed. Also all MPI processes need to call PetscViewerGetSubViewer() the same number of times
155: */
156: PetscCall(PetscViewerASCIIPrintf(viewer, "Submatrices (repetition 1):\n"));
157: for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
158: PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD, 1, (PetscObject *)submats, &gnsubdomains, gsubdomainnums));
159: PetscCall(PetscSortIntWithPermutation(nsubdomains, gsubdomainnums, gsubdomainperm));
160: for (gs = 0, s = 0; gs < gnsubdomains; ++gs) {
161: PetscViewer subviewer = NULL;
162: if (s < nsubdomains) {
163: PetscInt ss;
164: ss = gsubdomainperm[s];
165: if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
166: PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)submats[ss]), &subviewer));
167: PetscCall(MatView(submats[ss], subviewer));
168: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)submats[ss]), &subviewer));
169: ++s;
170: } else {
171: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
172: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
173: }
174: } else {
175: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
176: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
177: }
178: }
179: if (rep == 1) goto cleanup;
180: nsubdomains = 1;
181: PetscCall(MatCreateSubMatricesMPI(A, nsubdomains, rowis, colis, MAT_REUSE_MATRIX, &submats));
182: PetscCall(MatGetVecType(A, &vtype));
183: PetscCall(MatGetVecType(submats[0], &subtype));
184: PetscCall(PetscStrcmp(vtype, subtype, &flg));
185: PetscCheck(flg, subcomm, PETSC_ERR_PLIB, "Submatrix vector type %s differs from the parent vector type %s", subtype, vtype);
186: PetscCall(MatGetBlockSizes(submats[0], &i, &j));
187: PetscCheck(i == bs && j == bs, subcomm, PETSC_ERR_PLIB, "Submatrix block sizes %" PetscInt_FMT " and %" PetscInt_FMT " differ from the requested block size %" PetscInt_FMT, i, j, bs);
188: /*
189: Now view the matrices. To avoid deadlock when viewing a list of objects on different subcomms,
190: we need to obtain the global numbers of our local objects and wait for the corresponding global
191: number to be viewed.
192: */
193: PetscCall(PetscViewerASCIIPrintf(viewer, "Submatrices (repetition 2):\n"));
194: for (s = 0; s < nsubdomains; ++s) gsubdomainperm[s] = s;
195: PetscCall(PetscObjectsListGetGlobalNumbering(PETSC_COMM_WORLD, 1, (PetscObject *)submats, &gnsubdomains, gsubdomainnums));
196: PetscCall(PetscSortIntWithPermutation(nsubdomains, gsubdomainnums, gsubdomainperm));
197: for (gs = 0, s = 0; gs < gnsubdomains; ++gs) {
198: PetscViewer subviewer = NULL;
199: if (s < nsubdomains) {
200: PetscInt ss;
201: ss = gsubdomainperm[s];
202: if (gs == gsubdomainnums[ss]) { /* Global subdomain gs being viewed is my subdomain with local number ss. */
203: PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)submats[ss]), &subviewer));
204: PetscCall(MatView(submats[ss], subviewer));
205: PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)submats[ss]), &subviewer));
206: ++s;
207: } else {
208: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
209: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
210: }
211: } else {
212: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
213: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &subviewer));
214: }
215: }
216: cleanup:
217: for (k = 0; k < nsubdomains; ++k) PetscCall(MatDestroy(submats + k));
218: PetscCall(PetscFree(submats));
219: for (k = 0; k < nis; ++k) {
220: PetscCall(ISDestroy(rowis + k));
221: PetscCall(ISDestroy(colis + k));
222: }
223: PetscCall(MatDestroy(&A));
224: PetscCallMPI(MPI_Comm_free(&subcomm));
225: PetscCall(PetscFinalize());
226: return 0;
227: }
229: /*TEST
231: test:
232: nsize: 2
233: args: -total_subdomains 1
234: output_file: output/ex183_2_1.out
236: test:
237: suffix: 2
238: nsize: 3
239: args: -total_subdomains 2
240: output_file: output/ex183_3_2.out
242: test:
243: suffix: 3
244: nsize: 4
245: args: -total_subdomains 2
246: output_file: output/ex183_4_2.out
248: test:
249: suffix: 4
250: nsize: 6
251: args: -total_subdomains 2
252: output_file: output/ex183_6_2.out
254: test:
255: suffix: block_size
256: nsize: 4
257: args: -total_subdomains 2 -m 6 -mat_block_size {{2 3}} -permute_indices {{false true}}
258: filter: grep Submatrices
259: output_file: output/ex183_block_size.out
261: test:
262: suffix: cuda
263: nsize: 3
264: requires: cuda
265: filter: sed 's/aijcusparse/aij/g'
266: args: -total_subdomains 2 -mat_type aijcusparse
267: output_file: output/ex183_3_2.out
269: TEST*/