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: /*@C
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: /*@C
41: ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation
43: Collective
45: Input Parameter:
46: . coloring - the coloring object
48: Output Parameter:
49: . type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`
51: Level: intermediate
53: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringSetType()`
54: @*/
55: PetscErrorCode ISColoringGetType(ISColoring coloring, ISColoringType *type)
56: {
57: PetscFunctionBegin;
58: *type = coloring->ctype;
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: /*@
63: ISColoringDestroy - Destroys an `ISColoring` coloring context.
65: Collective
67: Input Parameter:
68: . iscoloring - the coloring context
70: Level: advanced
72: .seealso: `ISColoring`, `ISColoringView()`, `MatColoring`
73: @*/
74: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
75: {
76: PetscInt i;
78: PetscFunctionBegin;
79: if (!*iscoloring) PetscFunctionReturn(PETSC_SUCCESS);
80: PetscAssertPointer(*iscoloring, 1);
81: if (--(*iscoloring)->refct > 0) {
82: *iscoloring = NULL;
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: if ((*iscoloring)->is) {
87: for (i = 0; i < (*iscoloring)->n; i++) PetscCall(ISDestroy(&(*iscoloring)->is[i]));
88: PetscCall(PetscFree((*iscoloring)->is));
89: }
90: if ((*iscoloring)->allocated) PetscCall(PetscFree((*iscoloring)->colors));
91: PetscCall(PetscCommDestroy(&(*iscoloring)->comm));
92: PetscCall(PetscFree(*iscoloring));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@C
97: ISColoringViewFromOptions - Processes command line options to determine if/how an `ISColoring` object is to be viewed.
99: Collective
101: Input Parameters:
102: + obj - the `ISColoring` object
103: . bobj - prefix to use for viewing, or `NULL` to use prefix of `mat`
104: - optionname - option to activate viewing
106: Level: intermediate
108: Developer Notes:
109: This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`
111: .seealso: `ISColoring`, `ISColoringView()`
112: @*/
113: PetscErrorCode ISColoringViewFromOptions(ISColoring obj, PetscObject bobj, const char optionname[])
114: {
115: PetscViewer viewer;
116: PetscBool flg;
117: PetscViewerFormat format;
118: char *prefix;
120: PetscFunctionBegin;
121: prefix = bobj ? bobj->prefix : NULL;
122: PetscCall(PetscOptionsGetViewer(obj->comm, NULL, prefix, optionname, &viewer, &format, &flg));
123: if (flg) {
124: PetscCall(PetscViewerPushFormat(viewer, format));
125: PetscCall(ISColoringView(obj, viewer));
126: PetscCall(PetscViewerPopFormat(viewer));
127: PetscCall(PetscOptionsRestoreViewer(&viewer));
128: }
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@C
133: ISColoringView - Views an `ISColoring` coloring context.
135: Collective
137: Input Parameters:
138: + iscoloring - the coloring context
139: - viewer - the viewer
141: Level: advanced
143: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
144: @*/
145: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
146: {
147: PetscInt i;
148: PetscBool iascii;
149: IS *is;
151: PetscFunctionBegin;
152: PetscAssertPointer(iscoloring, 1);
153: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(iscoloring->comm, &viewer));
156: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
157: if (iascii) {
158: MPI_Comm comm;
159: PetscMPIInt size, rank;
161: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
162: PetscCallMPI(MPI_Comm_size(comm, &size));
163: PetscCallMPI(MPI_Comm_rank(comm, &rank));
164: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size));
165: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]));
166: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
167: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n));
168: PetscCall(PetscViewerFlush(viewer));
169: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
170: }
172: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is));
173: for (i = 0; i < iscoloring->n; i++) PetscCall(ISView(iscoloring->is[i], viewer));
174: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@C
179: ISColoringGetColors - Returns an array with the color for each local node
181: Not Collective
183: Input Parameter:
184: . iscoloring - the coloring context
186: Output Parameters:
187: + n - number of nodes
188: . nc - number of colors
189: - colors - color for each node
191: Level: advanced
193: Notes:
194: Do not free the `colors` array.
196: The `colors` array will only be valid for the lifetime of the `ISColoring`
198: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
199: @*/
200: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
201: {
202: PetscFunctionBegin;
203: PetscAssertPointer(iscoloring, 1);
205: if (n) *n = iscoloring->N;
206: if (nc) *nc = iscoloring->n;
207: if (colors) *colors = iscoloring->colors;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@C
212: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
214: Collective
216: Input Parameters:
217: + iscoloring - the coloring context
218: - 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
220: Output Parameters:
221: + nn - number of index sets in the coloring context
222: - isis - array of index sets
224: Level: advanced
226: Note:
227: If mode is `PETSC_USE_POINTER` then `ISColoringRestoreIS()` must be called when the `IS` are no longer needed
229: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetColoring()`, `ISColoringGetColors()`
230: @*/
231: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
232: {
233: PetscFunctionBegin;
234: PetscAssertPointer(iscoloring, 1);
236: if (nn) *nn = iscoloring->n;
237: if (isis) {
238: if (!iscoloring->is) {
239: PetscInt *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
240: ISColoringValue *colors = iscoloring->colors;
241: IS *is;
243: if (PetscDefined(USE_DEBUG)) {
244: for (i = 0; i < n; i++) PetscCheck(((PetscInt)colors[i]) < nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coloring is our of range index %d value %d number colors %d", (int)i, (int)colors[i], (int)nc);
245: }
247: /* generate the lists of nodes for each color */
248: PetscCall(PetscCalloc1(nc, &mcolors));
249: for (i = 0; i < n; i++) mcolors[colors[i]]++;
251: PetscCall(PetscMalloc1(nc, &ii));
252: PetscCall(PetscMalloc1(n, &ii[0]));
253: for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
254: PetscCall(PetscArrayzero(mcolors, nc));
256: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
257: PetscCallMPI(MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm));
258: base -= iscoloring->N;
259: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
260: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
261: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
262: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");
264: PetscCall(PetscMalloc1(nc, &is));
265: for (i = 0; i < nc; i++) PetscCall(ISCreateGeneral(iscoloring->comm, mcolors[i], ii[i], PETSC_COPY_VALUES, is + i));
267: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
268: *isis = is;
269: PetscCall(PetscFree(ii[0]));
270: PetscCall(PetscFree(ii));
271: PetscCall(PetscFree(mcolors));
272: } else {
273: *isis = iscoloring->is;
274: if (mode == PETSC_OWN_POINTER) iscoloring->is = NULL;
275: }
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@C
281: ISColoringRestoreIS - Restores the index sets extracted from the coloring context with `ISColoringGetIS()` using `PETSC_USE_POINTER`
283: Collective
285: Input Parameters:
286: + iscoloring - the coloring context
287: . mode - who retains ownership of the is
288: - is - array of index sets
290: Level: advanced
292: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
293: @*/
294: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
295: {
296: PetscFunctionBegin;
297: PetscAssertPointer(iscoloring, 1);
299: /* currently nothing is done here */
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*@
304: ISColoringCreate - Generates an `ISColoring` context from lists (provided by each MPI process) of colors for each node.
306: Collective
308: Input Parameters:
309: + comm - communicator for the processors creating the coloring
310: . ncolors - max color value
311: . n - number of nodes on this processor
312: . colors - array containing the colors for this MPI rank, color numbers begin at 0, for each local node
313: - mode - see `PetscCopyMode` for meaning of this flag.
315: Output Parameter:
316: . iscoloring - the resulting coloring data structure
318: Options Database Key:
319: . -is_coloring_view - Activates `ISColoringView()`
321: Level: advanced
323: Notes:
324: By default sets coloring type to `IS_COLORING_GLOBAL`
326: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
327: @*/
328: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
329: {
330: PetscMPIInt size, rank, tag;
331: PetscInt base, top, i;
332: PetscInt nc, ncwork;
333: MPI_Status status;
335: PetscFunctionBegin;
336: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
337: PetscCheck(ncolors <= PETSC_MAX_UINT16, 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_MAX_UINT16, PETSC_IS_COLORING_MAX, ncolors);
338: 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);
339: }
340: PetscCall(PetscNew(iscoloring));
341: PetscCall(PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag));
342: comm = (*iscoloring)->comm;
344: /* compute the number of the first node on my processor */
345: PetscCallMPI(MPI_Comm_size(comm, &size));
347: /* should use MPI_Scan() */
348: PetscCallMPI(MPI_Comm_rank(comm, &rank));
349: if (rank == 0) {
350: base = 0;
351: top = n;
352: } else {
353: PetscCallMPI(MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status));
354: top = base + n;
355: }
356: if (rank < size - 1) PetscCallMPI(MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm));
358: /* compute the total number of colors */
359: ncwork = 0;
360: for (i = 0; i < n; i++) {
361: if (ncwork < colors[i]) ncwork = colors[i];
362: }
363: ncwork++;
364: PetscCall(MPIU_Allreduce(&ncwork, &nc, 1, MPIU_INT, MPI_MAX, comm));
365: 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);
366: (*iscoloring)->n = nc;
367: (*iscoloring)->is = NULL;
368: (*iscoloring)->N = n;
369: (*iscoloring)->refct = 1;
370: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
371: if (mode == PETSC_COPY_VALUES) {
372: PetscCall(PetscMalloc1(n, &(*iscoloring)->colors));
373: PetscCall(PetscArraycpy((*iscoloring)->colors, colors, n));
374: (*iscoloring)->allocated = PETSC_TRUE;
375: } else if (mode == PETSC_OWN_POINTER) {
376: (*iscoloring)->colors = (ISColoringValue *)colors;
377: (*iscoloring)->allocated = PETSC_TRUE;
378: } else {
379: (*iscoloring)->colors = (ISColoringValue *)colors;
380: (*iscoloring)->allocated = PETSC_FALSE;
381: }
382: PetscCall(ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view"));
383: PetscCall(PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@
388: ISBuildTwoSided - Takes an `IS` that describes where each element will be mapped globally over all ranks.
389: Generates an `IS` that contains new numbers from remote or local on the `IS`.
391: Collective
393: Input Parameters:
394: + ito - an `IS` describes where each entry will be mapped. Negative target rank will be ignored
395: - toindx - an `IS` describes what indices should send. `NULL` means sending natural numbering
397: Output Parameter:
398: . rows - contains new numbers from remote or local
400: Level: advanced
402: Developer Notes:
403: This manual page is incomprehensible and still needs to be fixed
405: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
406: @*/
407: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
408: {
409: const PetscInt *ito_indices, *toindx_indices;
410: PetscInt *send_indices, rstart, *recv_indices, nrecvs, nsends;
411: PetscInt *tosizes, *fromsizes, i, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
412: PetscMPIInt *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
413: PetscLayout isrmap;
414: MPI_Comm comm;
415: PetscSF sf;
416: PetscSFNode *iremote;
418: PetscFunctionBegin;
419: PetscCall(PetscObjectGetComm((PetscObject)ito, &comm));
420: PetscCallMPI(MPI_Comm_size(comm, &size));
421: PetscCall(ISGetLocalSize(ito, &ito_ln));
422: PetscCall(ISGetLayout(ito, &isrmap));
423: PetscCall(PetscLayoutGetRange(isrmap, &rstart, NULL));
424: PetscCall(ISGetIndices(ito, &ito_indices));
425: PetscCall(PetscCalloc2(size, &tosizes_tmp, size + 1, &tooffsets_tmp));
426: for (i = 0; i < ito_ln; i++) {
427: if (ito_indices[i] < 0) continue;
428: 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);
429: tosizes_tmp[ito_indices[i]]++;
430: }
431: nto = 0;
432: for (i = 0; i < size; i++) {
433: tooffsets_tmp[i + 1] = tooffsets_tmp[i] + tosizes_tmp[i];
434: if (tosizes_tmp[i] > 0) nto++;
435: }
436: PetscCall(PetscCalloc2(nto, &toranks, 2 * nto, &tosizes));
437: nto = 0;
438: for (i = 0; i < size; i++) {
439: if (tosizes_tmp[i] > 0) {
440: toranks[nto] = i;
441: tosizes[2 * nto] = tosizes_tmp[i]; /* size */
442: tosizes[2 * nto + 1] = tooffsets_tmp[i]; /* offset */
443: nto++;
444: }
445: }
446: nsends = tooffsets_tmp[size];
447: PetscCall(PetscCalloc1(nsends, &send_indices));
448: if (toindx) PetscCall(ISGetIndices(toindx, &toindx_indices));
449: for (i = 0; i < ito_ln; i++) {
450: if (ito_indices[i] < 0) continue;
451: target_rank = ito_indices[i];
452: send_indices[tooffsets_tmp[target_rank]] = toindx ? toindx_indices[i] : (i + rstart);
453: tooffsets_tmp[target_rank]++;
454: }
455: if (toindx) PetscCall(ISRestoreIndices(toindx, &toindx_indices));
456: PetscCall(ISRestoreIndices(ito, &ito_indices));
457: PetscCall(PetscFree2(tosizes_tmp, tooffsets_tmp));
458: PetscCall(PetscCommBuildTwoSided(comm, 2, MPIU_INT, nto, toranks, tosizes, &nfrom, &fromranks, &fromsizes));
459: PetscCall(PetscFree2(toranks, tosizes));
460: PetscCall(PetscMalloc1(nfrom, &fromperm_newtoold));
461: for (i = 0; i < nfrom; i++) fromperm_newtoold[i] = i;
462: PetscCall(PetscSortMPIIntWithArray(nfrom, fromranks, fromperm_newtoold));
463: nrecvs = 0;
464: for (i = 0; i < nfrom; i++) nrecvs += fromsizes[i * 2];
465: PetscCall(PetscCalloc1(nrecvs, &recv_indices));
466: PetscCall(PetscMalloc1(nrecvs, &iremote));
467: nrecvs = 0;
468: for (i = 0; i < nfrom; i++) {
469: for (j = 0; j < fromsizes[2 * fromperm_newtoold[i]]; j++) {
470: iremote[nrecvs].rank = fromranks[i];
471: iremote[nrecvs++].index = fromsizes[2 * fromperm_newtoold[i] + 1] + j;
472: }
473: }
474: PetscCall(PetscSFCreate(comm, &sf));
475: PetscCall(PetscSFSetGraph(sf, nsends, nrecvs, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
476: PetscCall(PetscSFSetType(sf, PETSCSFBASIC));
477: /* how to put a prefix ? */
478: PetscCall(PetscSFSetFromOptions(sf));
479: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
480: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
481: PetscCall(PetscSFDestroy(&sf));
482: PetscCall(PetscFree(fromranks));
483: PetscCall(PetscFree(fromsizes));
484: PetscCall(PetscFree(fromperm_newtoold));
485: PetscCall(PetscFree(send_indices));
486: if (rows) {
487: PetscCall(PetscSortInt(nrecvs, recv_indices));
488: PetscCall(ISCreateGeneral(comm, nrecvs, recv_indices, PETSC_OWN_POINTER, rows));
489: } else {
490: PetscCall(PetscFree(recv_indices));
491: }
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@
496: ISPartitioningToNumbering - Takes an `IS' that represents a partitioning (the MPI rank that each local entry belongs to) and on each MPI process
497: generates an `IS` that contains a new global node number in the new ordering for each entry
499: Collective
501: Input Parameter:
502: . part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
504: Output Parameter:
505: . is - on each processor the index set that defines the global numbers
506: (in the new numbering) for all the nodes currently (before the partitioning)
507: on that processor
509: Level: advanced
511: Note:
512: The resulting `IS` tells where each local entry is mapped to in a new global ordering
514: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
515: @*/
516: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
517: {
518: MPI_Comm comm;
519: IS ndorder;
520: PetscInt i, np, npt, n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
521: const PetscInt *indices = NULL;
523: PetscFunctionBegin;
525: PetscAssertPointer(is, 2);
526: /* see if the partitioning comes from nested dissection */
527: PetscCall(PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder));
528: if (ndorder) {
529: PetscCall(PetscObjectReference((PetscObject)ndorder));
530: *is = ndorder;
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
535: /* count the number of partitions, i.e., virtual processors */
536: PetscCall(ISGetLocalSize(part, &n));
537: PetscCall(ISGetIndices(part, &indices));
538: np = 0;
539: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
540: PetscCall(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
541: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
543: /*
544: lsizes - number of elements of each partition on this particular processor
545: sums - total number of "previous" nodes for any particular partition
546: starts - global number of first element in each partition on this processor
547: */
548: PetscCall(PetscMalloc3(np, &lsizes, np, &starts, np, &sums));
549: PetscCall(PetscArrayzero(lsizes, np));
550: for (i = 0; i < n; i++) lsizes[indices[i]]++;
551: PetscCall(MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm));
552: PetscCallMPI(MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm));
553: for (i = 0; i < np; i++) starts[i] -= lsizes[i];
554: for (i = 1; i < np; i++) {
555: sums[i] += sums[i - 1];
556: starts[i] += sums[i - 1];
557: }
559: /*
560: For each local index give it the new global number
561: */
562: PetscCall(PetscMalloc1(n, &newi));
563: for (i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
564: PetscCall(PetscFree3(lsizes, starts, sums));
566: PetscCall(ISRestoreIndices(part, &indices));
567: PetscCall(ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is));
568: PetscCall(ISSetPermutation(*is));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@
573: ISPartitioningCount - Takes a `IS` that represents a partitioning (the MPI rank that each local entry belongs to) and determines the number of
574: resulting elements on each (partition) rank
576: Collective
578: Input Parameters:
579: + part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
580: - len - length of the array count, this is the total number of partitions
582: Output Parameter:
583: . count - array of length size, to contain the number of elements assigned
584: to each partition, where size is the number of partitions generated
585: (see notes below).
587: Level: advanced
589: Notes:
590: By default the number of partitions generated (and thus the length
591: of count) is the size of the communicator associated with `IS`,
592: but it can be set by `MatPartitioningSetNParts()`.
594: The resulting array of lengths can for instance serve as input of `PCBJacobiSetTotalBlocks()`.
596: If the partitioning has been obtained by `MatPartitioningApplyND()`, the returned count does not include the separators.
598: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
599: `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
600: @*/
601: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
602: {
603: MPI_Comm comm;
604: PetscInt i, n, *lsizes;
605: const PetscInt *indices;
606: PetscMPIInt npp;
608: PetscFunctionBegin;
609: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
610: if (len == PETSC_DEFAULT) {
611: PetscMPIInt size;
612: PetscCallMPI(MPI_Comm_size(comm, &size));
613: len = (PetscInt)size;
614: }
616: /* count the number of partitions */
617: PetscCall(ISGetLocalSize(part, &n));
618: PetscCall(ISGetIndices(part, &indices));
619: if (PetscDefined(USE_DEBUG)) {
620: PetscInt np = 0, npt;
621: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
622: PetscCall(MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm));
623: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
624: 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);
625: }
627: /*
628: lsizes - number of elements of each partition on this particular processor
629: sums - total number of "previous" nodes for any particular partition
630: starts - global number of first element in each partition on this processor
631: */
632: PetscCall(PetscCalloc1(len, &lsizes));
633: for (i = 0; i < n; i++) {
634: if (indices[i] > -1) lsizes[indices[i]]++;
635: }
636: PetscCall(ISRestoreIndices(part, &indices));
637: PetscCall(PetscMPIIntCast(len, &npp));
638: PetscCall(MPIU_Allreduce(lsizes, count, npp, MPIU_INT, MPI_SUM, comm));
639: PetscCall(PetscFree(lsizes));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: /*@
644: ISAllGather - Given an index set `IS` on each processor, generates a large
645: index set (same on each processor) by concatenating together each
646: processors index set.
648: Collective
650: Input Parameter:
651: . is - the distributed index set
653: Output Parameter:
654: . isout - the concatenated index set (same on all processors)
656: Level: intermediate
658: Notes:
659: `ISAllGather()` is clearly not scalable for large index sets.
661: The `IS` created on each processor must be created with a common
662: communicator (e.g., `PETSC_COMM_WORLD`). If the index sets were created
663: with `PETSC_COMM_SELF`, this routine will not work as expected, since
664: each process will generate its own new `IS` that consists only of
665: itself.
667: The communicator for this new `IS` is `PETSC_COMM_SELF`
669: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
670: @*/
671: PetscErrorCode ISAllGather(IS is, IS *isout)
672: {
673: PetscInt *indices, n, i, N, step, first;
674: const PetscInt *lindices;
675: MPI_Comm comm;
676: PetscMPIInt size, *sizes = NULL, *offsets = NULL, nn;
677: PetscBool stride;
679: PetscFunctionBegin;
681: PetscAssertPointer(isout, 2);
683: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
684: PetscCallMPI(MPI_Comm_size(comm, &size));
685: PetscCall(ISGetLocalSize(is, &n));
686: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride));
687: if (size == 1 && stride) { /* should handle parallel ISStride also */
688: PetscCall(ISStrideGetInfo(is, &first, &step));
689: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, first, step, isout));
690: } else {
691: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
693: PetscCall(PetscMPIIntCast(n, &nn));
694: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
695: offsets[0] = 0;
696: for (i = 1; i < size; i++) {
697: PetscInt s = offsets[i - 1] + sizes[i - 1];
698: PetscCall(PetscMPIIntCast(s, &offsets[i]));
699: }
700: N = offsets[size - 1] + sizes[size - 1];
702: PetscCall(PetscMalloc1(N, &indices));
703: PetscCall(ISGetIndices(is, &lindices));
704: PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm));
705: PetscCall(ISRestoreIndices(is, &lindices));
706: PetscCall(PetscFree2(sizes, offsets));
708: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout));
709: }
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: /*@C
714: ISAllGatherColors - Given a set of colors on each processor, generates a large
715: set (same on each processor) by concatenating together each processors colors
717: Collective
719: Input Parameters:
720: + comm - communicator to share the indices
721: . n - local size of set
722: - lindices - local colors
724: Output Parameters:
725: + outN - total number of indices
726: - outindices - all of the colors
728: Level: intermediate
730: Note:
731: `ISAllGatherColors()` is clearly not scalable for large index sets.
733: .seealso: `ISCOloringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
734: @*/
735: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue *lindices, PetscInt *outN, ISColoringValue *outindices[])
736: {
737: ISColoringValue *indices;
738: PetscInt i, N;
739: PetscMPIInt size, *offsets = NULL, *sizes = NULL, nn = n;
741: PetscFunctionBegin;
742: PetscCallMPI(MPI_Comm_size(comm, &size));
743: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
745: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
746: offsets[0] = 0;
747: for (i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
748: N = offsets[size - 1] + sizes[size - 1];
749: PetscCall(PetscFree2(sizes, offsets));
751: PetscCall(PetscMalloc1(N + 1, &indices));
752: PetscCallMPI(MPI_Allgatherv(lindices, (PetscMPIInt)n, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm));
754: *outindices = indices;
755: if (outN) *outN = N;
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: /*@
760: ISComplement - Given an index set `IS` generates the complement index set. That is
761: all indices that are NOT in the given set.
763: Collective
765: Input Parameters:
766: + is - the index set
767: . nmin - the first index desired in the local part of the complement
768: - 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`)
770: Output Parameter:
771: . isout - the complement
773: Level: intermediate
775: Notes:
776: The communicator for `isout` is the same as for the input `is`
778: For a parallel `is`, this will generate the local part of the complement on each process
780: To generate the entire complement (on each process) of a parallel `is`, first call `ISAllGather()` and then
781: call this routine.
783: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
784: @*/
785: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
786: {
787: const PetscInt *indices;
788: PetscInt n, i, j, unique, cnt, *nindices;
789: PetscBool sorted;
791: PetscFunctionBegin;
793: PetscAssertPointer(isout, 4);
794: PetscCheck(nmin >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be negative", nmin);
795: PetscCheck(nmin <= nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be greater than nmax %" PetscInt_FMT, nmin, nmax);
796: PetscCall(ISSorted(is, &sorted));
797: PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set must be sorted");
799: PetscCall(ISGetLocalSize(is, &n));
800: PetscCall(ISGetIndices(is, &indices));
801: if (PetscDefined(USE_DEBUG)) {
802: for (i = 0; i < n; i++) {
803: 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);
804: 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);
805: }
806: }
807: /* Count number of unique entries */
808: unique = (n > 0);
809: for (i = 0; i < n - 1; i++) {
810: if (indices[i + 1] != indices[i]) unique++;
811: }
812: PetscCall(PetscMalloc1(nmax - nmin - unique, &nindices));
813: cnt = 0;
814: for (i = nmin, j = 0; i < nmax; i++) {
815: if (j < n && i == indices[j]) do {
816: j++;
817: } while (j < n && i == indices[j]);
818: else nindices[cnt++] = i;
819: }
820: 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);
821: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout));
822: PetscCall(ISRestoreIndices(is, &indices));
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }