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