Actual source code: iscoloring.c
1: #include <petsc/private/isimpl.h>
2: #include <petscviewer.h>
3: #include <petscsf.h>
5: const char *const ISColoringTypes[] = {"global", "ghosted", "ISColoringType", "IS_COLORING_", NULL};
7: PetscErrorCode ISColoringReference(ISColoring coloring)
8: {
9: PetscFunctionBegin;
10: coloring->refct++;
11: PetscFunctionReturn(PETSC_SUCCESS);
12: }
14: /*@
15: ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation of a `Mat`
17: Collective
19: Input Parameters:
20: + coloring - the coloring object
21: - type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`
23: Level: intermediate
25: Notes:
26: `IS_COLORING_LOCAL` can lead to faster computations since parallel ghost point updates are not needed for each color
28: With `IS_COLORING_LOCAL` the coloring is in the numbering of the local vector, for `IS_COLORING_GLOBAL` it is in the numbering of the global vector
30: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringGetType()`
31: @*/
32: PetscErrorCode ISColoringSetType(ISColoring coloring, ISColoringType type)
33: {
34: PetscFunctionBegin;
35: coloring->ctype = type;
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: /*@
40: ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation
42: Collective
44: Input Parameter:
45: . coloring - the coloring object
47: Output Parameter:
48: . type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`
50: Level: intermediate
52: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringSetType()`
53: @*/
54: PetscErrorCode ISColoringGetType(ISColoring coloring, ISColoringType *type)
55: {
56: PetscFunctionBegin;
57: *type = coloring->ctype;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: /*@
62: ISColoringDestroy - Destroys an `ISColoring` coloring context.
64: Collective
66: Input Parameter:
67: . iscoloring - the coloring context
69: Level: advanced
71: .seealso: `ISColoring`, `ISColoringView()`, `MatColoring`
72: @*/
73: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
74: {
75: PetscInt i;
77: PetscFunctionBegin;
78: if (!*iscoloring) PetscFunctionReturn(PETSC_SUCCESS);
79: PetscAssertPointer(*iscoloring, 1);
80: if (--(*iscoloring)->refct > 0) {
81: *iscoloring = NULL;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: if ((*iscoloring)->is) {
86: for (i = 0; i < (*iscoloring)->n; i++) PetscCall(ISDestroy(&(*iscoloring)->is[i]));
87: PetscCall(PetscFree((*iscoloring)->is));
88: }
89: if ((*iscoloring)->allocated) PetscCall(PetscFree((*iscoloring)->colors));
90: PetscCall(PetscCommDestroy(&(*iscoloring)->comm));
91: PetscCall(PetscFree(*iscoloring));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@
96: ISColoringViewFromOptions - Processes command line options to determine if/how an `ISColoring` object is to be viewed.
98: Collective
100: Input Parameters:
101: + obj - the `ISColoring` object
102: . bobj - prefix to use for viewing, or `NULL` to use prefix of `mat`
103: - name - option to activate viewing
105: Options Database Key:
106: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
108: Level: intermediate
110: Developer Note:
111: This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`
113: .seealso: `ISColoring`, `ISColoringView()`, `PetscObjectViewFromOptions()`
114: @*/
115: PetscErrorCode ISColoringViewFromOptions(ISColoring obj, PetscObject bobj, const char name[])
116: {
117: PetscViewer viewer;
118: PetscBool flg;
119: PetscViewerFormat format;
120: char *prefix;
122: PetscFunctionBegin;
123: prefix = bobj ? bobj->prefix : NULL;
124: PetscCall(PetscOptionsCreateViewer(obj->comm, NULL, prefix, name, &viewer, &format, &flg));
125: if (flg) {
126: PetscCall(PetscViewerPushFormat(viewer, format));
127: PetscCall(ISColoringView(obj, viewer));
128: PetscCall(PetscViewerPopFormat(viewer));
129: PetscCall(PetscViewerDestroy(&viewer));
130: }
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: ISColoringView - Views an `ISColoring` coloring context.
137: Collective
139: Input Parameters:
140: + iscoloring - the coloring context
141: - viewer - the viewer
143: Level: advanced
145: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
146: @*/
147: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
148: {
149: PetscInt i;
150: PetscBool isascii;
151: IS *is;
153: PetscFunctionBegin;
154: PetscAssertPointer(iscoloring, 1);
155: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(iscoloring->comm, &viewer));
158: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
159: if (isascii) {
160: MPI_Comm comm;
161: PetscMPIInt size, rank;
163: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
164: PetscCallMPI(MPI_Comm_size(comm, &size));
165: PetscCallMPI(MPI_Comm_rank(comm, &rank));
166: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size));
167: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]));
168: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
169: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n));
170: PetscCall(PetscViewerFlush(viewer));
171: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
172: }
174: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is));
175: for (i = 0; i < iscoloring->n; i++) PetscCall(ISView(iscoloring->is[i], viewer));
176: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@C
181: ISColoringGetColors - Returns an array with the color for each local node
183: Not Collective
185: Input Parameter:
186: . iscoloring - the coloring context
188: Output Parameters:
189: + n - number of nodes
190: . nc - number of colors
191: - colors - color for each node
193: Level: advanced
195: Notes:
196: Do not free the `colors` array.
198: The `colors` array will only be valid for the lifetime of the `ISColoring`
200: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
201: @*/
202: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
203: {
204: PetscFunctionBegin;
205: PetscAssertPointer(iscoloring, 1);
207: if (n) *n = iscoloring->N;
208: if (nc) *nc = iscoloring->n;
209: if (colors) *colors = iscoloring->colors;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@C
214: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
216: Collective
218: Input Parameters:
219: + iscoloring - the coloring context
220: - mode - if this value is `PETSC_OWN_POINTER` then the caller owns the pointer and must free the array of `IS` and each `IS` in the array
222: Output Parameters:
223: + nn - number of index sets in the coloring context
224: - isis - array of index sets
226: Level: advanced
228: Note:
229: If mode is `PETSC_USE_POINTER` then `ISColoringRestoreIS()` must be called when the `IS` are no longer needed
231: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetColoring()`, `ISColoringGetColors()`
232: @*/
233: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
234: {
235: PetscFunctionBegin;
236: PetscAssertPointer(iscoloring, 1);
238: if (nn) *nn = iscoloring->n;
239: if (isis) {
240: if (!iscoloring->is) {
241: PetscInt *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
242: ISColoringValue *colors = iscoloring->colors;
243: IS *is;
245: if (PetscDefined(USE_DEBUG)) {
246: for (i = 0; i < n; i++) PetscCheck(((PetscInt)colors[i]) < nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coloring is our of range index %" PetscInt_FMT "value %d number colors %" PetscInt_FMT, i, (int)colors[i], nc);
247: }
249: /* generate the lists of nodes for each color */
250: PetscCall(PetscCalloc1(nc, &mcolors));
251: for (i = 0; i < n; i++) mcolors[colors[i]]++;
253: PetscCall(PetscMalloc1(nc, &ii));
254: PetscCall(PetscMalloc1(n, &ii[0]));
255: for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
256: PetscCall(PetscArrayzero(mcolors, nc));
258: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
259: PetscCallMPI(MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm));
260: base -= iscoloring->N;
261: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
262: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
263: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
264: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");
266: PetscCall(PetscMalloc1(nc, &is));
267: for (i = 0; i < nc; i++) PetscCall(ISCreateGeneral(iscoloring->comm, mcolors[i], ii[i], PETSC_COPY_VALUES, is + i));
269: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
270: *isis = is;
271: PetscCall(PetscFree(ii[0]));
272: PetscCall(PetscFree(ii));
273: PetscCall(PetscFree(mcolors));
274: } else {
275: *isis = iscoloring->is;
276: if (mode == PETSC_OWN_POINTER) iscoloring->is = NULL;
277: }
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: /*@C
283: ISColoringRestoreIS - Restores the index sets extracted from the coloring context with `ISColoringGetIS()` using `PETSC_USE_POINTER`
285: Collective
287: Input Parameters:
288: + iscoloring - the coloring context
289: . mode - who retains ownership of the is
290: - is - array of index sets
292: Level: advanced
294: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
295: @*/
296: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
297: {
298: PetscFunctionBegin;
299: PetscAssertPointer(iscoloring, 1);
301: /* currently nothing is done here */
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*@
306: ISColoringCreate - Generates an `ISColoring` context from lists (provided by each MPI process) of colors for each node.
308: Collective
310: Input Parameters:
311: + comm - communicator for the processors creating the coloring
312: . ncolors - max color value
313: . n - number of nodes on this processor
314: . colors - array containing the colors for this MPI rank, color numbers begin at 0, for each local node
315: - mode - see `PetscCopyMode` for meaning of this flag.
317: Output Parameter:
318: . iscoloring - the resulting coloring data structure
320: Options Database Key:
321: . -is_coloring_view - Activates `ISColoringView()`
323: Level: advanced
325: Notes:
326: By default sets coloring type to `IS_COLORING_GLOBAL`
328: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
329: @*/
330: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
331: {
332: PetscMPIInt size, rank, tag;
333: PetscInt base, top, i;
334: PetscInt nc, ncwork;
335: MPI_Status status;
337: PetscFunctionBegin;
338: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
339: PetscCheck(ncolors <= PETSC_UINT16_MAX, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds %d limit. This number is unrealistic. Perhaps a bug in code? Current max: %d user requested: %" PetscInt_FMT, PETSC_UINT16_MAX, PETSC_IS_COLORING_MAX, ncolors);
340: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short? Current max: %d user requested: %" PetscInt_FMT, PETSC_IS_COLORING_MAX, ncolors);
341: }
342: PetscCall(PetscNew(iscoloring));
343: PetscCall(PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag));
344: comm = (*iscoloring)->comm;
346: /* compute the number of the first node on my processor */
347: PetscCallMPI(MPI_Comm_size(comm, &size));
349: /* should use MPI_Scan() */
350: PetscCallMPI(MPI_Comm_rank(comm, &rank));
351: if (rank == 0) {
352: base = 0;
353: top = n;
354: } else {
355: PetscCallMPI(MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status));
356: top = base + n;
357: }
358: if (rank < size - 1) PetscCallMPI(MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm));
360: /* compute the total number of colors */
361: ncwork = 0;
362: for (i = 0; i < n; i++) {
363: if (ncwork < colors[i]) ncwork = colors[i];
364: }
365: ncwork++;
366: PetscCallMPI(MPIU_Allreduce(&ncwork, &nc, 1, MPIU_INT, MPI_MAX, comm));
367: PetscCheck(nc <= ncolors, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of colors passed in %" PetscInt_FMT " is less than the actual number of colors in array %" PetscInt_FMT, ncolors, nc);
368: (*iscoloring)->n = nc;
369: (*iscoloring)->is = NULL;
370: (*iscoloring)->N = n;
371: (*iscoloring)->refct = 1;
372: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
373: if (mode == PETSC_COPY_VALUES) {
374: PetscCall(PetscMalloc1(n, &(*iscoloring)->colors));
375: PetscCall(PetscArraycpy((*iscoloring)->colors, colors, n));
376: (*iscoloring)->allocated = PETSC_TRUE;
377: } else if (mode == PETSC_OWN_POINTER) {
378: (*iscoloring)->colors = (ISColoringValue *)colors;
379: (*iscoloring)->allocated = PETSC_TRUE;
380: } else {
381: (*iscoloring)->colors = (ISColoringValue *)colors;
382: (*iscoloring)->allocated = PETSC_FALSE;
383: }
384: PetscCall(ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view"));
385: PetscCall(PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: /*@
390: ISBuildTwoSided - Takes an `IS` that describes where each element will be mapped globally over all ranks.
391: Generates an `IS` that contains new numbers from remote or local on the `IS`.
393: Collective
395: Input Parameters:
396: + ito - an `IS` describes to which rank each entry will be mapped. Negative target rank will be ignored
397: - toindx - an `IS` describes what indices should send. `NULL` means sending natural numbering
399: Output Parameter:
400: . rows - contains new numbers from remote or local
402: Level: advanced
404: Developer Note:
405: This manual page is incomprehensible and still needs to be fixed
407: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
408: @*/
409: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
410: {
411: const PetscInt *ito_indices, *toindx_indices;
412: PetscInt *send_indices, rstart, *recv_indices, nrecvs, nsends;
413: PetscInt *tosizes, *fromsizes, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
414: PetscMPIInt *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
415: PetscLayout isrmap;
416: MPI_Comm comm;
417: PetscSF sf;
418: PetscSFNode *iremote;
420: PetscFunctionBegin;
421: PetscCall(PetscObjectGetComm((PetscObject)ito, &comm));
422: PetscCallMPI(MPI_Comm_size(comm, &size));
423: PetscCall(ISGetLocalSize(ito, &ito_ln));
424: PetscCall(ISGetLayout(ito, &isrmap));
425: PetscCall(PetscLayoutGetRange(isrmap, &rstart, NULL));
426: PetscCall(ISGetIndices(ito, &ito_indices));
427: PetscCall(PetscCalloc2(size, &tosizes_tmp, size + 1, &tooffsets_tmp));
428: for (PetscInt i = 0; i < ito_ln; i++) {
429: if (ito_indices[i] < 0) continue;
430: else PetscCheck(ito_indices[i] < size, comm, PETSC_ERR_ARG_OUTOFRANGE, "target rank %" PetscInt_FMT " is larger than communicator size %d ", ito_indices[i], size);
431: tosizes_tmp[ito_indices[i]]++;
432: }
433: nto = 0;
434: for (PetscMPIInt i = 0; i < size; i++) {
435: tooffsets_tmp[i + 1] = tooffsets_tmp[i] + tosizes_tmp[i];
436: if (tosizes_tmp[i] > 0) nto++;
437: }
438: PetscCall(PetscCalloc2(nto, &toranks, 2 * nto, &tosizes));
439: nto = 0;
440: for (PetscMPIInt i = 0; i < size; i++) {
441: if (tosizes_tmp[i] > 0) {
442: toranks[nto] = i;
443: tosizes[2 * nto] = tosizes_tmp[i]; /* size */
444: tosizes[2 * nto + 1] = tooffsets_tmp[i]; /* offset */
445: nto++;
446: }
447: }
448: nsends = tooffsets_tmp[size];
449: PetscCall(PetscCalloc1(nsends, &send_indices));
450: if (toindx) PetscCall(ISGetIndices(toindx, &toindx_indices));
451: for (PetscInt i = 0; i < ito_ln; i++) {
452: if (ito_indices[i] < 0) continue;
453: PetscCall(PetscMPIIntCast(ito_indices[i], &target_rank));
454: send_indices[tooffsets_tmp[target_rank]] = toindx ? toindx_indices[i] : (i + rstart);
455: tooffsets_tmp[target_rank]++;
456: }
457: if (toindx) PetscCall(ISRestoreIndices(toindx, &toindx_indices));
458: PetscCall(ISRestoreIndices(ito, &ito_indices));
459: PetscCall(PetscFree2(tosizes_tmp, tooffsets_tmp));
460: PetscCall(PetscCommBuildTwoSided(comm, 2, MPIU_INT, nto, toranks, tosizes, &nfrom, &fromranks, &fromsizes));
461: PetscCall(PetscFree2(toranks, tosizes));
462: PetscCall(PetscMalloc1(nfrom, &fromperm_newtoold));
463: for (PetscMPIInt i = 0; i < nfrom; i++) fromperm_newtoold[i] = i;
464: PetscCall(PetscSortMPIIntWithArray(nfrom, fromranks, fromperm_newtoold));
465: nrecvs = 0;
466: for (PetscMPIInt i = 0; i < nfrom; i++) nrecvs += fromsizes[i * 2];
467: PetscCall(PetscCalloc1(nrecvs, &recv_indices));
468: PetscCall(PetscMalloc1(nrecvs, &iremote));
469: nrecvs = 0;
470: for (PetscMPIInt i = 0; i < nfrom; i++) {
471: for (j = 0; j < fromsizes[2 * fromperm_newtoold[i]]; j++) {
472: iremote[nrecvs].rank = fromranks[i];
473: iremote[nrecvs++].index = fromsizes[2 * fromperm_newtoold[i] + 1] + j;
474: }
475: }
476: PetscCall(PetscSFCreate(comm, &sf));
477: PetscCall(PetscSFSetGraph(sf, nsends, nrecvs, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
478: PetscCall(PetscSFSetType(sf, PETSCSFBASIC));
479: /* how to put a prefix ? */
480: PetscCall(PetscSFSetFromOptions(sf));
481: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
482: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
483: PetscCall(PetscSFDestroy(&sf));
484: PetscCall(PetscFree(fromranks));
485: PetscCall(PetscFree(fromsizes));
486: PetscCall(PetscFree(fromperm_newtoold));
487: PetscCall(PetscFree(send_indices));
488: if (rows) {
489: PetscCall(PetscSortInt(nrecvs, recv_indices));
490: PetscCall(ISCreateGeneral(comm, nrecvs, recv_indices, PETSC_OWN_POINTER, rows));
491: } else {
492: PetscCall(PetscFree(recv_indices));
493: }
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*@
498: ISPartitioningToNumbering - Takes an `IS' that represents a partitioning (the MPI rank that each local entry belongs to) and on each MPI process
499: generates an `IS` that contains a new global node number in the new ordering for each entry
501: Collective
503: Input Parameter:
504: . part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
506: Output Parameter:
507: . is - on each processor the index set that defines the global numbers
508: (in the new numbering) for all the nodes currently (before the partitioning)
509: on that processor
511: Level: advanced
513: Note:
514: The resulting `IS` tells where each local entry is mapped to in a new global ordering
516: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
517: @*/
518: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
519: {
520: MPI_Comm comm;
521: IS ndorder;
522: PetscInt n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
523: const PetscInt *indices = NULL;
524: PetscMPIInt np, npt;
526: PetscFunctionBegin;
528: PetscAssertPointer(is, 2);
529: /* see if the partitioning comes from nested dissection */
530: PetscCall(PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder));
531: if (ndorder) {
532: PetscCall(PetscObjectReference((PetscObject)ndorder));
533: *is = ndorder;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
538: /* count the number of partitions, i.e., virtual processors */
539: PetscCall(ISGetLocalSize(part, &n));
540: PetscCall(ISGetIndices(part, &indices));
541: np = 0;
542: for (PetscInt i = 0; i < n; i++) PetscCall(PetscMPIIntCast(PetscMax(np, indices[i]), &np));
543: PetscCallMPI(MPIU_Allreduce(&np, &npt, 1, MPI_INT, MPI_MAX, comm));
544: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
546: /*
547: lsizes - number of elements of each partition on this particular processor
548: sums - total number of "previous" nodes for any particular partition
549: starts - global number of first element in each partition on this processor
550: */
551: PetscCall(PetscMalloc3(np, &lsizes, np, &starts, np, &sums));
552: PetscCall(PetscArrayzero(lsizes, np));
553: for (PetscInt i = 0; i < n; i++) lsizes[indices[i]]++;
554: PetscCallMPI(MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm));
555: PetscCallMPI(MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm));
556: for (PetscMPIInt i = 0; i < np; i++) starts[i] -= lsizes[i];
557: for (PetscMPIInt i = 1; i < np; i++) {
558: sums[i] += sums[i - 1];
559: starts[i] += sums[i - 1];
560: }
562: /*
563: For each local index give it the new global number
564: */
565: PetscCall(PetscMalloc1(n, &newi));
566: for (PetscInt i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
567: PetscCall(PetscFree3(lsizes, starts, sums));
569: PetscCall(ISRestoreIndices(part, &indices));
570: PetscCall(ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is));
571: PetscCall(ISSetPermutation(*is));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*@
576: ISPartitioningCount - Takes a `IS` that represents a partitioning (the MPI rank that each local entry belongs to) and determines the number of
577: resulting elements on each (partition) rank
579: Collective
581: Input Parameters:
582: + part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
583: - len - length of the array count, this is the total number of partitions
585: Output Parameter:
586: . count - array of length size, to contain the number of elements assigned
587: to each partition, where size is the number of partitions generated
588: (see notes below).
590: Level: advanced
592: Notes:
593: By default the number of partitions generated (and thus the length
594: of count) is the size of the communicator associated with `IS`,
595: but it can be set by `MatPartitioningSetNParts()`.
597: The resulting array of lengths can for instance serve as input of `PCBJacobiSetTotalBlocks()`.
599: If the partitioning has been obtained by `MatPartitioningApplyND()`, the returned count does not include the separators.
601: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
602: `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
603: @*/
604: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
605: {
606: MPI_Comm comm;
607: PetscInt i, n, *lsizes;
608: const PetscInt *indices;
610: PetscFunctionBegin;
611: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
612: if (len == PETSC_DEFAULT) {
613: PetscMPIInt size;
615: PetscCallMPI(MPI_Comm_size(comm, &size));
616: len = size;
617: }
619: /* count the number of partitions */
620: PetscCall(ISGetLocalSize(part, &n));
621: PetscCall(ISGetIndices(part, &indices));
622: if (PetscDefined(USE_DEBUG)) {
623: PetscInt np = 0, npt;
624: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
625: PetscCallMPI(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
626: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
627: PetscCheck(np <= len, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Length of count array %" PetscInt_FMT " is less than number of partitions %" PetscInt_FMT, len, np);
628: }
630: /*
631: lsizes - number of elements of each partition on this particular processor
632: sums - total number of "previous" nodes for any particular partition
633: starts - global number of first element in each partition on this processor
634: */
635: PetscCall(PetscCalloc1(len, &lsizes));
636: for (i = 0; i < n; i++) {
637: if (indices[i] > -1) lsizes[indices[i]]++;
638: }
639: PetscCall(ISRestoreIndices(part, &indices));
640: PetscCallMPI(MPIU_Allreduce(lsizes, count, len, MPIU_INT, MPI_SUM, comm));
641: PetscCall(PetscFree(lsizes));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /*@
646: ISAllGather - Given an index set `IS` on each processor, generates a large
647: index set (same on each processor) by concatenating together each
648: processors index set.
650: Collective
652: Input Parameter:
653: . is - the distributed index set
655: Output Parameter:
656: . isout - the concatenated index set (same on all processors)
658: Level: intermediate
660: Notes:
661: `ISAllGather()` is clearly not scalable for large index sets.
663: The `IS` created on each processor must be created with a common
664: communicator (e.g., `PETSC_COMM_WORLD`). If the index sets were created
665: with `PETSC_COMM_SELF`, this routine will not work as expected, since
666: each process will generate its own new `IS` that consists only of
667: itself.
669: The communicator for this new `IS` is `PETSC_COMM_SELF`
671: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
672: @*/
673: PetscErrorCode ISAllGather(IS is, IS *isout)
674: {
675: PetscInt *indices, n, i, N, step, first;
676: const PetscInt *lindices;
677: MPI_Comm comm;
678: PetscMPIInt size, *sizes = NULL, *offsets = NULL, nn;
679: PetscBool stride;
681: PetscFunctionBegin;
683: PetscAssertPointer(isout, 2);
685: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
686: PetscCallMPI(MPI_Comm_size(comm, &size));
687: PetscCall(ISGetLocalSize(is, &n));
688: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride));
689: if (size == 1 && stride) { /* should handle parallel ISStride also */
690: PetscCall(ISStrideGetInfo(is, &first, &step));
691: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, first, step, isout));
692: } else {
693: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
695: PetscCall(PetscMPIIntCast(n, &nn));
696: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
697: offsets[0] = 0;
698: for (i = 1; i < size; i++) {
699: PetscInt s = offsets[i - 1] + sizes[i - 1];
700: PetscCall(PetscMPIIntCast(s, &offsets[i]));
701: }
702: N = offsets[size - 1] + sizes[size - 1];
704: PetscCall(PetscMalloc1(N, &indices));
705: PetscCall(ISGetIndices(is, &lindices));
706: PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm));
707: PetscCall(ISRestoreIndices(is, &lindices));
708: PetscCall(PetscFree2(sizes, offsets));
710: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout));
711: }
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: /*@C
716: ISAllGatherColors - Given a set of colors on each processor, generates a large
717: set (same on each processor) by concatenating together each processors colors
719: Collective
721: Input Parameters:
722: + comm - communicator to share the indices
723: . n - local size of set
724: - lindices - local colors
726: Output Parameters:
727: + outN - total number of indices
728: - outindices - all of the colors
730: Level: intermediate
732: Note:
733: `ISAllGatherColors()` is clearly not scalable for large index sets.
735: .seealso: `ISColoringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
736: @*/
737: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue lindices[], PetscInt *outN, ISColoringValue *outindices[])
738: {
739: ISColoringValue *indices;
740: PetscInt N;
741: PetscMPIInt size, *offsets = NULL, *sizes = NULL, nn;
743: PetscFunctionBegin;
744: PetscCall(PetscMPIIntCast(n, &nn));
745: PetscCallMPI(MPI_Comm_size(comm, &size));
746: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
748: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
749: offsets[0] = 0;
750: for (PetscMPIInt i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
751: N = offsets[size - 1] + sizes[size - 1];
752: PetscCall(PetscFree2(sizes, offsets));
754: PetscCall(PetscMalloc1(N + 1, &indices));
755: PetscCallMPI(MPI_Allgatherv(lindices, nn, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm));
757: *outindices = indices;
758: if (outN) *outN = N;
759: PetscFunctionReturn(PETSC_SUCCESS);
760: }
762: /*@
763: ISComplement - Given an index set `IS` generates the complement index set. That is
764: all indices that are NOT in the given set.
766: Collective
768: Input Parameters:
769: + is - the index set
770: . nmin - the first index desired in the local part of the complement
771: - nmax - the largest index desired in the local part of the complement (note that all indices in `is` must be greater or equal to `nmin` and less than `nmax`)
773: Output Parameter:
774: . isout - the complement
776: Level: intermediate
778: Notes:
779: The communicator for `isout` is the same as for the input `is`
781: For a parallel `is`, this will generate the local part of the complement on each process
783: To generate the entire complement (on each process) of a parallel `is`, first call `ISAllGather()` and then
784: call this routine.
786: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
787: @*/
788: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
789: {
790: const PetscInt *indices;
791: PetscInt n, i, j, unique, cnt, *nindices;
792: PetscBool sorted;
794: PetscFunctionBegin;
796: PetscAssertPointer(isout, 4);
797: PetscCheck(nmin >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be negative", nmin);
798: PetscCheck(nmin <= nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be greater than nmax %" PetscInt_FMT, nmin, nmax);
799: PetscCall(ISSorted(is, &sorted));
800: PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set must be sorted");
802: PetscCall(ISGetLocalSize(is, &n));
803: PetscCall(ISGetIndices(is, &indices));
804: if (PetscDefined(USE_DEBUG)) {
805: for (i = 0; i < n; i++) {
806: PetscCheck(indices[i] >= nmin, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT "'s value %" PetscInt_FMT " is smaller than minimum given %" PetscInt_FMT, i, indices[i], nmin);
807: PetscCheck(indices[i] < nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Index %" PetscInt_FMT "'s value %" PetscInt_FMT " is larger than maximum given %" PetscInt_FMT, i, indices[i], nmax);
808: }
809: }
810: /* Count number of unique entries */
811: unique = (n > 0);
812: for (i = 0; i < n - 1; i++) {
813: if (indices[i + 1] != indices[i]) unique++;
814: }
815: PetscCall(PetscMalloc1(nmax - nmin - unique, &nindices));
816: cnt = 0;
817: for (i = nmin, j = 0; i < nmax; i++) {
818: if (j < n && i == indices[j]) do {
819: j++;
820: } while (j < n && i == indices[j]);
821: else nindices[cnt++] = i;
822: }
823: PetscCheck(cnt == nmax - nmin - unique, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of entries found in complement %" PetscInt_FMT " does not match expected %" PetscInt_FMT, cnt, nmax - nmin - unique);
824: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout));
825: PetscCall(ISSetInfo(*isout, IS_SORTED, IS_GLOBAL, PETSC_FALSE, PETSC_TRUE));
826: PetscCall(ISRestoreIndices(is, &indices));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }