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_INT_MAX;
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: PetscBool f;
123: const PetscInt *i1, *i2;
124: PetscInt n1, n2, n3, p1, p2, *iout;
126: PetscFunctionBegin;
129: PetscCheckSameComm(is1, 1, is2, 2);
131: PetscCall(ISSorted(is1, &f));
132: PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 1 is not sorted");
133: PetscCall(ISSorted(is2, &f));
134: PetscCheck(f, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Arg 2 is not sorted");
136: PetscCall(ISGetLocalSize(is1, &n1));
137: PetscCall(ISGetLocalSize(is2, &n2));
138: if (!n2) {
139: PetscCall(ISDuplicate(is1, is3));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
142: PetscCall(ISGetIndices(is1, &i1));
143: PetscCall(ISGetIndices(is2, &i2));
145: p1 = 0;
146: p2 = 0;
147: n3 = 0;
148: do {
149: if (p1 == n1) { /* cleanup for is2 */
150: n3 += n2 - p2;
151: break;
152: } else {
153: while (p2 < n2 && i2[p2] < i1[p1]) {
154: n3++;
155: p2++;
156: }
157: if (p2 == n2) {
158: /* cleanup for is1 */
159: n3 += n1 - p1;
160: break;
161: } else {
162: if (i2[p2] == i1[p1]) {
163: n3++;
164: p1++;
165: p2++;
166: }
167: }
168: }
169: if (p2 == n2) {
170: /* cleanup for is1 */
171: n3 += n1 - p1;
172: break;
173: } else {
174: while (p1 < n1 && i1[p1] < i2[p2]) {
175: n3++;
176: p1++;
177: }
178: if (p1 == n1) {
179: /* clean up for is2 */
180: n3 += n2 - p2;
181: break;
182: } else {
183: if (i1[p1] == i2[p2]) {
184: n3++;
185: p1++;
186: p2++;
187: }
188: }
189: }
190: } while (p1 < n1 || p2 < n2);
192: if (n3 == n1) { /* no new elements to be added */
193: PetscCall(ISRestoreIndices(is1, &i1));
194: PetscCall(ISRestoreIndices(is2, &i2));
195: PetscCall(ISDuplicate(is1, is3));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
198: PetscCall(PetscMalloc1(n3, &iout));
200: p1 = 0;
201: p2 = 0;
202: n3 = 0;
203: do {
204: if (p1 == n1) { /* cleanup for is2 */
205: while (p2 < n2) iout[n3++] = i2[p2++];
206: break;
207: } else {
208: while (p2 < n2 && i2[p2] < i1[p1]) iout[n3++] = i2[p2++];
209: if (p2 == n2) { /* cleanup for is1 */
210: while (p1 < n1) iout[n3++] = i1[p1++];
211: break;
212: } else {
213: if (i2[p2] == i1[p1]) {
214: iout[n3++] = i1[p1++];
215: p2++;
216: }
217: }
218: }
219: if (p2 == n2) { /* cleanup for is1 */
220: while (p1 < n1) iout[n3++] = i1[p1++];
221: break;
222: } else {
223: while (p1 < n1 && i1[p1] < i2[p2]) iout[n3++] = i1[p1++];
224: if (p1 == n1) { /* clean up for is2 */
225: while (p2 < n2) iout[n3++] = i2[p2++];
226: break;
227: } else {
228: if (i1[p1] == i2[p2]) {
229: iout[n3++] = i1[p1++];
230: p2++;
231: }
232: }
233: }
234: } while (p1 < n1 || p2 < n2);
236: PetscCall(ISRestoreIndices(is1, &i1));
237: PetscCall(ISRestoreIndices(is2, &i2));
238: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is1), n3, iout, PETSC_OWN_POINTER, is3));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@
243: ISExpand - Computes the union of two index sets, by concatenating 2 lists and
244: removing duplicates.
246: Collective
248: Input Parameters:
249: + is1 - first index set
250: - is2 - index values to be added
252: Output Parameter:
253: . isout - `is1` + `is2` The index set `is2` is appended to `is1` removing duplicates
255: Level: intermediate
257: Notes:
258: Negative values are removed from the lists. This requires O(imax-imin)
259: memory and O(imax-imin) work, where imin and imax are the bounds on the
260: indices in `is1` and `is2`.
262: `is1` and `is2` do not need to be sorted.
264: The operations are performed separately on each MPI rank
266: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISIntersect()`
267: @*/
268: PetscErrorCode ISExpand(IS is1, IS is2, IS *isout)
269: {
270: PetscInt i, n1, n2, imin, imax, nout, *iout;
271: const PetscInt *i1, *i2;
272: PetscBool sorted1 = PETSC_TRUE, sorted2 = PETSC_TRUE;
273: PetscBT mask;
274: MPI_Comm comm;
276: PetscFunctionBegin;
279: PetscAssertPointer(isout, 3);
281: PetscCheck(is1 || is2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both arguments cannot be NULL");
282: if (!is1) {
283: PetscCall(ISDuplicate(is2, isout));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
286: if (!is2) {
287: PetscCall(ISDuplicate(is1, isout));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
290: PetscCall(ISGetIndices(is1, &i1));
291: PetscCall(ISGetLocalSize(is1, &n1));
292: PetscCall(ISGetIndices(is2, &i2));
293: PetscCall(ISGetLocalSize(is2, &n2));
295: /* Create a bit mask array to contain required values */
296: if (n1 || n2) {
297: imin = PETSC_INT_MAX;
298: imax = 0;
299: if (n1) {
300: PetscCall(ISSorted(is1, &sorted1));
301: if (sorted1 && i1[0] >= 0) imin = i1[0], imax = i1[n1 - 1];
302: else
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: }
309: if (n2) {
310: PetscCall(ISSorted(is2, &sorted2));
311: if (sorted2 && i2[0] >= 0) imin = PetscMin(imin, i2[0]), imax = PetscMax(imax, i2[n2 - 1]);
312: else
313: for (i = 0; i < n2; i++) {
314: if (i2[i] < 0) continue;
315: imin = PetscMin(imin, i2[i]);
316: imax = PetscMax(imax, i2[i]);
317: }
318: }
319: } else imin = imax = 0;
321: PetscCall(PetscMalloc1(n1 + n2, &iout));
322: nout = 0;
323: PetscCall(PetscBTCreate(imax - imin, &mask));
324: /* Put the values from is1 */
325: for (i = 0; i < n1; i++) {
326: if (i1[i] < 0) continue;
327: if (!PetscBTLookupSet(mask, i1[i] - imin)) iout[nout++] = i1[i];
328: }
329: n1 = -nout;
330: PetscCall(ISRestoreIndices(is1, &i1));
331: /* Put the values from is2 */
332: for (i = 0; i < n2; i++) {
333: if (i2[i] < 0) continue;
334: if (!PetscBTLookupSet(mask, i2[i] - imin)) iout[nout++] = i2[i];
335: }
336: PetscCall(ISRestoreIndices(is2, &i2));
338: /* create the new IS containing the sum */
339: PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
340: PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
341: /* no entries of is2 (resp. is1) was inserted, so if is1 (resp. is2) is sorted, then so is isout */
342: if ((-n1 == nout && sorted1) || (n1 == 0 && sorted2)) PetscCall(ISSetInfo(*isout, IS_SORTED, IS_LOCAL, PETSC_FALSE, PETSC_TRUE));
343: PetscCall(PetscBTDestroy(&mask));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: ISIntersect - Computes the intersection of two index sets, by sorting and comparing.
350: Collective
352: Input Parameters:
353: + is1 - first index set
354: - is2 - second index set
356: Output Parameter:
357: . isout - the sorted intersection of `is1` and `is2`
359: Level: intermediate
361: Notes:
362: Negative values are removed from the lists. This requires O(min(is1,is2))
363: memory and O(max(is1,is2)log(max(is1,is2))) work
365: `is1` and `is2` do not need to be sorted.
367: The operations are performed separately on each MPI rank
369: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISConcatenate()`
370: @*/
371: PetscErrorCode ISIntersect(IS is1, IS is2, IS *isout)
372: {
373: PetscInt i, n1, n2, nout, *iout;
374: const PetscInt *i1, *i2;
375: IS is1sorted = NULL, is2sorted = NULL;
376: PetscBool sorted, lsorted;
377: MPI_Comm comm;
379: PetscFunctionBegin;
382: PetscCheckSameComm(is1, 1, is2, 2);
383: PetscAssertPointer(isout, 3);
384: PetscCall(PetscObjectGetComm((PetscObject)is1, &comm));
386: PetscCall(ISGetLocalSize(is1, &n1));
387: PetscCall(ISGetLocalSize(is2, &n2));
388: if (n1 < n2) {
389: IS tempis = is1;
390: PetscInt ntemp = n1;
392: is1 = is2;
393: is2 = tempis;
394: n1 = n2;
395: n2 = ntemp;
396: }
397: PetscCall(ISSorted(is1, &lsorted));
398: PetscCallMPI(MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm));
399: if (!sorted) {
400: PetscCall(ISDuplicate(is1, &is1sorted));
401: PetscCall(ISSort(is1sorted));
402: PetscCall(ISGetIndices(is1sorted, &i1));
403: } else {
404: is1sorted = is1;
405: PetscCall(PetscObjectReference((PetscObject)is1));
406: PetscCall(ISGetIndices(is1, &i1));
407: }
408: PetscCall(ISSorted(is2, &lsorted));
409: PetscCallMPI(MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm));
410: if (!sorted) {
411: PetscCall(ISDuplicate(is2, &is2sorted));
412: PetscCall(ISSort(is2sorted));
413: PetscCall(ISGetIndices(is2sorted, &i2));
414: } else {
415: is2sorted = is2;
416: PetscCall(PetscObjectReference((PetscObject)is2));
417: PetscCall(ISGetIndices(is2, &i2));
418: }
420: PetscCall(PetscMalloc1(n2, &iout));
422: for (i = 0, nout = 0; i < n2; i++) {
423: PetscInt key = i2[i];
424: PetscInt loc;
426: PetscCall(ISLocate(is1sorted, key, &loc));
427: if (loc >= 0) {
428: if (!nout || iout[nout - 1] < key) iout[nout++] = key;
429: }
430: }
431: PetscCall(PetscRealloc(nout * sizeof(PetscInt), &iout));
433: /* create the new IS containing the sum */
434: PetscCall(ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout));
435: PetscCall(ISSetInfo(*isout, IS_SORTED, IS_GLOBAL, PETSC_FALSE, PETSC_TRUE));
436: PetscCall(ISRestoreIndices(is2sorted, &i2));
437: PetscCall(ISDestroy(&is2sorted));
438: PetscCall(ISRestoreIndices(is1sorted, &i1));
439: PetscCall(ISDestroy(&is1sorted));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
444: {
445: PetscFunctionBegin;
446: *isect = NULL;
447: if (is2 && is1) {
448: char composeStr[33] = {0};
449: PetscObjectId is2id;
451: PetscCall(PetscObjectGetId((PetscObject)is2, &is2id));
452: PetscCall(PetscSNPrintf(composeStr, 32, "ISIntersect_Caching_%" PetscInt64_FMT, is2id));
453: PetscCall(PetscObjectQuery((PetscObject)is1, composeStr, (PetscObject *)isect));
454: if (*isect == NULL) {
455: PetscCall(ISIntersect(is1, is2, isect));
456: PetscCall(PetscObjectCompose((PetscObject)is1, composeStr, (PetscObject)*isect));
457: } else {
458: PetscCall(PetscObjectReference((PetscObject)*isect));
459: }
460: }
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: /*@
465: ISConcatenate - Forms a new `IS` by locally concatenating the indices from an `IS` list without reordering.
467: Collective
469: Input Parameters:
470: + comm - communicator of the concatenated `IS`.
471: . len - size of islist array (nonnegative)
472: - islist - array of index sets
474: Output Parameter:
475: . isout - The concatenated index set; empty, if `len` == 0.
477: Level: intermediate
479: Notes:
480: The semantics of calling this on comm imply that the comms of the members of `islist` also contain this rank.
482: .seealso: [](sec_scatter), `IS`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISIntersect()`
483: @*/
484: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
485: {
486: PetscInt i, n, N;
487: const PetscInt *iidx;
488: PetscInt *idx;
490: PetscFunctionBegin;
491: if (len) PetscAssertPointer(islist, 3);
492: if (PetscDefined(USE_DEBUG)) {
493: for (i = 0; i < len; ++i)
495: }
496: PetscAssertPointer(isout, 4);
497: if (!len) {
498: PetscCall(ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
501: PetscCheck(len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Negative array length: %" PetscInt_FMT, len);
502: N = 0;
503: for (i = 0; i < len; ++i) {
504: if (islist[i]) {
505: PetscCall(ISGetLocalSize(islist[i], &n));
506: N += n;
507: }
508: }
509: PetscCall(PetscMalloc1(N, &idx));
510: N = 0;
511: for (i = 0; i < len; ++i) {
512: if (islist[i]) {
513: PetscCall(ISGetLocalSize(islist[i], &n));
514: if (n) {
515: PetscCall(ISGetIndices(islist[i], &iidx));
516: PetscCall(PetscArraycpy(idx + N, iidx, n));
517: PetscCall(ISRestoreIndices(islist[i], &iidx));
518: N += n;
519: }
520: }
521: }
522: PetscCall(ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@
527: ISListToPair - Convert an `IS` list to a pair of `IS` of equal length defining an equivalent integer multimap.
528: Each `IS` in `islist` is assigned an integer j so that all of the indices of that `IS` are
529: mapped to j.
531: Collective
533: Input Parameters:
534: + comm - `MPI_Comm`
535: . listlen - `IS` list length
536: - islist - `IS` list
538: Output Parameters:
539: + xis - domain `IS`
540: - yis - range `IS`
542: Level: developer
544: Notes:
545: The global integers assigned to the `IS` of `islist` might not correspond to the
546: local numbers of the `IS` on that list, but the two *orderings* are the same: the global
547: integers assigned to the `IS` on the local list form a strictly increasing sequence.
549: The `IS` in `islist` can belong to subcommunicators of `comm`, and the subcommunicators
550: on the input `IS` list are assumed to be in a "deadlock-free" order.
552: Local lists of `PetscObject`s (or their subcomms) on a comm are "deadlock-free" if subcomm1
553: precedes subcomm2 on any local list, then it precedes subcomm2 on all ranks.
554: Equivalently, the local numbers of the subcomms on each local list are drawn from some global
555: numbering. This is ensured, for example, by `ISPairToList()`.
557: .seealso: `IS`, `ISPairToList()`
558: @*/
559: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
560: {
561: PetscInt ncolors, *colors, leni, len, *xinds, *yinds, k;
562: const PetscInt *indsi;
564: PetscFunctionBegin;
565: PetscCall(PetscMalloc1(listlen, &colors));
566: PetscCall(PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject *)islist, &ncolors, colors));
567: len = 0;
568: for (PetscInt i = 0; i < listlen; ++i) {
569: PetscCall(ISGetLocalSize(islist[i], &leni));
570: len += leni;
571: }
572: PetscCall(PetscMalloc1(len, &xinds));
573: PetscCall(PetscMalloc1(len, &yinds));
574: k = 0;
575: for (PetscInt i = 0; i < listlen; ++i) {
576: PetscCall(ISGetLocalSize(islist[i], &leni));
577: PetscCall(ISGetIndices(islist[i], &indsi));
578: for (PetscInt j = 0; j < leni; ++j) {
579: xinds[k] = indsi[j];
580: yinds[k] = colors[i];
581: ++k;
582: }
583: }
584: PetscCall(PetscFree(colors));
585: PetscCall(ISCreateGeneral(comm, len, xinds, PETSC_OWN_POINTER, xis));
586: PetscCall(ISCreateGeneral(comm, len, yinds, PETSC_OWN_POINTER, yis));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /*@C
591: ISPairToList - Convert an `IS` pair encoding an integer map to a list of `IS`.
593: Collective
595: Input Parameters:
596: + xis - domain `IS`
597: - yis - range `IS`, the maximum value must be less than `PETSC_MPI_INT_MAX`
599: Output Parameters:
600: + listlen - length of `islist`
601: - islist - list of `IS`s breaking up indices by color
603: Level: developer
605: Notes:
606: Each `IS` in `islist` contains the preimage for each index on `yis`.
607: The `IS` in `islist` are constructed on the subcommunicators of the input `IS`
608: pair. Each subcommunicator corresponds to the preimage of some index j -- this subcomm
609: contains exactly the MPI processes that assign some indices i to j. This is essentially the inverse
610: of `ISListToPair()`.
612: `xis` and `yis` must be of the same length and have congruent communicators.
614: The resulting `IS` have subcommunicators in a "deadlock-free" order (see `ISListToPair()`).
616: .seealso: `IS`, `ISListToPair()`
617: @*/
618: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
619: {
620: IS indis = xis, coloris = yis;
621: PetscInt *inds, *colors, llen, ilen, lstart, lend, lcount, l;
622: PetscMPIInt rank, size, llow, lhigh, low, high, color, subsize;
623: const PetscInt *ccolors, *cinds;
624: MPI_Comm comm, subcomm;
626: PetscFunctionBegin;
629: PetscCheckSameComm(xis, 1, yis, 2);
630: PetscAssertPointer(listlen, 3);
631: PetscAssertPointer(islist, 4);
632: PetscCall(PetscObjectGetComm((PetscObject)xis, &comm));
633: PetscCallMPI(MPI_Comm_rank(comm, &rank));
634: PetscCallMPI(MPI_Comm_rank(comm, &size));
635: /* Extract, copy and sort the local indices and colors on the color. */
636: PetscCall(ISGetLocalSize(coloris, &llen));
637: PetscCall(ISGetLocalSize(indis, &ilen));
638: PetscCheck(llen == ilen, comm, PETSC_ERR_ARG_SIZ, "Incompatible IS sizes: %" PetscInt_FMT " and %" PetscInt_FMT, ilen, llen);
639: PetscCall(ISGetIndices(coloris, &ccolors));
640: PetscCall(ISGetIndices(indis, &cinds));
641: PetscCall(PetscMalloc2(ilen, &inds, llen, &colors));
642: PetscCall(PetscArraycpy(inds, cinds, ilen));
643: PetscCall(PetscArraycpy(colors, ccolors, llen));
644: PetscCall(PetscSortIntWithArray(llen, colors, inds));
645: /* Determine the global extent of colors. */
646: llow = 0;
647: lhigh = -1;
648: lstart = 0;
649: lcount = 0;
650: while (lstart < llen) {
651: lend = lstart + 1;
652: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
653: PetscCall(PetscMPIIntCast(PetscMin(llow, colors[lstart]), &llow));
654: PetscCall(PetscMPIIntCast(PetscMax(lhigh, colors[lstart]), &lhigh));
655: ++lcount;
656: }
657: PetscCallMPI(MPIU_Allreduce(&llow, &low, 1, MPI_INT, MPI_MIN, comm));
658: PetscCallMPI(MPIU_Allreduce(&lhigh, &high, 1, MPI_INT, MPI_MAX, comm));
659: *listlen = 0;
660: if (low <= high) {
661: if (lcount > 0) {
662: *listlen = lcount;
663: if (!*islist) PetscCall(PetscMalloc1(lcount, islist));
664: }
665: /*
666: Traverse all possible global colors, and participate in the subcommunicators
667: for the locally-supported colors.
668: */
669: lcount = 0;
670: lstart = 0;
671: lend = 0;
672: for (l = low; l <= high; ++l) {
673: /*
674: Find the range of indices with the same color, which is not smaller than l.
675: Observe that, since colors is sorted, and is a subsequence of [low,high],
676: as soon as we find a new color, it is >= l.
677: */
678: if (lstart < llen) {
679: /* The start of the next locally-owned color is identified. Now look for the end. */
680: if (lstart == lend) {
681: lend = lstart + 1;
682: while (lend < llen && colors[lend] == colors[lstart]) ++lend;
683: }
684: /* Now check whether the identified color segment matches l. */
685: 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);
686: }
687: color = (PetscMPIInt)(colors[lstart] == l);
688: /* Check whether a proper subcommunicator exists. */
689: PetscCallMPI(MPIU_Allreduce(&color, &subsize, 1, MPI_INT, MPI_SUM, comm));
691: if (subsize == 1) subcomm = PETSC_COMM_SELF;
692: else if (subsize == size) subcomm = comm;
693: else {
694: /* a proper communicator is necessary, so we create it. */
695: PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm));
696: }
697: if (colors[lstart] == l) {
698: /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
699: PetscCall(ISCreateGeneral(subcomm, lend - lstart, inds + lstart, PETSC_COPY_VALUES, *islist + lcount));
700: /* Position lstart at the beginning of the next local color. */
701: lstart = lend;
702: /* Increment the counter of the local colors split off into an IS. */
703: ++lcount;
704: }
705: if (subsize > 0 && subsize < size) {
706: /*
707: Irrespective of color, destroy the split off subcomm:
708: a subcomm used in the IS creation above is duplicated
709: into a proper PETSc comm.
710: */
711: PetscCallMPI(MPI_Comm_free(&subcomm));
712: }
713: } /* for (l = low; l < high; ++l) */
714: } /* if (low <= high) */
715: PetscCall(PetscFree2(inds, colors));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: ISEmbed - Embed `IS` `a` into `IS` `b` by finding the locations in `b` that have the same indices as in `a`.
721: If `c` is the `IS` of these locations, we have `a = b*c`, regarded as a composition of the
722: corresponding `ISLocalToGlobalMapping`.
724: Not Collective
726: Input Parameters:
727: + a - `IS` to embed
728: . b - `IS` to embed into
729: - drop - flag indicating whether to drop indices of `a` that are not in `b`.
731: Output Parameter:
732: . c - local embedding indices
734: Level: developer
736: Notes:
737: If some of the global indices of `a` are not among the indices of `b`, the embedding is impossible. The local indices of `a`
738: corresponding to these global indices are either mapped to -1 (if `!drop`) or are omitted (if `drop`). In the former
739: case the size of `c` is the same as that of `a`, in the latter case the size of `c` may be smaller.
741: The resulting `IS` is sequential, since the index substitution it encodes is purely local.
743: .seealso: `IS`, `ISLocalToGlobalMapping`
744: @*/
745: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
746: {
747: ISLocalToGlobalMapping ltog;
748: ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
749: PetscInt alen, clen, *cindices, *cindices2;
750: const PetscInt *aindices;
752: PetscFunctionBegin;
755: PetscAssertPointer(c, 4);
756: PetscCall(ISLocalToGlobalMappingCreateIS(b, <og));
757: PetscCall(ISGetLocalSize(a, &alen));
758: PetscCall(ISGetIndices(a, &aindices));
759: PetscCall(PetscMalloc1(alen, &cindices));
760: if (!drop) gtoltype = IS_GTOLM_MASK;
761: PetscCall(ISGlobalToLocalMappingApply(ltog, gtoltype, alen, aindices, &clen, cindices));
762: PetscCall(ISRestoreIndices(a, &aindices));
763: PetscCall(ISLocalToGlobalMappingDestroy(<og));
764: if (clen != alen) {
765: cindices2 = cindices;
766: PetscCall(PetscMalloc1(clen, &cindices));
767: PetscCall(PetscArraycpy(cindices, cindices2, clen));
768: PetscCall(PetscFree(cindices2));
769: }
770: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c));
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: /*@
775: ISSortPermutation - calculate the permutation of the indices into a nondecreasing order.
777: Not Collective
779: Input Parameters:
780: + f - `IS` to sort
781: - always - build the permutation even when `f`'s indices are nondecreasing.
783: Output Parameter:
784: . h - permutation or `NULL`, if `f` is nondecreasing and `always` == `PETSC_FALSE`.
786: Level: advanced
788: Notes:
789: Indices in `f` are unchanged. f[h[i]] is the i-th smallest f index.
791: If always == `PETSC_FALSE`, an extra check is performed to see whether
792: the `f` indices are nondecreasing. `h` is built on `PETSC_COMM_SELF`, since
793: the permutation has a local meaning only.
795: .seealso: `IS`, `ISLocalToGlobalMapping`, `ISSort()`
796: @*/
797: PetscErrorCode ISSortPermutation(IS f, PetscBool always, IS *h)
798: {
799: const PetscInt *findices;
800: PetscInt fsize, *hindices, i;
801: PetscBool isincreasing;
803: PetscFunctionBegin;
805: PetscAssertPointer(h, 3);
806: PetscCall(ISGetLocalSize(f, &fsize));
807: PetscCall(ISGetIndices(f, &findices));
808: *h = NULL;
809: if (!always) {
810: isincreasing = PETSC_TRUE;
811: for (i = 1; i < fsize; ++i) {
812: if (findices[i] <= findices[i - 1]) {
813: isincreasing = PETSC_FALSE;
814: break;
815: }
816: }
817: if (isincreasing) {
818: PetscCall(ISRestoreIndices(f, &findices));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
821: }
822: PetscCall(PetscMalloc1(fsize, &hindices));
823: for (i = 0; i < fsize; ++i) hindices[i] = i;
824: PetscCall(PetscSortIntWithPermutation(fsize, findices, hindices));
825: PetscCall(ISRestoreIndices(f, &findices));
826: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, fsize, hindices, PETSC_OWN_POINTER, h));
827: PetscCall(ISSetInfo(*h, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, PETSC_TRUE));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }