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, &ltog));
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(&ltog));
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: }