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: /*@
8: ISColoringReference - Increases the reference count of an `ISColoring` object by one
10: Logically collective
12: Input Parameter:
13: . coloring - the `ISColoring` object
15: Level: developer
17: Note:
18: The reference count is decreased by a matching call to `ISColoringDestroy()`.
20: .seealso: `ISColoring`, `ISColoringCreate()`, `ISColoringDestroy()`
21: @*/
22: PetscErrorCode ISColoringReference(ISColoring coloring)
23: {
24: PetscFunctionBegin;
25: coloring->refct++;
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: /*@
30: ISColoringSetType - indicates if the coloring is for the local representation (including ghost points) or the global representation of a `Mat`
32: Collective
34: Input Parameters:
35: + coloring - the coloring object
36: - type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`
38: Level: intermediate
40: Notes:
41: `IS_COLORING_LOCAL` can lead to faster computations since parallel ghost point updates are not needed for each color
43: 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
45: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringGetType()`
46: @*/
47: PetscErrorCode ISColoringSetType(ISColoring coloring, ISColoringType type)
48: {
49: PetscFunctionBegin;
50: coloring->ctype = type;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@
55: ISColoringGetType - gets if the coloring is for the local representation (including ghost points) or the global representation
57: Collective
59: Input Parameter:
60: . coloring - the coloring object
62: Output Parameter:
63: . type - either `IS_COLORING_LOCAL` or `IS_COLORING_GLOBAL`
65: Level: intermediate
67: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringSetType()`
68: @*/
69: PetscErrorCode ISColoringGetType(ISColoring coloring, ISColoringType *type)
70: {
71: PetscFunctionBegin;
72: *type = coloring->ctype;
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: ISColoringDestroy - Destroys an `ISColoring` coloring context.
79: Collective
81: Input Parameter:
82: . iscoloring - the coloring context
84: Level: advanced
86: .seealso: `ISColoring`, `ISColoringView()`, `MatColoring`
87: @*/
88: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
89: {
90: PetscInt i;
92: PetscFunctionBegin;
93: if (!*iscoloring) PetscFunctionReturn(PETSC_SUCCESS);
94: PetscAssertPointer(*iscoloring, 1);
95: if (--(*iscoloring)->refct > 0) {
96: *iscoloring = NULL;
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: if ((*iscoloring)->is) {
101: for (i = 0; i < (*iscoloring)->n; i++) PetscCall(ISDestroy(&(*iscoloring)->is[i]));
102: PetscCall(PetscFree((*iscoloring)->is));
103: }
104: if ((*iscoloring)->allocated) PetscCall(PetscFree((*iscoloring)->colors));
105: PetscCall(PetscCommDestroy(&(*iscoloring)->comm));
106: PetscCall(PetscFree(*iscoloring));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: /*@
111: ISColoringViewFromOptions - Processes command line options to determine if/how an `ISColoring` object is to be viewed.
113: Collective
115: Input Parameters:
116: + obj - the `ISColoring` object
117: . bobj - object that provides the options database key prefix to use for viewing, or `NULL` to use prefix of `obj`
118: - name - option to activate viewing
120: Options Database Key:
121: . -name viewer_specification - See `PetscOptionsCreateViewer()` for the values of `viewer_specification`
123: Level: intermediate
125: Note:
126: This checks the options database, creates the viewer on-the-fly, uses it and then destroys it. Hence it should not be called in heavily used routines,
127: rather `PetscOptionsCreateViewer()` should be used to construct the viewer once which can then be utilized in the heavily used routine.
129: Developer Note:
130: This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`
132: .seealso: `ISColoring`, `ISColoringView()`, `PetscObjectViewFromOptions()`, `PetscOptionsCreateViewer()`
133: @*/
134: PetscErrorCode ISColoringViewFromOptions(ISColoring obj, PetscObject bobj, const char name[])
135: {
136: PetscViewer viewer;
137: PetscBool flg;
138: PetscViewerFormat format;
139: char *prefix;
141: PetscFunctionBegin;
142: prefix = bobj ? bobj->prefix : NULL;
143: PetscCall(PetscOptionsCreateViewer(obj->comm, NULL, prefix, name, &viewer, &format, &flg));
144: if (flg) {
145: PetscCall(PetscViewerPushFormat(viewer, format));
146: PetscCall(ISColoringView(obj, viewer));
147: PetscCall(PetscViewerPopFormat(viewer));
148: PetscCall(PetscViewerDestroy(&viewer));
149: }
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /*@
154: ISColoringView - Views an `ISColoring` coloring context.
156: Collective
158: Input Parameters:
159: + iscoloring - the coloring context
160: - viewer - the viewer
162: Level: advanced
164: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
165: @*/
166: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
167: {
168: PetscInt i;
169: PetscBool isascii;
170: IS *is;
172: PetscFunctionBegin;
173: PetscAssertPointer(iscoloring, 1);
174: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(iscoloring->comm, &viewer));
177: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
178: if (isascii) {
179: MPI_Comm comm;
180: PetscMPIInt size, rank;
182: PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
183: PetscCallMPI(MPI_Comm_size(comm, &size));
184: PetscCallMPI(MPI_Comm_rank(comm, &rank));
185: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size));
186: PetscCall(PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]));
187: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
188: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n));
189: PetscCall(PetscViewerFlush(viewer));
190: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
191: }
193: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is));
194: for (i = 0; i < iscoloring->n; i++) PetscCall(ISView(iscoloring->is[i], viewer));
195: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: ISColoringGetColors - Returns an array with the color for each local node
202: Not Collective
204: Input Parameter:
205: . iscoloring - the coloring context
207: Output Parameters:
208: + n - number of nodes
209: . nc - number of colors
210: - colors - color for each node
212: Level: advanced
214: Notes:
215: Do not free the `colors` array.
217: The `colors` array will only be valid for the lifetime of the `ISColoring`
219: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
220: @*/
221: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
222: {
223: PetscFunctionBegin;
224: PetscAssertPointer(iscoloring, 1);
226: if (n) *n = iscoloring->N;
227: if (nc) *nc = iscoloring->n;
228: if (colors) *colors = iscoloring->colors;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*@
233: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
235: Collective
237: Input Parameters:
238: + iscoloring - the coloring context
239: - 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
241: Output Parameters:
242: + nn - number of index sets in the coloring context
243: - isis - array of index sets
245: Level: advanced
247: Note:
248: If mode is `PETSC_USE_POINTER` then `ISColoringRestoreIS()` must be called when the `IS` are no longer needed
250: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetType()`, `ISColoringGetColors()`
251: @*/
252: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
253: {
254: PetscFunctionBegin;
255: PetscAssertPointer(iscoloring, 1);
257: if (nn) *nn = iscoloring->n;
258: if (isis) {
259: if (!iscoloring->is) {
260: PetscInt *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
261: ISColoringValue *colors = iscoloring->colors;
262: IS *is;
264: if (PetscDefined(USE_DEBUG)) {
265: 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);
266: }
268: /* generate the lists of nodes for each color */
269: PetscCall(PetscCalloc1(nc, &mcolors));
270: for (i = 0; i < n; i++) mcolors[colors[i]]++;
272: PetscCall(PetscMalloc1(nc, &ii));
273: PetscCall(PetscMalloc1(n, &ii[0]));
274: for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
275: PetscCall(PetscArrayzero(mcolors, nc));
277: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
278: PetscCallMPI(MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm));
279: base -= iscoloring->N;
280: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
281: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
282: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
283: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");
285: PetscCall(PetscMalloc1(nc, &is));
286: for (i = 0; i < nc; i++) PetscCall(ISCreateGeneral(iscoloring->comm, mcolors[i], ii[i], PETSC_COPY_VALUES, is + i));
288: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
289: *isis = is;
290: PetscCall(PetscFree(ii[0]));
291: PetscCall(PetscFree(ii));
292: PetscCall(PetscFree(mcolors));
293: } else {
294: *isis = iscoloring->is;
295: if (mode == PETSC_OWN_POINTER) iscoloring->is = NULL;
296: }
297: }
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: /*@
302: ISColoringRestoreIS - Restores the index sets extracted from the coloring context with `ISColoringGetIS()` using `PETSC_USE_POINTER`
304: Collective
306: Input Parameters:
307: + iscoloring - the coloring context
308: . mode - who retains ownership of the is
309: - is - array of index sets
311: Level: advanced
313: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
314: @*/
315: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
316: {
317: PetscFunctionBegin;
318: PetscAssertPointer(iscoloring, 1);
320: /* currently nothing is done here */
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: ISColoringCreate - Generates an `ISColoring` context from lists (provided by each MPI process) of colors for each node.
327: Collective
329: Input Parameters:
330: + comm - communicator for the processors creating the coloring
331: . ncolors - maximum color value
332: . n - number of nodes on this processor
333: . colors - array containing the colors for this MPI process, color numbers begin at 0, for each local node
334: - mode - see `PetscCopyMode` for meaning of this flag.
336: Output Parameter:
337: . iscoloring - the resulting coloring data structure
339: Options Database Key:
340: . -is_coloring_view - Activates `ISColoringView()`
342: Level: advanced
344: Notes:
345: By default sets coloring type to `IS_COLORING_GLOBAL`
347: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
348: @*/
349: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
350: {
351: PetscMPIInt size, rank, tag;
352: PetscInt base, top, i;
353: PetscInt nc;
354: MPI_Status status;
356: PetscFunctionBegin;
357: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
358: 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);
359: 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);
360: }
361: PetscCall(PetscNew(iscoloring));
362: PetscCall(PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag));
363: comm = (*iscoloring)->comm;
365: /* compute the number of the first node on my processor */
366: PetscCallMPI(MPI_Comm_size(comm, &size));
368: /* should use MPI_Scan() */
369: PetscCallMPI(MPI_Comm_rank(comm, &rank));
370: if (rank == 0) {
371: base = 0;
372: top = n;
373: } else {
374: PetscCallMPI(MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status));
375: top = base + n;
376: }
377: if (rank < size - 1) PetscCallMPI(MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm));
379: /* compute the total number of colors */
380: nc = 0;
381: for (i = 0; i < n; i++) {
382: if (nc < colors[i]) nc = colors[i];
383: }
384: nc++;
385: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &nc, 1, MPIU_INT, MPI_MAX, comm));
386: 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);
387: (*iscoloring)->n = nc;
388: (*iscoloring)->is = NULL;
389: (*iscoloring)->N = n;
390: (*iscoloring)->refct = 1;
391: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
392: if (mode == PETSC_COPY_VALUES) {
393: PetscCall(PetscMalloc1(n, &(*iscoloring)->colors));
394: PetscCall(PetscArraycpy((*iscoloring)->colors, colors, n));
395: (*iscoloring)->allocated = PETSC_TRUE;
396: } else if (mode == PETSC_OWN_POINTER) {
397: (*iscoloring)->colors = (ISColoringValue *)colors;
398: (*iscoloring)->allocated = PETSC_TRUE;
399: } else {
400: (*iscoloring)->colors = (ISColoringValue *)colors;
401: (*iscoloring)->allocated = PETSC_FALSE;
402: }
403: PetscCall(ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view"));
404: PetscCall(PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: ISBuildTwoSided - Takes an `IS` that describes where each element will be mapped globally over all ranks.
410: Generates an `IS` that contains new numbers from remote or local on the `IS`.
412: Collective
414: Input Parameters:
415: + ito - an `IS` describes to which rank each entry will be mapped. Negative target rank will be ignored
416: - toindx - an `IS` describes what indices should send. `NULL` means sending natural numbering
418: Output Parameter:
419: . rows - contains new numbers from remote or local
421: Level: advanced
423: Developer Note:
424: This manual page is incomprehensible and still needs to be fixed
426: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
427: @*/
428: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
429: {
430: const PetscInt *ito_indices, *toindx_indices;
431: PetscInt *send_indices, rstart, *recv_indices, nrecvs, nsends;
432: PetscInt *tosizes, *fromsizes, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
433: PetscMPIInt *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
434: PetscLayout isrmap;
435: MPI_Comm comm;
436: PetscSF sf;
437: PetscSFNode *iremote;
439: PetscFunctionBegin;
440: PetscCall(PetscObjectGetComm((PetscObject)ito, &comm));
441: PetscCallMPI(MPI_Comm_size(comm, &size));
442: PetscCall(ISGetLocalSize(ito, &ito_ln));
443: PetscCall(ISGetLayout(ito, &isrmap));
444: PetscCall(PetscLayoutGetRange(isrmap, &rstart, NULL));
445: PetscCall(ISGetIndices(ito, &ito_indices));
446: PetscCall(PetscCalloc2(size, &tosizes_tmp, size + 1, &tooffsets_tmp));
447: for (PetscInt i = 0; i < ito_ln; i++) {
448: if (ito_indices[i] < 0) continue;
449: 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);
450: tosizes_tmp[ito_indices[i]]++;
451: }
452: nto = 0;
453: for (PetscMPIInt i = 0; i < size; i++) {
454: tooffsets_tmp[i + 1] = tooffsets_tmp[i] + tosizes_tmp[i];
455: if (tosizes_tmp[i] > 0) nto++;
456: }
457: PetscCall(PetscCalloc2(nto, &toranks, 2 * nto, &tosizes));
458: nto = 0;
459: for (PetscMPIInt i = 0; i < size; i++) {
460: if (tosizes_tmp[i] > 0) {
461: toranks[nto] = i;
462: tosizes[2 * nto] = tosizes_tmp[i]; /* size */
463: tosizes[2 * nto + 1] = tooffsets_tmp[i]; /* offset */
464: nto++;
465: }
466: }
467: nsends = tooffsets_tmp[size];
468: PetscCall(PetscCalloc1(nsends, &send_indices));
469: if (toindx) PetscCall(ISGetIndices(toindx, &toindx_indices));
470: for (PetscInt i = 0; i < ito_ln; i++) {
471: if (ito_indices[i] < 0) continue;
472: PetscCall(PetscMPIIntCast(ito_indices[i], &target_rank));
473: send_indices[tooffsets_tmp[target_rank]] = toindx ? toindx_indices[i] : (i + rstart);
474: tooffsets_tmp[target_rank]++;
475: }
476: if (toindx) PetscCall(ISRestoreIndices(toindx, &toindx_indices));
477: PetscCall(ISRestoreIndices(ito, &ito_indices));
478: PetscCall(PetscFree2(tosizes_tmp, tooffsets_tmp));
479: PetscCall(PetscCommBuildTwoSided(comm, 2, MPIU_INT, nto, toranks, tosizes, &nfrom, &fromranks, &fromsizes));
480: PetscCall(PetscFree2(toranks, tosizes));
481: PetscCall(PetscMalloc1(nfrom, &fromperm_newtoold));
482: for (PetscMPIInt i = 0; i < nfrom; i++) fromperm_newtoold[i] = i;
483: PetscCall(PetscSortMPIIntWithArray(nfrom, fromranks, fromperm_newtoold));
484: nrecvs = 0;
485: for (PetscMPIInt i = 0; i < nfrom; i++) nrecvs += fromsizes[i * 2];
486: PetscCall(PetscCalloc1(nrecvs, &recv_indices));
487: PetscCall(PetscMalloc1(nrecvs, &iremote));
488: nrecvs = 0;
489: for (PetscMPIInt i = 0; i < nfrom; i++) {
490: for (j = 0; j < fromsizes[2 * fromperm_newtoold[i]]; j++) {
491: iremote[nrecvs].rank = fromranks[i];
492: iremote[nrecvs++].index = fromsizes[2 * fromperm_newtoold[i] + 1] + j;
493: }
494: }
495: PetscCall(PetscSFCreate(comm, &sf));
496: PetscCall(PetscSFSetGraph(sf, nsends, nrecvs, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
497: PetscCall(PetscSFSetType(sf, PETSCSFBASIC));
498: /* how to put a prefix ? */
499: PetscCall(PetscSFSetFromOptions(sf));
500: PetscCall(PetscSFBcastBegin(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
501: PetscCall(PetscSFBcastEnd(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE));
502: PetscCall(PetscSFDestroy(&sf));
503: PetscCall(PetscFree(fromranks));
504: PetscCall(PetscFree(fromsizes));
505: PetscCall(PetscFree(fromperm_newtoold));
506: PetscCall(PetscFree(send_indices));
507: if (rows) {
508: PetscCall(PetscSortInt(nrecvs, recv_indices));
509: PetscCall(ISCreateGeneral(comm, nrecvs, recv_indices, PETSC_OWN_POINTER, rows));
510: } else {
511: PetscCall(PetscFree(recv_indices));
512: }
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: /*@
517: ISPartitioningToNumbering - Takes an `IS' that represents a partitioning (the MPI rank that each local entry belongs to) and on each MPI process
518: generates an `IS` that contains a new global node number in the new ordering for each entry
520: Collective
522: Input Parameter:
523: . part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
525: Output Parameter:
526: . is - on each processor the index set that defines the global numbers
527: (in the new numbering) for all the nodes currently (before the partitioning)
528: on that processor
530: Level: advanced
532: Note:
533: The resulting `IS` tells where each local entry is mapped to in a new global ordering
535: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
536: @*/
537: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
538: {
539: MPI_Comm comm;
540: IS ndorder;
541: PetscInt n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
542: const PetscInt *indices = NULL;
543: PetscMPIInt np;
545: PetscFunctionBegin;
547: PetscAssertPointer(is, 2);
548: /* see if the partitioning comes from nested dissection */
549: PetscCall(PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder));
550: if (ndorder) {
551: PetscCall(PetscObjectReference((PetscObject)ndorder));
552: *is = ndorder;
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
557: /* count the number of partitions, i.e., virtual processors */
558: PetscCall(ISGetLocalSize(part, &n));
559: PetscCall(ISGetIndices(part, &indices));
560: np = 0;
561: for (PetscInt i = 0; i < n; i++) PetscCall(PetscMPIIntCast(PetscMax(np, indices[i]), &np));
562: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &np, 1, MPI_INT, MPI_MAX, comm));
563: np++; /* so that it looks like a MPI_Comm_size output */
565: /*
566: lsizes - number of elements of each partition on this particular processor
567: sums - total number of "previous" nodes for any particular partition
568: starts - global number of first element in each partition on this processor
569: */
570: PetscCall(PetscMalloc3(np, &lsizes, np, &starts, np, &sums));
571: PetscCall(PetscArrayzero(lsizes, np));
572: for (PetscInt i = 0; i < n; i++) lsizes[indices[i]]++;
573: PetscCallMPI(MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm));
574: PetscCallMPI(MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm));
575: for (PetscMPIInt i = 0; i < np; i++) starts[i] -= lsizes[i];
576: for (PetscMPIInt i = 1; i < np; i++) {
577: sums[i] += sums[i - 1];
578: starts[i] += sums[i - 1];
579: }
581: /*
582: For each local index give it the new global number
583: */
584: PetscCall(PetscMalloc1(n, &newi));
585: for (PetscInt i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
586: PetscCall(PetscFree3(lsizes, starts, sums));
588: PetscCall(ISRestoreIndices(part, &indices));
589: PetscCall(ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is));
590: PetscCall(ISSetPermutation(*is));
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: /*@
595: ISPartitioningCount - Takes a `IS` that represents a partitioning (the MPI rank that each local entry belongs to) and determines the number of
596: resulting elements on each (partition) rank
598: Collective
600: Input Parameters:
601: + part - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
602: - len - length of the array count, this is the total number of partitions
604: Output Parameter:
605: . count - array of length size, to contain the number of elements assigned
606: to each partition, where size is the number of partitions generated
607: (see notes below).
609: Level: advanced
611: Notes:
612: By default the number of partitions generated (and thus the length
613: of count) is the size of the communicator associated with `IS`,
614: but it can be set by `MatPartitioningSetNParts()`.
616: The resulting array of lengths can for instance serve as input of `PCBJacobiSetTotalBlocks()`.
618: If the partitioning has been obtained by `MatPartitioningApplyND()`, the returned count does not include the separators.
620: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
621: `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
622: @*/
623: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
624: {
625: MPI_Comm comm;
626: PetscInt i, n;
627: const PetscInt *indices;
629: PetscFunctionBegin;
630: PetscCall(PetscObjectGetComm((PetscObject)part, &comm));
631: if (len == PETSC_DEFAULT) {
632: PetscMPIInt size;
634: PetscCallMPI(MPI_Comm_size(comm, &size));
635: len = size;
636: }
638: /* count the number of partitions */
639: PetscCall(ISGetLocalSize(part, &n));
640: PetscCall(ISGetIndices(part, &indices));
641: if (PetscDefined(USE_DEBUG)) {
642: PetscInt np = 0;
643: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
644: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &np, 1, MPIU_INT, MPI_MAX, comm));
645: np++; /* so that it looks like a MPI_Comm_size output */
646: 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);
647: }
649: /*
650: lsizes - number of elements of each partition on this particular processor
651: sums - total number of "previous" nodes for any particular partition
652: starts - global number of first element in each partition on this processor
653: */
654: PetscCall(PetscArrayzero(count, len));
655: for (i = 0; i < n; i++) {
656: if (indices[i] > -1) count[indices[i]]++;
657: }
658: PetscCall(ISRestoreIndices(part, &indices));
659: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, count, len, MPIU_INT, MPI_SUM, comm));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*@
664: ISAllGather - Given an index set `IS` on each processor, generates a large
665: index set (same on each processor) by concatenating together each
666: processors index set.
668: Collective
670: Input Parameter:
671: . is - the distributed index set
673: Output Parameter:
674: . isout - the concatenated index set (same on all processors)
676: Level: intermediate
678: Notes:
679: `ISAllGather()` is clearly not scalable for large index sets.
681: The `IS` created on each processor must be created with a common
682: communicator (e.g., `PETSC_COMM_WORLD`). If the index sets were created
683: with `PETSC_COMM_SELF`, this routine will not work as expected, since
684: each process will generate its own new `IS` that consists only of
685: itself.
687: The communicator for this new `IS` is `PETSC_COMM_SELF`
689: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
690: @*/
691: PetscErrorCode ISAllGather(IS is, IS *isout)
692: {
693: PetscInt *indices, n, i, N, step, first;
694: const PetscInt *lindices;
695: MPI_Comm comm;
696: PetscMPIInt size, *sizes = NULL, *offsets = NULL, nn;
697: PetscBool stride;
699: PetscFunctionBegin;
701: PetscAssertPointer(isout, 2);
703: PetscCall(PetscObjectGetComm((PetscObject)is, &comm));
704: PetscCallMPI(MPI_Comm_size(comm, &size));
705: PetscCall(ISGetLocalSize(is, &n));
706: PetscCall(PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride));
707: if (size == 1 && stride) { /* should handle parallel ISStride also */
708: PetscCall(ISStrideGetInfo(is, &first, &step));
709: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, first, step, isout));
710: } else {
711: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
713: PetscCall(PetscMPIIntCast(n, &nn));
714: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
715: offsets[0] = 0;
716: for (i = 1; i < size; i++) {
717: PetscInt s = offsets[i - 1] + sizes[i - 1];
718: PetscCall(PetscMPIIntCast(s, &offsets[i]));
719: }
720: N = offsets[size - 1] + sizes[size - 1];
722: PetscCall(PetscMalloc1(N, &indices));
723: PetscCall(ISGetIndices(is, &lindices));
724: PetscCallMPI(MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm));
725: PetscCall(ISRestoreIndices(is, &lindices));
726: PetscCall(PetscFree2(sizes, offsets));
728: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout));
729: }
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: /*@
734: ISAllGatherColors - Given a set of colors on each processor, generates a large
735: set (same on each processor) by concatenating together each processors colors
737: Collective
739: Input Parameters:
740: + comm - communicator to share the indices
741: . n - local size of set
742: - lindices - local colors
744: Output Parameters:
745: + outN - total number of indices
746: - outindices - all of the colors
748: Level: intermediate
750: Note:
751: `ISAllGatherColors()` is clearly not scalable for large index sets.
753: .seealso: `ISColoringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
754: @*/
755: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue lindices[], PetscInt *outN, ISColoringValue *outindices[])
756: {
757: ISColoringValue *indices;
758: PetscInt N;
759: PetscMPIInt size, *offsets = NULL, *sizes = NULL, nn;
761: PetscFunctionBegin;
762: PetscCall(PetscMPIIntCast(n, &nn));
763: PetscCallMPI(MPI_Comm_size(comm, &size));
764: PetscCall(PetscMalloc2(size, &sizes, size, &offsets));
766: PetscCallMPI(MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm));
767: offsets[0] = 0;
768: for (PetscMPIInt i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
769: N = offsets[size - 1] + sizes[size - 1];
770: PetscCall(PetscFree2(sizes, offsets));
772: PetscCall(PetscMalloc1(N + 1, &indices));
773: PetscCallMPI(MPI_Allgatherv(lindices, nn, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm));
775: *outindices = indices;
776: if (outN) *outN = N;
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: /*@
781: ISComplement - Given an index set `IS` generates the complement index set. That is
782: all indices that are NOT in the given set.
784: Collective
786: Input Parameters:
787: + is - the index set
788: . nmin - the first index desired in the local part of the complement
789: - 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`)
791: Output Parameter:
792: . isout - the complement
794: Level: intermediate
796: Notes:
797: The communicator for `isout` is the same as for the input `is`
799: For a parallel `is`, this will generate the local part of the complement on each process
801: To generate the entire complement (on each process) of a parallel `is`, first call `ISAllGather()` and then
802: call this routine.
804: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
805: @*/
806: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
807: {
808: const PetscInt *indices;
809: PetscInt n, i, j, unique, cnt, *nindices;
810: PetscBool sorted;
812: PetscFunctionBegin;
814: PetscAssertPointer(isout, 4);
815: PetscCheck(nmin >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be negative", nmin);
816: PetscCheck(nmin <= nmax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nmin %" PetscInt_FMT " cannot be greater than nmax %" PetscInt_FMT, nmin, nmax);
817: PetscCall(ISSorted(is, &sorted));
818: PetscCheck(sorted, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Index set must be sorted");
820: PetscCall(ISGetLocalSize(is, &n));
821: PetscCall(ISGetIndices(is, &indices));
822: if (PetscDefined(USE_DEBUG)) {
823: for (i = 0; i < n; i++) {
824: 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);
825: 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);
826: }
827: }
828: /* Count number of unique entries */
829: unique = (n > 0);
830: for (i = 0; i < n - 1; i++) {
831: if (indices[i + 1] != indices[i]) unique++;
832: }
833: PetscCall(PetscMalloc1(nmax - nmin - unique, &nindices));
834: cnt = 0;
835: for (i = nmin, j = 0; i < nmax; i++) {
836: if (j < n && i == indices[j]) do {
837: j++;
838: } while (j < n && i == indices[j]);
839: else nindices[cnt++] = i;
840: }
841: 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);
842: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout));
843: PetscCall(ISSetInfo(*isout, IS_SORTED, IS_GLOBAL, PETSC_FALSE, PETSC_TRUE));
844: PetscCall(ISRestoreIndices(is, &indices));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }