Actual source code: iscoloring.c
2: #include <petsc/private/isimpl.h>
3: #include <petscviewer.h>
4: #include <petscsf.h>
6: const char *const ISColoringTypes[] = {"global", "ghosted", "ISColoringType", "IS_COLORING_", NULL};
8: PetscErrorCode ISColoringReference(ISColoring coloring)
9: {
10: coloring->refct++;
11: return 0;
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: Notes:
24: `IS_COLORING_LOCAL` can lead to faster computations since parallel ghost point updates are not needed for each color
26: 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
28: Level: intermediate
30: .seealso: `MatFDColoringCreate()`, `ISColoring`, `ISColoringType`, `ISColoringCreate()`, `IS_COLORING_LOCAL`, `IS_COLORING_GLOBAL`, `ISColoringGetType()`
31: @*/
32: PetscErrorCode ISColoringSetType(ISColoring coloring, ISColoringType type)
33: {
34: coloring->ctype = type;
35: return 0;
36: }
38: /*@C
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: *type = coloring->ctype;
57: return 0;
58: }
60: /*@
61: ISColoringDestroy - Destroys an `ISColoring` coloring context.
63: Collective
65: Input Parameter:
66: . iscoloring - the coloring context
68: Level: advanced
70: .seealso: `ISColoring`, `ISColoringView()`, `MatColoring`
71: @*/
72: PetscErrorCode ISColoringDestroy(ISColoring *iscoloring)
73: {
74: PetscInt i;
76: if (!*iscoloring) return 0;
78: if (--(*iscoloring)->refct > 0) {
79: *iscoloring = NULL;
80: return 0;
81: }
83: if ((*iscoloring)->is) {
84: for (i = 0; i < (*iscoloring)->n; i++) ISDestroy(&(*iscoloring)->is[i]);
85: PetscFree((*iscoloring)->is);
86: }
87: if ((*iscoloring)->allocated) PetscFree((*iscoloring)->colors);
88: PetscCommDestroy(&(*iscoloring)->comm);
89: PetscFree((*iscoloring));
90: return 0;
91: }
93: /*@C
94: ISColoringViewFromOptions - Processes command line options to determine if/how an `ISColoring` object is to be viewed.
96: Collective
98: Input Parameters:
99: + obj - the `ISColoring` object
100: . prefix - prefix to use for viewing, or NULL to use prefix of 'mat'
101: - optionname - option to activate viewing
103: Level: intermediate
105: Developer Note:
106: This cannot use `PetscObjectViewFromOptions()` because `ISColoring` is not a `PetscObject`
108: .seealso: `ISColoring`, `ISColoringView()`
109: @*/
110: PetscErrorCode ISColoringViewFromOptions(ISColoring obj, PetscObject bobj, const char optionname[])
111: {
112: PetscViewer viewer;
113: PetscBool flg;
114: PetscViewerFormat format;
115: char *prefix;
117: prefix = bobj ? bobj->prefix : NULL;
118: PetscOptionsGetViewer(obj->comm, NULL, prefix, optionname, &viewer, &format, &flg);
119: if (flg) {
120: PetscViewerPushFormat(viewer, format);
121: ISColoringView(obj, viewer);
122: PetscViewerPopFormat(viewer);
123: PetscViewerDestroy(&viewer);
124: }
125: return 0;
126: }
128: /*@C
129: ISColoringView - Views an `ISColoring` coloring context.
131: Collective
133: Input Parameters:
134: + iscoloring - the coloring context
135: - viewer - the viewer
137: Level: advanced
139: .seealso: `ISColoring()`, `ISColoringViewFromOptions()`, `ISColoringDestroy()`, `ISColoringGetIS()`, `MatColoring`
140: @*/
141: PetscErrorCode ISColoringView(ISColoring iscoloring, PetscViewer viewer)
142: {
143: PetscInt i;
144: PetscBool iascii;
145: IS *is;
148: if (!viewer) PetscViewerASCIIGetStdout(iscoloring->comm, &viewer);
151: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
152: if (iascii) {
153: MPI_Comm comm;
154: PetscMPIInt size, rank;
156: PetscObjectGetComm((PetscObject)viewer, &comm);
157: MPI_Comm_size(comm, &size);
158: MPI_Comm_rank(comm, &rank);
159: PetscViewerASCIIPrintf(viewer, "ISColoring Object: %d MPI processes\n", size);
160: PetscViewerASCIIPrintf(viewer, "ISColoringType: %s\n", ISColoringTypes[iscoloring->ctype]);
161: PetscViewerASCIIPushSynchronized(viewer);
162: PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Number of colors %" PetscInt_FMT "\n", rank, iscoloring->n);
163: PetscViewerFlush(viewer);
164: PetscViewerASCIIPopSynchronized(viewer);
165: }
167: ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &is);
168: for (i = 0; i < iscoloring->n; i++) ISView(iscoloring->is[i], viewer);
169: ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &is);
170: return 0;
171: }
173: /*@C
174: ISColoringGetColors - Returns an array with the color for each local node
176: Not Collective
178: Input Parameter:
179: . iscoloring - the coloring context
181: Output Parameters:
182: + n - number of nodes
183: . nc - number of colors
184: - colors - color for each node
186: Level: advanced
188: Notes:
189: Do not free the colors[] array.
191: The colors[] array will only be valid for the lifetime of the `ISColoring`
193: .seealso: `ISColoring`, `ISColoringValue`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetIS()`
194: @*/
195: PetscErrorCode ISColoringGetColors(ISColoring iscoloring, PetscInt *n, PetscInt *nc, const ISColoringValue **colors)
196: {
199: if (n) *n = iscoloring->N;
200: if (nc) *nc = iscoloring->n;
201: if (colors) *colors = iscoloring->colors;
202: return 0;
203: }
205: /*@C
206: ISColoringGetIS - Extracts index sets from the coloring context. Each is contains the nodes of one color
208: Collective
210: Input Parameters:
211: + iscoloring - the coloring context
212: - 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
214: Output Parameters:
215: + nn - number of index sets in the coloring context
216: - is - array of index sets
218: Level: advanced
220: Note:
221: If mode is `PETSC_USE_POINTER` then `ISColoringRestoreIS()` must be called when the `IS` are no longer needed
223: .seealso: `ISColoring`, `IS`, `ISColoringRestoreIS()`, `ISColoringView()`, `ISColoringGetColoring()`, `ISColoringGetColors()`
224: @*/
225: PetscErrorCode ISColoringGetIS(ISColoring iscoloring, PetscCopyMode mode, PetscInt *nn, IS *isis[])
226: {
229: if (nn) *nn = iscoloring->n;
230: if (isis) {
231: if (!iscoloring->is) {
232: PetscInt *mcolors, **ii, nc = iscoloring->n, i, base, n = iscoloring->N;
233: ISColoringValue *colors = iscoloring->colors;
234: IS *is;
236: if (PetscDefined(USE_DEBUG)) {
238: }
240: /* generate the lists of nodes for each color */
241: PetscCalloc1(nc, &mcolors);
242: for (i = 0; i < n; i++) mcolors[colors[i]]++;
244: PetscMalloc1(nc, &ii);
245: PetscMalloc1(n, &ii[0]);
246: for (i = 1; i < nc; i++) ii[i] = ii[i - 1] + mcolors[i - 1];
247: PetscArrayzero(mcolors, nc);
249: if (iscoloring->ctype == IS_COLORING_GLOBAL) {
250: MPI_Scan(&iscoloring->N, &base, 1, MPIU_INT, MPI_SUM, iscoloring->comm);
251: base -= iscoloring->N;
252: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i + base; /* global idx */
253: } else if (iscoloring->ctype == IS_COLORING_LOCAL) {
254: for (i = 0; i < n; i++) ii[colors[i]][mcolors[colors[i]]++] = i; /* local idx */
255: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Not provided for this ISColoringType type");
257: PetscMalloc1(nc, &is);
258: for (i = 0; i < nc; i++) ISCreateGeneral(iscoloring->comm, mcolors[i], ii[i], PETSC_COPY_VALUES, is + i);
260: if (mode != PETSC_OWN_POINTER) iscoloring->is = is;
261: *isis = is;
262: PetscFree(ii[0]);
263: PetscFree(ii);
264: PetscFree(mcolors);
265: } else {
266: *isis = iscoloring->is;
267: }
268: }
269: return 0;
270: }
272: /*@C
273: ISColoringRestoreIS - Restores the index sets extracted from the coloring context with `ISColoringGetIS()` using `PETSC_USE_POINTER`
275: Collective
277: Input Parameters:
278: + iscoloring - the coloring context
279: . mode - who retains ownership of the is
280: - is - array of index sets
282: Level: advanced
284: .seealso: `ISColoring()`, `IS`, `ISColoringGetIS()`, `ISColoringView()`, `PetscCopyMode`
285: @*/
286: PetscErrorCode ISColoringRestoreIS(ISColoring iscoloring, PetscCopyMode mode, IS *is[])
287: {
290: /* currently nothing is done here */
291: return 0;
292: }
294: /*@
295: ISColoringCreate - Generates an `ISColoring` context from lists (provided by each MPI rank) of colors for each node.
297: Collective
299: Input Parameters:
300: + comm - communicator for the processors creating the coloring
301: . ncolors - max color value
302: . n - number of nodes on this processor
303: . colors - array containing the colors for this MPI rank, color numbers begin at 0, for each local node
304: - mode - see `PetscCopyMode` for meaning of this flag.
306: Output Parameter:
307: . iscoloring - the resulting coloring data structure
309: Options Database Key:
310: . -is_coloring_view - Activates `ISColoringView()`
312: Level: advanced
314: Notes:
315: By default sets coloring type to `IS_COLORING_GLOBAL`
317: .seealso: `ISColoring`, `ISColoringValue`, `MatColoringCreate()`, `ISColoringView()`, `ISColoringDestroy()`, `ISColoringSetType()`
318: @*/
319: PetscErrorCode ISColoringCreate(MPI_Comm comm, PetscInt ncolors, PetscInt n, const ISColoringValue colors[], PetscCopyMode mode, ISColoring *iscoloring)
320: {
321: PetscMPIInt size, rank, tag;
322: PetscInt base, top, i;
323: PetscInt nc, ncwork;
324: MPI_Status status;
326: if (ncolors != PETSC_DECIDE && ncolors > IS_COLORING_MAX) {
328: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Max color value exceeds limit. Perhaps reconfigure PETSc with --with-is-color-value-type=short?\n Current max: %d user requested: %" PetscInt_FMT, PETSC_IS_COLORING_MAX, ncolors);
329: }
330: PetscNew(iscoloring);
331: PetscCommDuplicate(comm, &(*iscoloring)->comm, &tag);
332: comm = (*iscoloring)->comm;
334: /* compute the number of the first node on my processor */
335: MPI_Comm_size(comm, &size);
337: /* should use MPI_Scan() */
338: MPI_Comm_rank(comm, &rank);
339: if (rank == 0) {
340: base = 0;
341: top = n;
342: } else {
343: MPI_Recv(&base, 1, MPIU_INT, rank - 1, tag, comm, &status);
344: top = base + n;
345: }
346: if (rank < size - 1) MPI_Send(&top, 1, MPIU_INT, rank + 1, tag, comm);
348: /* compute the total number of colors */
349: ncwork = 0;
350: for (i = 0; i < n; i++) {
351: if (ncwork < colors[i]) ncwork = colors[i];
352: }
353: ncwork++;
354: MPIU_Allreduce(&ncwork, &nc, 1, MPIU_INT, MPI_MAX, comm);
356: (*iscoloring)->n = nc;
357: (*iscoloring)->is = NULL;
358: (*iscoloring)->N = n;
359: (*iscoloring)->refct = 1;
360: (*iscoloring)->ctype = IS_COLORING_GLOBAL;
361: if (mode == PETSC_COPY_VALUES) {
362: PetscMalloc1(n, &(*iscoloring)->colors);
363: PetscArraycpy((*iscoloring)->colors, colors, n);
364: (*iscoloring)->allocated = PETSC_TRUE;
365: } else if (mode == PETSC_OWN_POINTER) {
366: (*iscoloring)->colors = (ISColoringValue *)colors;
367: (*iscoloring)->allocated = PETSC_TRUE;
368: } else {
369: (*iscoloring)->colors = (ISColoringValue *)colors;
370: (*iscoloring)->allocated = PETSC_FALSE;
371: }
372: ISColoringViewFromOptions(*iscoloring, NULL, "-is_coloring_view");
373: PetscInfo(0, "Number of colors %" PetscInt_FMT "\n", nc);
374: return 0;
375: }
377: /*@
378: ISBuildTwoSided - Takes an `IS` that describes where each element will be mapped globally over all ranks.
379: Generates an `IS` that contains new numbers from remote or local on the `IS`.
381: Collective
383: Input Parameters:
384: + ito - an `IS` describes where each entry will be mapped. Negative target rank will be ignored
385: - toindx - an `IS` describes what indices should send. NULL means sending natural numbering
387: Output Parameter:
388: . rows - contains new numbers from remote or local
390: Level: advanced
392: Developer Note:
393: This manual page is incomprehensible and needs to be fixed
395: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `ISPartitioningToNumbering()`, `ISPartitioningCount()`
396: @*/
397: PetscErrorCode ISBuildTwoSided(IS ito, IS toindx, IS *rows)
398: {
399: const PetscInt *ito_indices, *toindx_indices;
400: PetscInt *send_indices, rstart, *recv_indices, nrecvs, nsends;
401: PetscInt *tosizes, *fromsizes, i, j, *tosizes_tmp, *tooffsets_tmp, ito_ln;
402: PetscMPIInt *toranks, *fromranks, size, target_rank, *fromperm_newtoold, nto, nfrom;
403: PetscLayout isrmap;
404: MPI_Comm comm;
405: PetscSF sf;
406: PetscSFNode *iremote;
408: PetscObjectGetComm((PetscObject)ito, &comm);
409: MPI_Comm_size(comm, &size);
410: ISGetLocalSize(ito, &ito_ln);
411: ISGetLayout(ito, &isrmap);
412: PetscLayoutGetRange(isrmap, &rstart, NULL);
413: ISGetIndices(ito, &ito_indices);
414: PetscCalloc2(size, &tosizes_tmp, size + 1, &tooffsets_tmp);
415: for (i = 0; i < ito_ln; i++) {
416: if (ito_indices[i] < 0) continue;
418: tosizes_tmp[ito_indices[i]]++;
419: }
420: nto = 0;
421: for (i = 0; i < size; i++) {
422: tooffsets_tmp[i + 1] = tooffsets_tmp[i] + tosizes_tmp[i];
423: if (tosizes_tmp[i] > 0) nto++;
424: }
425: PetscCalloc2(nto, &toranks, 2 * nto, &tosizes);
426: nto = 0;
427: for (i = 0; i < size; i++) {
428: if (tosizes_tmp[i] > 0) {
429: toranks[nto] = i;
430: tosizes[2 * nto] = tosizes_tmp[i]; /* size */
431: tosizes[2 * nto + 1] = tooffsets_tmp[i]; /* offset */
432: nto++;
433: }
434: }
435: nsends = tooffsets_tmp[size];
436: PetscCalloc1(nsends, &send_indices);
437: if (toindx) ISGetIndices(toindx, &toindx_indices);
438: for (i = 0; i < ito_ln; i++) {
439: if (ito_indices[i] < 0) continue;
440: target_rank = ito_indices[i];
441: send_indices[tooffsets_tmp[target_rank]] = toindx ? toindx_indices[i] : (i + rstart);
442: tooffsets_tmp[target_rank]++;
443: }
444: if (toindx) ISRestoreIndices(toindx, &toindx_indices);
445: ISRestoreIndices(ito, &ito_indices);
446: PetscFree2(tosizes_tmp, tooffsets_tmp);
447: PetscCommBuildTwoSided(comm, 2, MPIU_INT, nto, toranks, tosizes, &nfrom, &fromranks, &fromsizes);
448: PetscFree2(toranks, tosizes);
449: PetscMalloc1(nfrom, &fromperm_newtoold);
450: for (i = 0; i < nfrom; i++) fromperm_newtoold[i] = i;
451: PetscSortMPIIntWithArray(nfrom, fromranks, fromperm_newtoold);
452: nrecvs = 0;
453: for (i = 0; i < nfrom; i++) nrecvs += fromsizes[i * 2];
454: PetscCalloc1(nrecvs, &recv_indices);
455: PetscMalloc1(nrecvs, &iremote);
456: nrecvs = 0;
457: for (i = 0; i < nfrom; i++) {
458: for (j = 0; j < fromsizes[2 * fromperm_newtoold[i]]; j++) {
459: iremote[nrecvs].rank = fromranks[i];
460: iremote[nrecvs++].index = fromsizes[2 * fromperm_newtoold[i] + 1] + j;
461: }
462: }
463: PetscSFCreate(comm, &sf);
464: PetscSFSetGraph(sf, nsends, nrecvs, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);
465: PetscSFSetType(sf, PETSCSFBASIC);
466: /* how to put a prefix ? */
467: PetscSFSetFromOptions(sf);
468: PetscSFBcastBegin(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE);
469: PetscSFBcastEnd(sf, MPIU_INT, send_indices, recv_indices, MPI_REPLACE);
470: PetscSFDestroy(&sf);
471: PetscFree(fromranks);
472: PetscFree(fromsizes);
473: PetscFree(fromperm_newtoold);
474: PetscFree(send_indices);
475: if (rows) {
476: PetscSortInt(nrecvs, recv_indices);
477: ISCreateGeneral(comm, nrecvs, recv_indices, PETSC_OWN_POINTER, rows);
478: } else {
479: PetscFree(recv_indices);
480: }
481: return 0;
482: }
484: /*@
485: ISPartitioningToNumbering - Takes an `IS' that represents a partitioning (the MPI rank that each local entry belongs to) and on each MPI rank
486: generates an `IS` that contains a new global node number in the new ordering for each entry
488: Collective
490: Input Parameters:
491: . partitioning - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
493: Output Parameter:
494: . is - on each processor the index set that defines the global numbers
495: (in the new numbering) for all the nodes currently (before the partitioning)
496: on that processor
498: Level: advanced
500: Note:
501: The resulting `IS` tells where each local entry is mapped to in a new global ordering
503: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningCount()`
504: @*/
505: PetscErrorCode ISPartitioningToNumbering(IS part, IS *is)
506: {
507: MPI_Comm comm;
508: IS ndorder;
509: PetscInt i, np, npt, n, *starts = NULL, *sums = NULL, *lsizes = NULL, *newi = NULL;
510: const PetscInt *indices = NULL;
514: /* see if the partitioning comes from nested dissection */
515: PetscObjectQuery((PetscObject)part, "_petsc_matpartitioning_ndorder", (PetscObject *)&ndorder);
516: if (ndorder) {
517: PetscObjectReference((PetscObject)ndorder);
518: *is = ndorder;
519: return 0;
520: }
522: PetscObjectGetComm((PetscObject)part, &comm);
523: /* count the number of partitions, i.e., virtual processors */
524: ISGetLocalSize(part, &n);
525: ISGetIndices(part, &indices);
526: np = 0;
527: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
528: MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm);
529: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
531: /*
532: lsizes - number of elements of each partition on this particular processor
533: sums - total number of "previous" nodes for any particular partition
534: starts - global number of first element in each partition on this processor
535: */
536: PetscMalloc3(np, &lsizes, np, &starts, np, &sums);
537: PetscArrayzero(lsizes, np);
538: for (i = 0; i < n; i++) lsizes[indices[i]]++;
539: MPIU_Allreduce(lsizes, sums, np, MPIU_INT, MPI_SUM, comm);
540: MPI_Scan(lsizes, starts, np, MPIU_INT, MPI_SUM, comm);
541: for (i = 0; i < np; i++) starts[i] -= lsizes[i];
542: for (i = 1; i < np; i++) {
543: sums[i] += sums[i - 1];
544: starts[i] += sums[i - 1];
545: }
547: /*
548: For each local index give it the new global number
549: */
550: PetscMalloc1(n, &newi);
551: for (i = 0; i < n; i++) newi[i] = starts[indices[i]]++;
552: PetscFree3(lsizes, starts, sums);
554: ISRestoreIndices(part, &indices);
555: ISCreateGeneral(comm, n, newi, PETSC_OWN_POINTER, is);
556: ISSetPermutation(*is);
557: return 0;
558: }
560: /*@
561: ISPartitioningCount - Takes a `IS` that represents a partitioning (the MPI rank that each local entry belongs to) and determines the number of
562: resulting elements on each (partition) rank
564: Collective
566: Input Parameters:
567: + partitioning - a partitioning as generated by `MatPartitioningApply()` or `MatPartitioningApplyND()`
568: - len - length of the array count, this is the total number of partitions
570: Output Parameter:
571: . count - array of length size, to contain the number of elements assigned
572: to each partition, where size is the number of partitions generated
573: (see notes below).
575: Level: advanced
577: Notes:
578: By default the number of partitions generated (and thus the length
579: of count) is the size of the communicator associated with `IS`,
580: but it can be set by `MatPartitioningSetNParts()`.
582: The resulting array of lengths can for instance serve as input of `PCBJacobiSetTotalBlocks()`.
584: If the partitioning has been obtained by `MatPartitioningApplyND()`, the returned count does not include the separators.
586: .seealso: [](sec_scatter), `IS`, `MatPartitioningCreate()`, `AOCreateBasic()`, `ISPartitioningToNumbering()`,
587: `MatPartitioningSetNParts()`, `MatPartitioningApply()`, `MatPartitioningApplyND()`
588: @*/
589: PetscErrorCode ISPartitioningCount(IS part, PetscInt len, PetscInt count[])
590: {
591: MPI_Comm comm;
592: PetscInt i, n, *lsizes;
593: const PetscInt *indices;
594: PetscMPIInt npp;
596: PetscObjectGetComm((PetscObject)part, &comm);
597: if (len == PETSC_DEFAULT) {
598: PetscMPIInt size;
599: MPI_Comm_size(comm, &size);
600: len = (PetscInt)size;
601: }
603: /* count the number of partitions */
604: ISGetLocalSize(part, &n);
605: ISGetIndices(part, &indices);
606: if (PetscDefined(USE_DEBUG)) {
607: PetscInt np = 0, npt;
608: for (i = 0; i < n; i++) np = PetscMax(np, indices[i]);
609: MPIU_Allreduce(&np, &npt, 1, MPIU_INT, MPI_MAX, comm);
610: np = npt + 1; /* so that it looks like a MPI_Comm_size output */
612: }
614: /*
615: lsizes - number of elements of each partition on this particular processor
616: sums - total number of "previous" nodes for any particular partition
617: starts - global number of first element in each partition on this processor
618: */
619: PetscCalloc1(len, &lsizes);
620: for (i = 0; i < n; i++) {
621: if (indices[i] > -1) lsizes[indices[i]]++;
622: }
623: ISRestoreIndices(part, &indices);
624: PetscMPIIntCast(len, &npp);
625: MPIU_Allreduce(lsizes, count, npp, MPIU_INT, MPI_SUM, comm);
626: PetscFree(lsizes);
627: return 0;
628: }
630: /*@
631: ISAllGather - Given an index set `IS` on each processor, generates a large
632: index set (same on each processor) by concatenating together each
633: processors index set.
635: Collective
637: Input Parameter:
638: . is - the distributed index set
640: Output Parameter:
641: . isout - the concatenated index set (same on all processors)
643: Level: intermediate
645: Notes:
646: `ISAllGather()` is clearly not scalable for large index sets.
648: The `IS` created on each processor must be created with a common
649: communicator (e.g., `PETSC_COMM_WORLD`). If the index sets were created
650: with `PETSC_COMM_SELF`, this routine will not work as expected, since
651: each process will generate its own new `IS` that consists only of
652: itself.
654: The communicator for this new `IS` is `PETSC_COMM_SELF`
656: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`
657: @*/
658: PetscErrorCode ISAllGather(IS is, IS *isout)
659: {
660: PetscInt *indices, n, i, N, step, first;
661: const PetscInt *lindices;
662: MPI_Comm comm;
663: PetscMPIInt size, *sizes = NULL, *offsets = NULL, nn;
664: PetscBool stride;
669: PetscObjectGetComm((PetscObject)is, &comm);
670: MPI_Comm_size(comm, &size);
671: ISGetLocalSize(is, &n);
672: PetscObjectTypeCompare((PetscObject)is, ISSTRIDE, &stride);
673: if (size == 1 && stride) { /* should handle parallel ISStride also */
674: ISStrideGetInfo(is, &first, &step);
675: ISCreateStride(PETSC_COMM_SELF, n, first, step, isout);
676: } else {
677: PetscMalloc2(size, &sizes, size, &offsets);
679: PetscMPIIntCast(n, &nn);
680: MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm);
681: offsets[0] = 0;
682: for (i = 1; i < size; i++) {
683: PetscInt s = offsets[i - 1] + sizes[i - 1];
684: PetscMPIIntCast(s, &offsets[i]);
685: }
686: N = offsets[size - 1] + sizes[size - 1];
688: PetscMalloc1(N, &indices);
689: ISGetIndices(is, &lindices);
690: MPI_Allgatherv((void *)lindices, nn, MPIU_INT, indices, sizes, offsets, MPIU_INT, comm);
691: ISRestoreIndices(is, &lindices);
692: PetscFree2(sizes, offsets);
694: ISCreateGeneral(PETSC_COMM_SELF, N, indices, PETSC_OWN_POINTER, isout);
695: }
696: return 0;
697: }
699: /*@C
700: ISAllGatherColors - Given a a set of colors on each processor, generates a large
701: set (same on each processor) by concatenating together each processors colors
703: Collective
705: Input Parameters:
706: + comm - communicator to share the indices
707: . n - local size of set
708: - lindices - local colors
710: Output Parameters:
711: + outN - total number of indices
712: - outindices - all of the colors
714: Level: intermediate
716: Note:
717: `ISAllGatherColors()` is clearly not scalable for large index sets.
719: .seealso: `ISCOloringValue`, `ISColoring()`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
720: @*/
721: PetscErrorCode ISAllGatherColors(MPI_Comm comm, PetscInt n, ISColoringValue *lindices, PetscInt *outN, ISColoringValue *outindices[])
722: {
723: ISColoringValue *indices;
724: PetscInt i, N;
725: PetscMPIInt size, *offsets = NULL, *sizes = NULL, nn = n;
727: MPI_Comm_size(comm, &size);
728: PetscMalloc2(size, &sizes, size, &offsets);
730: MPI_Allgather(&nn, 1, MPI_INT, sizes, 1, MPI_INT, comm);
731: offsets[0] = 0;
732: for (i = 1; i < size; i++) offsets[i] = offsets[i - 1] + sizes[i - 1];
733: N = offsets[size - 1] + sizes[size - 1];
734: PetscFree2(sizes, offsets);
736: PetscMalloc1(N + 1, &indices);
737: MPI_Allgatherv(lindices, (PetscMPIInt)n, MPIU_COLORING_VALUE, indices, sizes, offsets, MPIU_COLORING_VALUE, comm);
739: *outindices = indices;
740: if (outN) *outN = N;
741: return 0;
742: }
744: /*@
745: ISComplement - Given an index set `IS` generates the complement index set. That is all
746: all indices that are NOT in the given set.
748: Collective
750: Input Parameters:
751: + is - the index set
752: . nmin - the first index desired in the local part of the complement
753: - 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)
755: Output Parameter:
756: . isout - the complement
758: Level: intermediate
760: Notes:
761: The communicator for this new `IS` is the same as for the input IS
763: For a parallel `IS`, this will generate the local part of the complement on each process
765: To generate the entire complement (on each process) of a parallel `IS`, first call `ISAllGather()` and then
766: call this routine.
768: .seealso: [](sec_scatter), `IS`, `ISCreateGeneral()`, `ISCreateStride()`, `ISCreateBlock()`, `ISAllGather()`
769: @*/
770: PetscErrorCode ISComplement(IS is, PetscInt nmin, PetscInt nmax, IS *isout)
771: {
772: const PetscInt *indices;
773: PetscInt n, i, j, unique, cnt, *nindices;
774: PetscBool sorted;
780: ISSorted(is, &sorted);
783: ISGetLocalSize(is, &n);
784: ISGetIndices(is, &indices);
785: if (PetscDefined(USE_DEBUG)) {
786: for (i = 0; i < n; i++) {
789: }
790: }
791: /* Count number of unique entries */
792: unique = (n > 0);
793: for (i = 0; i < n - 1; i++) {
794: if (indices[i + 1] != indices[i]) unique++;
795: }
796: PetscMalloc1(nmax - nmin - unique, &nindices);
797: cnt = 0;
798: for (i = nmin, j = 0; i < nmax; i++) {
799: if (j < n && i == indices[j]) do {
800: j++;
801: } while (j < n && i == indices[j]);
802: else nindices[cnt++] = i;
803: }
805: ISCreateGeneral(PetscObjectComm((PetscObject)is), cnt, nindices, PETSC_OWN_POINTER, isout);
806: ISRestoreIndices(is, &indices);
807: return 0;
808: }