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