Actual source code: isdiff.c
1: #include <petsc/private/isimpl.h>
2: #include <petsc/private/sectionimpl.h>
3: #include <petscbt.h>
5: /*@
6: ISDifference - Computes the difference between two index sets.
8: Collective
10: Input Parameters:
11: + is1 - first index, to have items removed from it
12: - is2 - index values to be removed
14: Output Parameter:
15: . isout - is1 - is2
17: Level: intermediate
19: Notes:
20: Negative values are removed from the lists. `is2` may have values
21: that are not in `is1`.
23: This computation requires O(imax-imin) memory and O(imax-imin)
24: work, where imin and imax are the bounds on the indices in is1.
26: If `is2` is `NULL`, the result is the same as for an empty `IS`, i.e., a duplicate of `is1`.
28: The difference is computed separately on each MPI rank
30: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISSum()`, `ISExpand()`
31: @*/
32: PetscErrorCode ISDifference(IS is1, IS is2, IS *isout)
33: {
34: PetscInt i, n1, n2, imin, imax, nout, *iout;
35: const PetscInt *i1, *i2;
36: PetscBT mask;
37: MPI_Comm comm;
39: PetscFunctionBegin;
41: PetscAssertPointer(isout, 3);
42: if (!is2) {
43: PetscCall(ISDuplicate(is1, isout));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
48: PetscCall(ISGetIndices(is1, &i1));
49: PetscCall(ISGetLocalSize(is1, &n1));
51: /* Create a bit mask array to contain required values */
52: if (n1) {
53: imin = PETSC_MAX_INT;
54: imax = 0;
55: for (i = 0; i < n1; i++) {
56: if (i1[i] < 0) continue;
57: imin = PetscMin(imin, i1[i]);
58: imax = PetscMax(imax, i1[i]);
59: }
60: } else imin = imax = 0;
62: PetscCall(PetscBTCreate(imax - imin, &mask));
63: /* Put the values from is1 */
64: for (i = 0; i < n1; i++) {
65: if (i1[i] < 0) continue;
66: PetscCall(PetscBTSet(mask, i1[i] - imin));
67: }
68: PetscCall(ISRestoreIndices(is1, &i1));
69: /* Remove the values from is2 */
70: PetscCall(ISGetIndices(is2, &i2));
71: PetscCall(ISGetLocalSize(is2, &n2));
72: for (i = 0; i < n2; i++) {
73: if (i2[i] < imin || i2[i] > imax) continue;
74: PetscCall(PetscBTClear(mask, i2[i] - imin));
75: }
76: PetscCall(ISRestoreIndices(is2, &i2));
78: /* Count the number in the difference */
79: nout = 0;
80: for (i = 0; i < imax - imin + 1; i++) {
81: if (PetscBTLookup(mask, i)) nout++;
82: }
84: /* create the new IS containing the difference */
85: PetscCall(PetscMalloc1(nout, &iout));
86: nout = 0;
87: for (i = 0; i < imax - imin + 1; i++) {
88: if (PetscBTLookup(mask, i)) iout[nout++] = i + imin;
89: }
90: PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
91: PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
93: PetscCall(PetscBTDestroy(&mask));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: /*@
98: ISSum - Computes the sum (union) of two index sets.
100: Only sequential version (at the moment)
102: Input Parameters:
103: + is1 - index set to be extended
104: - is2 - index values to be added
106: Output Parameter:
107: . is3 - the sum; this can not be `is1` or `is2`
109: Level: intermediate
111: Notes:
112: If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;
114: Both index sets need to be sorted on input.
116: The sum is computed separately on each MPI rank
118: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISExpand()`
119: @*/
120: PetscErrorCode ISSum(IS is1, IS is2, IS *is3)
121: {
122: MPI_Comm comm;
123: PetscBool f;
124: PetscMPIInt size;
125: const PetscInt *i1, *i2;
126: PetscInt n1, n2, n3, p1, p2, *iout;
128: PetscFunctionBegin;
131: PetscCall(PetscObjectGetComm((PetscObject)(is1), &comm));
132: PetscCallMPI(MPI_Comm_size(comm, &size));
133: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Currently only for uni-processor IS");
135: PetscCall(ISSorted(is1, &f));
136: PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 1 is not sorted");
137: PetscCall(ISSorted(is2, &f));
138: PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 2 is not sorted");
140: PetscCall(ISGetLocalSize(is1, &n1));
141: PetscCall(ISGetLocalSize(is2, &n2));
142: if (!n2) {
143: PetscCall(ISDuplicate(is1, is3));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
146: PetscCall(ISGetIndices(is1, &i1));
147: PetscCall(ISGetIndices(is2, &i2));
149: p1 = 0;
150: p2 = 0;
151: n3 = 0;
152: do {
153: if (p1 == n1) { /* cleanup for is2 */
154: n3 += n2 - p2;
155: break;
156: } else {
157: while (p2 < n2 && i2[p2] < i1[p1]) {
158: n3++;
159: p2++;
160: }
161: if (p2 == n2) {
162: /* cleanup for is1 */
163: n3 += n1 - p1;
164: break;
165: } else {
166: if (i2[p2] == i1[p1]) {
167: n3++;
168: p1++;
169: p2++;
170: }
171: }
172: }
173: if (p2 == n2) {
174: /* cleanup for is1 */
175: n3 += n1 - p1;
176: break;
177: } else {
178: while (p1 < n1 && i1[p1] < i2[p2]) {
179: n3++;
180: p1++;
181: }
182: if (p1 == n1) {
183: /* clean up for is2 */
184: n3 += n2 - p2;
185: break;
186: } else {
187: if (i1[p1] == i2[p2]) {
188: n3++;
189: p1++;
190: p2++;
191: }
192: }
193: }
194: } while (p1 < n1 || p2 < n2);
196: if (n3 == n1) { /* no new elements to be added */
197: PetscCall(ISRestoreIndices(is1, &i1));
198: PetscCall(ISRestoreIndices(is2, &i2));
199: PetscCall(ISDuplicate(is1, is3));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
202: PetscCall(PetscMalloc1(n3, &iout));
204: p1 = 0;
205: p2 = 0;
206: n3 = 0;
207: do {
208: if (p1 == n1) { /* cleanup for is2 */
209: while (p2 < n2) iout[n3++] = i2[p2++];
210: break;
211: } else {
212: while (p2 < n2 && i2[p2] < i1[p1]) iout[n3++] = i2[p2++];
213: if (p2 == n2) { /* cleanup for is1 */
214: while (p1 < n1) iout[n3++] = i1[p1++];
215: break;
216: } else {
217: if (i2[p2] == i1[p1]) {
218: iout[n3++] = i1[p1++];
219: p2++;
220: }
221: }
222: }
223: if (p2 == n2) { /* cleanup for is1 */
224: while (p1 < n1) iout[n3++] = i1[p1++];
225: break;
226: } else {
227: while (p1 < n1 && i1[p1] < i2[p2]) iout[n3++] = i1[p1++];
228: if (p1 == n1) { /* clean up for is2 */
229: while (p2 < n2) iout[n3++] = i2[p2++];
230: break;
231: } else {
232: if (i1[p1] == i2[p2]) {
233: iout[n3++] = i1[p1++];
234: p2++;
235: }
236: }
237: }
238: } while (p1 < n1 || p2 < n2);
240: PetscCall(ISRestoreIndices(is1, &i1));
241: PetscCall(ISRestoreIndices(is2, &i2));
242: PetscCall(ISCreateGeneral(comm, n3, iout, PETSC_OWN_POINTER, is3));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*@
247: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
248: removing duplicates.
250: Collective
252: Input Parameters:
253: + is1 - first index set
254: - is2 - index values to be added
256: Output Parameter:
257: . isout - `is1` + `is2` The index set `is2` is appended to `is1` removing duplicates
259: Level: intermediate
261: Notes:
262: Negative values are removed from the lists. This requires O(imax-imin)
263: memory and O(imax-imin) work, where imin and imax are the bounds on the
264: indices in `is1` and `is2`.
266: `is1` and `is2` do not need to be sorted.
268: The operations are performed separately on each MPI rank
270: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISIntersect()`
271: @*/
272: PetscErrorCode ISExpand(IS is1, IS is2, IS *isout)
273: {
274: PetscInt i, n1, n2, imin, imax, nout, *iout;
275: const PetscInt *i1, *i2;
276: PetscBT mask;
277: MPI_Comm comm;
279: PetscFunctionBegin;
282: PetscAssertPointer(isout, 3);
284: PetscCheck(is1 || is2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
285: if (!is1) {
286: PetscCall(ISDuplicate(is2, isout));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
289: if (!is2) {
290: PetscCall(ISDuplicate(is1, isout));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
293: PetscCall(ISGetIndices(is1, &i1));
294: PetscCall(ISGetLocalSize(is1, &n1));
295: PetscCall(ISGetIndices(is2, &i2));
296: PetscCall(ISGetLocalSize(is2, &n2));
298: /* Create a bit mask array to contain required values */
299: if (n1 || n2) {
300: imin = PETSC_MAX_INT;
301: imax = 0;
302: for (i = 0; i < n1; i++) {
303: if (i1[i] < 0) continue;
304: imin = PetscMin(imin, i1[i]);
305: imax = PetscMax(imax, i1[i]);
306: }
307: for (i = 0; i < n2; i++) {
308: if (i2[i] < 0) continue;
309: imin = PetscMin(imin, i2[i]);
310: imax = PetscMax(imax, i2[i]);
311: }
312: } else imin = imax = 0;
314: PetscCall(PetscMalloc1(n1 + n2, &iout));
315: nout = 0;
316: PetscCall(PetscBTCreate(imax - imin, &mask));
317: /* Put the values from is1 */
318: for (i = 0; i < n1; i++) {
319: if (i1[i] < 0) continue;
320: if (!PetscBTLookupSet(mask, i1[i] - imin)) iout[nout++] = i1[i];
321: }
322: PetscCall(ISRestoreIndices(is1, &i1));
323: /* Put the values from is2 */
324: for (i = 0; i < n2; i++) {
325: if (i2[i] < 0) continue;
326: if (!PetscBTLookupSet(mask, i2[i] - imin)) iout[nout++] = i2[i];
327: }
328: PetscCall(ISRestoreIndices(is2, &i2));
330: /* create the new IS containing the sum */
331: PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
332: PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
334: PetscCall(PetscBTDestroy(&mask));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: /*@
339: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
341: Collective
343: Input Parameters:
344: + is1 - first index set
345: - is2 - second index set
347: Output Parameter:
348: . isout - the sorted intersection of `is1` and `is2`
350: Level: intermediate
352: Notes:
353: Negative values are removed from the lists. This requires O(min(is1,is2))
354: memory and O(max(is1,is2)log(max(is1,is2))) work
356: `is1` and `is2` do not need to be sorted.
358: The operations are performed separately on each MPI rank
360: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISConcatenate()`
361: @*/
362: PetscErrorCode ISIntersect(IS is1, IS is2, IS *isout)
363: {
364: PetscInt i, n1, n2, nout, *iout;
365: const PetscInt *i1, *i2;
366: IS is1sorted = NULL, is2sorted = NULL;
367: PetscBool sorted, lsorted;
368: MPI_Comm comm;
370: PetscFunctionBegin;
373: PetscCheckSameComm(is1, 1, is2, 2);
374: PetscAssertPointer(isout, 3);
375: PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
377: PetscCall(ISGetLocalSize(is1, &n1));
378: PetscCall(ISGetLocalSize(is2, &n2));
379: if (n1 < n2) {
380: IS tempis = is1;
381: PetscInt ntemp = n1;
383: is1 = is2;
384: is2 = tempis;
385: n1 = n2;
386: n2 = ntemp;
387: }
388: PetscCall(ISSorted(is1, &lsorted));
389: PetscCall(MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm));
390: if (!sorted) {
391: PetscCall(ISDuplicate(is1, &is1sorted));
392: PetscCall(ISSort(is1sorted));
393: PetscCall(ISGetIndices(is1sorted, &i1));
394: } else {
395: is1sorted = is1;
396: PetscCall(PetscObjectReference((PetscObject)is1));
397: PetscCall(ISGetIndices(is1, &i1));
398: }
399: PetscCall(ISSorted(is2, &lsorted));
400: PetscCall(MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm));
401: if (!sorted) {
402: PetscCall(ISDuplicate(is2, &is2sorted));
403: PetscCall(ISSort(is2sorted));
404: PetscCall(ISGetIndices(is2sorted, &i2));
405: } else {
406: is2sorted = is2;
407: PetscCall(PetscObjectReference((PetscObject)is2));
408: PetscCall(ISGetIndices(is2, &i2));
409: }
411: PetscCall(PetscMalloc1(n2, &iout));
413: for (i = 0, nout = 0; i < n2; i++) {
414: PetscInt key = i2[i];
415: PetscInt loc;
417: PetscCall(ISLocate(is1sorted, key, &loc));
418: if (loc >= 0) {
419: if (!nout || iout[nout - 1] < key) iout[nout++] = key;
420: }
421: }
422: PetscCall(PetscRealloc(nout * sizeof(PetscInt), &iout));
424: /* create the new IS containing the sum */
425: PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
427: PetscCall(ISRestoreIndices(is2sorted, &i2));
428: PetscCall(ISDestroy(&is2sorted));
429: PetscCall(ISRestoreIndices(is1sorted, &i1));
430: PetscCall(ISDestroy(&is1sorted));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
435: {
436: PetscFunctionBegin;
437: *isect = NULL;
438: if (is2 && is1) {
439: char composeStr[33] = {0};
440: PetscObjectId is2id;
442: PetscCall(PetscObjectGetId((PetscObject)is2, &is2id));
443: PetscCall(PetscSNPrintf(composeStr, 32, "ISIntersect_Caching_%" PetscInt64_FMT, is2id));
444: PetscCall(PetscObjectQuery((PetscObject)is1, composeStr, (PetscObject *)isect));
445: if (*isect == NULL) {
446: PetscCall(ISIntersect(is1, is2, isect));
447: PetscCall(PetscObjectCompose((PetscObject)is1, composeStr, (PetscObject)*isect));
448: } else {
449: PetscCall(PetscObjectReference((PetscObject)*isect));
450: }
451: }
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: /*@
456: ISConcatenate - Forms a new `IS` by locally concatenating the indices from an `IS` list without reordering.
458: Collective
460: Input Parameters:
461: + comm - communicator of the concatenated `IS`.
462: . len - size of islist array (nonnegative)
463: - islist - array of index sets
465: Output Parameter:
466: . isout - The concatenated index set; empty, if `len` == 0.
468: Level: intermediate
470: Notes:
471: The semantics of calling this on comm imply that the comms of the members of `islist` also contain this rank.
473: .seealso: [](sec_scatter), `IS`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISIntersect()`
474: @*/
475: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
476: {
477: PetscInt i, n, N;
478: const PetscInt *iidx;
479: PetscInt *idx;
481: PetscFunctionBegin;
482: if (len) PetscAssertPointer(islist, 3);
483: if (PetscDefined(USE_DEBUG)) {
484: for (i = 0; i < len; ++i)
486: }
487: PetscAssertPointer(isout, 4);
488: if (!len) {
489: PetscCall(ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
492: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %" PetscInt_FMT, len);
493: N = 0;
494: for (i = 0; i < len; ++i) {
495: if (islist[i]) {
496: PetscCall(ISGetLocalSize(islist[i], &n));
497: N += n;
498: }
499: }
500: PetscCall(PetscMalloc1(N, &idx));
501: N = 0;
502: for (i = 0; i < len; ++i) {
503: if (islist[i]) {
504: PetscCall(ISGetLocalSize(islist[i], &n));
505: if (n) {
506: PetscCall(ISGetIndices(islist[i], &iidx));
507: PetscCall(PetscArraycpy(idx + N, iidx, n));
508: PetscCall(ISRestoreIndices(islist[i], &iidx));
509: N += n;
510: }
511: }
512: }
513: PetscCall(ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: /*@
518: ISListToPair - Convert an `IS` list to a pair of `IS` of equal length defining an equivalent integer multimap.
519: Each `IS` on the input list is assigned an integer j so that all of the indices of that `IS` are
520: mapped to j.
522: Collective
524: Input Parameters:
525: + comm - `MPI_Comm`
526: . listlen - `IS` list length
527: - islist - `IS` list
529: Output Parameters:
530: + xis - domain `IS`
531: - yis - range `IS`
533: Level: developer
535: Notes:
536: The global integers assigned to the `IS` of the local input list might not correspond to the
537: local numbers of the `IS` on that list, but the two *orderings* are the same: the global
538: integers assigned to the `IS` on the local list form a strictly increasing sequence.
540: The `IS` on the input list can belong to subcommunicators of comm, and the subcommunicators
541: on the input `IS` list are assumed to be in a "deadlock-free" order.
543: Local lists of `PetscObject`s (or their subcomms) on a comm are "deadlock-free" if subcomm1
544: precedes subcomm2 on any local list, then it precedes subcomm2 on all ranks.
545: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
546: numbering. This is ensured, for example, by `ISPairToList()`.
548: .seealso: `IS`, `ISPairToList()`
549: @*/
550: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
551: {
552: PetscInt ncolors, *colors, i, leni, len, *xinds, *yinds, k, j;
553: const PetscInt *indsi;
555: PetscFunctionBegin;
556: PetscCall(PetscMalloc1(listlen, &colors));
557: PetscCall(PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject *)islist, &ncolors, colors));
558: len = 0;
559: for (i = 0; i < listlen; ++i) {
560: PetscCall(ISGetLocalSize(islist[i], &leni));
561: len += leni;
562: }
563: PetscCall(PetscMalloc1(len, &xinds));
564: PetscCall(PetscMalloc1(len, &yinds));
565: k = 0;
566: for (i = 0; i < listlen; ++i) {
567: PetscCall(ISGetLocalSize(islist[i], &leni));
568: PetscCall(ISGetIndices(islist[i], &indsi));
569: for (j = 0; j < leni; ++j) {
570: xinds[k] = indsi[j];
571: yinds[k] = colors[i];
572: ++k;
573: }
574: }
575: PetscCall(PetscFree(colors));
576: PetscCall(ISCreateGeneral(comm, len, xinds, PETSC_OWN_POINTER, xis));
577: PetscCall(ISCreateGeneral(comm, len, yinds, PETSC_OWN_POINTER, yis));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@
582: ISPairToList - Convert an `IS` pair encoding an integer map to a list of `IS`.
584: Collective
586: Input Parameters:
587: + xis - domain `IS`
588: - yis - range `IS`
590: Output Parameters:
591: + listlen - length of `islist`
592: - islist - list of `IS`s breaking up indis by color
594: Level: developer
596: Notes:
597: Each `IS` on the output list contains the preimage for each index on the second input
598: `IS`. The `IS` on the output list are constructed on the subcommunicators of the input `IS`
599: pair. Each subcommunicator corresponds to the preimage of some index j -- this subcomm
600: contains exactly the MPI processes that assign some indices i to j. This is essentially the inverse
601: of `ISListToPair()`.
603: `xis` and `yis` must be of the same length and have congruent communicators.
605: The resulting `IS` have subcommunicators in a "deadlock-free" order (see `ISListToPair()`).
607: .seealso: `IS`, `ISListToPair()`
608: @*/
609: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
610: {
611: IS indis = xis, coloris = yis;
612: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount, l;
613: PetscMPIInt rank, size, llow, lhigh, low, high, color, subsize;
614: const PetscInt *ccolors, *cinds;
615: MPI_Comm comm, subcomm;
617: PetscFunctionBegin;
620: PetscCheckSameComm(xis, 1, yis, 2);
621: PetscAssertPointer(listlen, 3);
622: PetscAssertPointer(islist, 4);
623: PetscCall(PetscObjectGetComm((PetscObject)xis, &comm));
624: PetscCallMPI(MPI_Comm_rank(comm, &rank));
625: PetscCallMPI(MPI_Comm_rank(comm, &size));
626: /* Extract, copy and sort the local indices and colors on the color. */
627: PetscCall(ISGetLocalSize(coloris, &llen));
628: PetscCall(ISGetLocalSize(indis, &ilen));
629: PetscCheck(llen == ilen, comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %" PetscInt_FMT " and %" PetscInt_FMT, ilen, llen);
630: PetscCall(ISGetIndices(coloris, &ccolors));
631: PetscCall(ISGetIndices(indis, &cinds));
632: PetscCall(PetscMalloc2(ilen, &inds, llen, &colors));
633: PetscCall(PetscArraycpy(inds, cinds, ilen));
634: PetscCall(PetscArraycpy(colors, ccolors, llen));
635: PetscCall(PetscSortIntWithArray(llen, colors, inds));
636: /* Determine the global extent of colors. */
637: llow = 0;
638: lhigh = -1;
639: lstart = 0;
640: lcount = 0;
641: while (lstart < llen) {
642: lend = lstart + 1;
643: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
644: llow = PetscMin(llow, colors[lstart]);
645: lhigh = PetscMax(lhigh, colors[lstart]);
646: ++lcount;
647: }
648: PetscCall(MPIU_Allreduce(&llow, &low, 1, MPI_INT, MPI_MIN, comm));
649: PetscCall(MPIU_Allreduce(&lhigh, &high, 1, MPI_INT, MPI_MAX, comm));
650: *listlen = 0;
651: if (low <= high) {
652: if (lcount > 0) {
653: *listlen = lcount;
654: if (!*islist) PetscCall(PetscMalloc1(lcount, islist));
655: }
656: /*
657: Traverse all possible global colors, and participate in the subcommunicators
658: for the locally-supported colors.
659: */
660: lcount = 0;
661: lstart = 0;
662: lend = 0;
663: for (l = low; l <= high; ++l) {
664: /*
665: Find the range of indices with the same color, which is not smaller than l.
666: Observe that, since colors is sorted, and is a subsequence of [low,high],
667: as soon as we find a new color, it is >= l.
668: */
669: if (lstart < llen) {
670: /* The start of the next locally-owned color is identified. Now look for the end. */
671: if (lstart == lend) {
672: lend = lstart + 1;
673: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
674: }
675: /* Now check whether the identified color segment matches l. */
676: PetscCheck(colors[lstart] >= l, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Locally owned color %" PetscInt_FMT " at location %" PetscInt_FMT " is < than the next global color %" PetscInt_FMT, colors[lstart], lcount, l);
677: }
678: color = (PetscMPIInt)(colors[lstart] == l);
679: /* Check whether a proper subcommunicator exists. */
680: PetscCall(MPIU_Allreduce(&color, &subsize, 1, MPI_INT, MPI_SUM, comm));
682: if (subsize == 1) subcomm = PETSC_COMM_SELF;
683: else if (subsize == size) subcomm = comm;
684: else {
685: /* a proper communicator is necessary, so we create it. */
686: PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm));
687: }
688: if (colors[lstart] == l) {
689: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
690: PetscCall(ISCreateGeneral(subcomm, lend - lstart, inds + lstart, PETSC_COPY_VALUES, *islist + lcount));
691: /* Position lstart at the beginning of the next local color. */
692: lstart = lend;
693: /* Increment the counter of the local colors split off into an IS. */
694: ++lcount;
695: }
696: if (subsize > 0 && subsize < size) {
697: /*
698: Irrespective of color, destroy the split off subcomm:
699: a subcomm used in the IS creation above is duplicated
700: into a proper PETSc comm.
701: */
702: PetscCallMPI(MPI_Comm_free(&subcomm));
703: }
704: } /* for (l = low; l < high; ++l) */
705: } /* if (low <= high) */
706: PetscCall(PetscFree2(inds, colors));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: ISEmbed - Embed `IS` `a` into `IS` `b` by finding the locations in `b` that have the same indices as in `a`.
712: If `c` is the `IS` of these locations, we have `a = b*c`, regarded as a composition of the
713: corresponding `ISLocalToGlobalMapping`.
715: Not Collective
717: Input Parameters:
718: + a - `IS` to embed
719: . b - `IS` to embed into
720: - drop - flag indicating whether to drop indices of `a` that are not in `b`.
722: Output Parameter:
723: . c - local embedding indices
725: Level: developer
727: Notes:
728: If some of the global indices of `a` are not among the indices of `b`, the embedding is impossible. The local indices of `a`
729: corresponding to these global indices are either mapped to -1 (if `!drop`) or are omitted (if `drop`). In the former
730: case the size of `c` is the same as that of `a`, in the latter case the size of `c` may be smaller.
732: The resulting `IS` is sequential, since the index substitution it encodes is purely local.
734: .seealso: `IS`, `ISLocalToGlobalMapping`
735: @*/
736: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
737: {
738: ISLocalToGlobalMapping ltog;
739: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
740: PetscInt alen, clen, *cindices, *cindices2;
741: const PetscInt *aindices;
743: PetscFunctionBegin;
746: PetscAssertPointer(c, 4);
747: PetscCall(ISLocalToGlobalMappingCreateIS(b, <og));
748: PetscCall(ISGetLocalSize(a, &alen));
749: PetscCall(ISGetIndices(a, &aindices));
750: PetscCall(PetscMalloc1(alen, &cindices));
751: if (!drop) gtoltype = IS_GTOLM_MASK;
752: PetscCall(ISGlobalToLocalMappingApply(ltog, gtoltype, alen, aindices, &clen, cindices));
753: PetscCall(ISRestoreIndices(a, &aindices));
754: PetscCall(ISLocalToGlobalMappingDestroy(<og));
755: if (clen != alen) {
756: cindices2 = cindices;
757: PetscCall(PetscMalloc1(clen, &cindices));
758: PetscCall(PetscArraycpy(cindices, cindices2, clen));
759: PetscCall(PetscFree(cindices2));
760: }
761: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c));
762: PetscFunctionReturn(PETSC_SUCCESS);
763: }
765: /*@
766: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
768: Not Collective
770: Input Parameters:
771: + f - `IS` to sort
772: - always - build the permutation even when `f`'s indices are nondecreasing.
774: Output Parameter:
775: . h - permutation or `NULL`, if `f` is nondecreasing and `always` == `PETSC_FALSE`.
777: Level: advanced
779: Notes:
780: Indices in `f` are unchanged. f[h[i]] is the i-th smallest f index.
782: If always == `PETSC_FALSE`, an extra check is performed to see whether
783: the `f` indices are nondecreasing. `h` is built on `PETSC_COMM_SELF`, since
784: the permutation has a local meaning only.
786: .seealso: `IS`, `ISLocalToGlobalMapping`, `ISSort()`
787: @*/
788: PetscErrorCode ISSortPermutation(IS f, PetscBool always, IS *h)
789: {
790: const PetscInt *findices;
791: PetscInt fsize, *hindices, i;
792: PetscBool isincreasing;
794: PetscFunctionBegin;
796: PetscAssertPointer(h, 3);
797: PetscCall(ISGetLocalSize(f, &fsize));
798: PetscCall(ISGetIndices(f, &findices));
799: *h = NULL;
800: if (!always) {
801: isincreasing = PETSC_TRUE;
802: for (i = 1; i < fsize; ++i) {
803: if (findices[i] <= findices[i - 1]) {
804: isincreasing = PETSC_FALSE;
805: break;
806: }
807: }
808: if (isincreasing) {
809: PetscCall(ISRestoreIndices(f, &findices));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
812: }
813: PetscCall(PetscMalloc1(fsize, &hindices));
814: for (i = 0; i < fsize; ++i) hindices[i] = i;
815: PetscCall(PetscSortIntWithPermutation(fsize, findices, hindices));
816: PetscCall(ISRestoreIndices(f, &findices));
817: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, fsize, hindices, PETSC_OWN_POINTER, h));
818: PetscCall(ISSetInfo(*h, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, PETSC_TRUE));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }