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 on is1

 11:    Input Parameters:
 12: +  is1 - first index, to have items removed from it
 13: -  is2 - index values to be removed

 15:    Output Parameters:
 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. This requires O(imax-imin) memory and O(imax-imin)
 23:    work, where imin and imax are the bounds on the indices in is1.

 25:    If is2 is NULL, the result is the same as for an empty IS, i.e., a duplicate of is1.

 27:    The difference is computed separately on each MPI rank

 29: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISSum()`, `ISExpand()`
 30: @*/
 31: PetscErrorCode ISDifference(IS is1, IS is2, IS *isout)
 32: {
 33:   PetscInt        i, n1, n2, imin, imax, nout, *iout;
 34:   const PetscInt *i1, *i2;
 35:   PetscBT         mask;
 36:   MPI_Comm        comm;

 40:   if (!is2) {
 41:     ISDuplicate(is1, isout);
 42:     return 0;
 43:   }

 46:   ISGetIndices(is1, &i1);
 47:   ISGetLocalSize(is1, &n1);

 49:   /* Create a bit mask array to contain required values */
 50:   if (n1) {
 51:     imin = PETSC_MAX_INT;
 52:     imax = 0;
 53:     for (i = 0; i < n1; i++) {
 54:       if (i1[i] < 0) continue;
 55:       imin = PetscMin(imin, i1[i]);
 56:       imax = PetscMax(imax, i1[i]);
 57:     }
 58:   } else imin = imax = 0;

 60:   PetscBTCreate(imax - imin, &mask);
 61:   /* Put the values from is1 */
 62:   for (i = 0; i < n1; i++) {
 63:     if (i1[i] < 0) continue;
 64:     PetscBTSet(mask, i1[i] - imin);
 65:   }
 66:   ISRestoreIndices(is1, &i1);
 67:   /* Remove the values from is2 */
 68:   ISGetIndices(is2, &i2);
 69:   ISGetLocalSize(is2, &n2);
 70:   for (i = 0; i < n2; i++) {
 71:     if (i2[i] < imin || i2[i] > imax) continue;
 72:     PetscBTClear(mask, i2[i] - imin);
 73:   }
 74:   ISRestoreIndices(is2, &i2);

 76:   /* Count the number in the difference */
 77:   nout = 0;
 78:   for (i = 0; i < imax - imin + 1; i++) {
 79:     if (PetscBTLookup(mask, i)) nout++;
 80:   }

 82:   /* create the new IS containing the difference */
 83:   PetscMalloc1(nout, &iout);
 84:   nout = 0;
 85:   for (i = 0; i < imax - imin + 1; i++) {
 86:     if (PetscBTLookup(mask, i)) iout[nout++] = i + imin;
 87:   }
 88:   PetscObjectGetComm((PetscObject)is1, &comm);
 89:   ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);

 91:   PetscBTDestroy(&mask);
 92:   return 0;
 93: }

 95: /*@
 96:    ISSum - Computes the sum (union) of two index sets.

 98:    Only sequential version (at the moment)

100:    Input Parameters:
101: +  is1 - index set to be extended
102: -  is2 - index values to be added

104:    Output Parameter:
105: .   is3 - the sum; this can not be is1 or is2

107:    Level: intermediate

109:    Notes:
110:    If n1 and n2 are the sizes of the sets, this takes O(n1+n2) time;

112:    Both index sets need to be sorted on input.

114:    The sum is computed separately on each MPI rank

116: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISExpand()`
117: @*/
118: PetscErrorCode ISSum(IS is1, IS is2, IS *is3)
119: {
120:   MPI_Comm        comm;
121:   PetscBool       f;
122:   PetscMPIInt     size;
123:   const PetscInt *i1, *i2;
124:   PetscInt        n1, n2, n3, p1, p2, *iout;

128:   PetscObjectGetComm((PetscObject)(is1), &comm);
129:   MPI_Comm_size(comm, &size);

132:   ISSorted(is1, &f);
134:   ISSorted(is2, &f);

137:   ISGetLocalSize(is1, &n1);
138:   ISGetLocalSize(is2, &n2);
139:   if (!n2) {
140:     ISDuplicate(is1, is3);
141:     return 0;
142:   }
143:   ISGetIndices(is1, &i1);
144:   ISGetIndices(is2, &i2);

146:   p1 = 0;
147:   p2 = 0;
148:   n3 = 0;
149:   do {
150:     if (p1 == n1) { /* cleanup for is2 */
151:       n3 += n2 - p2;
152:       break;
153:     } else {
154:       while (p2 < n2 && i2[p2] < i1[p1]) {
155:         n3++;
156:         p2++;
157:       }
158:       if (p2 == n2) {
159:         /* cleanup for is1 */
160:         n3 += n1 - p1;
161:         break;
162:       } else {
163:         if (i2[p2] == i1[p1]) {
164:           n3++;
165:           p1++;
166:           p2++;
167:         }
168:       }
169:     }
170:     if (p2 == n2) {
171:       /* cleanup for is1 */
172:       n3 += n1 - p1;
173:       break;
174:     } else {
175:       while (p1 < n1 && i1[p1] < i2[p2]) {
176:         n3++;
177:         p1++;
178:       }
179:       if (p1 == n1) {
180:         /* clean up for is2 */
181:         n3 += n2 - p2;
182:         break;
183:       } else {
184:         if (i1[p1] == i2[p2]) {
185:           n3++;
186:           p1++;
187:           p2++;
188:         }
189:       }
190:     }
191:   } while (p1 < n1 || p2 < n2);

193:   if (n3 == n1) { /* no new elements to be added */
194:     ISRestoreIndices(is1, &i1);
195:     ISRestoreIndices(is2, &i2);
196:     ISDuplicate(is1, is3);
197:     return 0;
198:   }
199:   PetscMalloc1(n3, &iout);

201:   p1 = 0;
202:   p2 = 0;
203:   n3 = 0;
204:   do {
205:     if (p1 == n1) { /* cleanup for is2 */
206:       while (p2 < n2) iout[n3++] = i2[p2++];
207:       break;
208:     } else {
209:       while (p2 < n2 && i2[p2] < i1[p1]) iout[n3++] = i2[p2++];
210:       if (p2 == n2) { /* cleanup for is1 */
211:         while (p1 < n1) iout[n3++] = i1[p1++];
212:         break;
213:       } else {
214:         if (i2[p2] == i1[p1]) {
215:           iout[n3++] = i1[p1++];
216:           p2++;
217:         }
218:       }
219:     }
220:     if (p2 == n2) { /* cleanup for is1 */
221:       while (p1 < n1) iout[n3++] = i1[p1++];
222:       break;
223:     } else {
224:       while (p1 < n1 && i1[p1] < i2[p2]) iout[n3++] = i1[p1++];
225:       if (p1 == n1) { /* clean up for is2 */
226:         while (p2 < n2) iout[n3++] = i2[p2++];
227:         break;
228:       } else {
229:         if (i1[p1] == i2[p2]) {
230:           iout[n3++] = i1[p1++];
231:           p2++;
232:         }
233:       }
234:     }
235:   } while (p1 < n1 || p2 < n2);

237:   ISRestoreIndices(is1, &i1);
238:   ISRestoreIndices(is2, &i2);
239:   ISCreateGeneral(comm, n3, iout, PETSC_OWN_POINTER, is3);
240:   return 0;
241: }

243: /*@
244:    ISExpand - Computes the union of two index sets, by concatenating 2 lists and
245:    removing duplicates.

247:    Collective on is1

249:    Input Parameters:
250: +  is1 - first index set
251: -  is2 - index values to be added

253:    Output Parameters:
254: .  isout - is1 + is2 The index set is2 is appended to is1 removing duplicates

256:    Level: intermediate

258:    Notes:
259:    Negative values are removed from the lists. This requires O(imax-imin)
260:    memory and O(imax-imin) work, where imin and imax are the bounds on the
261:    indices in is1 and is2.

263:    The IS's do not need to be sorted.

265:    The operations are performed separately on each MPI rank

267: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISIntersect()`
268: @*/
269: PetscErrorCode ISExpand(IS is1, IS is2, IS *isout)
270: {
271:   PetscInt        i, n1, n2, imin, imax, nout, *iout;
272:   const PetscInt *i1, *i2;
273:   PetscBT         mask;
274:   MPI_Comm        comm;


281:   if (!is1) {
282:     ISDuplicate(is2, isout);
283:     return 0;
284:   }
285:   if (!is2) {
286:     ISDuplicate(is1, isout);
287:     return 0;
288:   }
289:   ISGetIndices(is1, &i1);
290:   ISGetLocalSize(is1, &n1);
291:   ISGetIndices(is2, &i2);
292:   ISGetLocalSize(is2, &n2);

294:   /* Create a bit mask array to contain required values */
295:   if (n1 || n2) {
296:     imin = PETSC_MAX_INT;
297:     imax = 0;
298:     for (i = 0; i < n1; i++) {
299:       if (i1[i] < 0) continue;
300:       imin = PetscMin(imin, i1[i]);
301:       imax = PetscMax(imax, i1[i]);
302:     }
303:     for (i = 0; i < n2; i++) {
304:       if (i2[i] < 0) continue;
305:       imin = PetscMin(imin, i2[i]);
306:       imax = PetscMax(imax, i2[i]);
307:     }
308:   } else imin = imax = 0;

310:   PetscMalloc1(n1 + n2, &iout);
311:   nout = 0;
312:   PetscBTCreate(imax - imin, &mask);
313:   /* Put the values from is1 */
314:   for (i = 0; i < n1; i++) {
315:     if (i1[i] < 0) continue;
316:     if (!PetscBTLookupSet(mask, i1[i] - imin)) iout[nout++] = i1[i];
317:   }
318:   ISRestoreIndices(is1, &i1);
319:   /* Put the values from is2 */
320:   for (i = 0; i < n2; i++) {
321:     if (i2[i] < 0) continue;
322:     if (!PetscBTLookupSet(mask, i2[i] - imin)) iout[nout++] = i2[i];
323:   }
324:   ISRestoreIndices(is2, &i2);

326:   /* create the new IS containing the sum */
327:   PetscObjectGetComm((PetscObject)is1, &comm);
328:   ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);

330:   PetscBTDestroy(&mask);
331:   return 0;
332: }

334: /*@
335:    ISIntersect - Computes the intersection of two index sets, by sorting and comparing.

337:    Collective on is1

339:    Input Parameters:
340: +  is1 - first index set
341: -  is2 - second index set

343:    Output Parameters:
344: .  isout - the sorted intersection of is1 and is2

346:    Level: intermediate

348:    Notes:
349:    Negative values are removed from the lists. This requires O(min(is1,is2))
350:    memory and O(max(is1,is2)log(max(is1,is2))) work

352:    The IS's do not need to be sorted.

354:    The operations are performed separately on each MPI rank

356: .seealso: [](sec_scatter), `IS`, `ISDestroy()`, `ISView()`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISConcatenate()`
357: @*/
358: PetscErrorCode ISIntersect(IS is1, IS is2, IS *isout)
359: {
360:   PetscInt        i, n1, n2, nout, *iout;
361:   const PetscInt *i1, *i2;
362:   IS              is1sorted = NULL, is2sorted = NULL;
363:   PetscBool       sorted, lsorted;
364:   MPI_Comm        comm;

370:   PetscObjectGetComm((PetscObject)is1, &comm);

372:   ISGetLocalSize(is1, &n1);
373:   ISGetLocalSize(is2, &n2);
374:   if (n1 < n2) {
375:     IS       tempis = is1;
376:     PetscInt ntemp  = n1;

378:     is1 = is2;
379:     is2 = tempis;
380:     n1  = n2;
381:     n2  = ntemp;
382:   }
383:   ISSorted(is1, &lsorted);
384:   MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm);
385:   if (!sorted) {
386:     ISDuplicate(is1, &is1sorted);
387:     ISSort(is1sorted);
388:     ISGetIndices(is1sorted, &i1);
389:   } else {
390:     is1sorted = is1;
391:     PetscObjectReference((PetscObject)is1);
392:     ISGetIndices(is1, &i1);
393:   }
394:   ISSorted(is2, &lsorted);
395:   MPIU_Allreduce(&lsorted, &sorted, 1, MPIU_BOOL, MPI_LAND, comm);
396:   if (!sorted) {
397:     ISDuplicate(is2, &is2sorted);
398:     ISSort(is2sorted);
399:     ISGetIndices(is2sorted, &i2);
400:   } else {
401:     is2sorted = is2;
402:     PetscObjectReference((PetscObject)is2);
403:     ISGetIndices(is2, &i2);
404:   }

406:   PetscMalloc1(n2, &iout);

408:   for (i = 0, nout = 0; i < n2; i++) {
409:     PetscInt key = i2[i];
410:     PetscInt loc;

412:     ISLocate(is1sorted, key, &loc);
413:     if (loc >= 0) {
414:       if (!nout || iout[nout - 1] < key) iout[nout++] = key;
415:     }
416:   }
417:   PetscRealloc(nout * sizeof(PetscInt), &iout);

419:   /* create the new IS containing the sum */
420:   ISCreateGeneral(comm, nout, iout, PETSC_OWN_POINTER, isout);

422:   ISRestoreIndices(is2sorted, &i2);
423:   ISDestroy(&is2sorted);
424:   ISRestoreIndices(is1sorted, &i1);
425:   ISDestroy(&is1sorted);
426:   return 0;
427: }

429: PetscErrorCode ISIntersect_Caching_Internal(IS is1, IS is2, IS *isect)
430: {
431:   *isect = NULL;
432:   if (is2 && is1) {
433:     char          composeStr[33] = {0};
434:     PetscObjectId is2id;

436:     PetscObjectGetId((PetscObject)is2, &is2id);
437:     PetscSNPrintf(composeStr, 32, "ISIntersect_Caching_%" PetscInt64_FMT, is2id);
438:     PetscObjectQuery((PetscObject)is1, composeStr, (PetscObject *)isect);
439:     if (*isect == NULL) {
440:       ISIntersect(is1, is2, isect);
441:       PetscObjectCompose((PetscObject)is1, composeStr, (PetscObject)*isect);
442:     } else {
443:       PetscObjectReference((PetscObject)*isect);
444:     }
445:   }
446:   return 0;
447: }

449: /*@
450:    ISConcatenate - Forms a new `IS` by locally concatenating the indices from an `IS` list without reordering.

452:    Collective.

454:    Input Parameters:
455: +  comm    - communicator of the concatenated IS.
456: .  len     - size of islist array (nonnegative)
457: -  islist  - array of index sets

459:    Output Parameters:
460: .  isout   - The concatenated index set; empty, if len == 0.

462:    Level: intermediate

464:    Notes:
465:     The semantics of calling this on comm imply that the comms of the members of islist also contain this rank.

467: .seealso: [](sec_scatter), `IS`, `ISDifference()`, `ISSum()`, `ISExpand()`, `ISIntersect()`
468: @*/
469: PetscErrorCode ISConcatenate(MPI_Comm comm, PetscInt len, const IS islist[], IS *isout)
470: {
471:   PetscInt        i, n, N;
472:   const PetscInt *iidx;
473:   PetscInt       *idx;

476:   if (PetscDefined(USE_DEBUG)) {
477:     for (i = 0; i < len; ++i)
479:   }
481:   if (!len) {
482:     ISCreateGeneral(comm, 0, NULL, PETSC_OWN_POINTER, isout);
483:     return 0;
484:   }
486:   N = 0;
487:   for (i = 0; i < len; ++i) {
488:     if (islist[i]) {
489:       ISGetLocalSize(islist[i], &n);
490:       N += n;
491:     }
492:   }
493:   PetscMalloc1(N, &idx);
494:   N = 0;
495:   for (i = 0; i < len; ++i) {
496:     if (islist[i]) {
497:       ISGetLocalSize(islist[i], &n);
498:       ISGetIndices(islist[i], &iidx);
499:       PetscArraycpy(idx + N, iidx, n);
500:       ISRestoreIndices(islist[i], &iidx);
501:       N += n;
502:     }
503:   }
504:   ISCreateGeneral(comm, N, idx, PETSC_OWN_POINTER, isout);
505:   return 0;
506: }

508: /*@
509:    ISListToPair  - Convert an `IS` list to a pair of `IS` of equal length defining an equivalent integer multimap.
510:                    Each `IS` on the input list is assigned an integer j so that all of the indices of that `IS` are
511:                    mapped to j.

513:    Collective.

515:    Input Parameters:
516: +  comm    -  `MPI_Comm`
517: .  listlen -  `IS` list length
518: -  islist  -  `IS` list

520:    Output Parameters:
521: +  xis -  domain `IS`
522: -  yis -  range  `I`S

524:   Level: developer

526:   Notes:
527:   The global integers assigned to the `IS` of the local input list might not correspond to the
528:   local numbers of the `IS` on that list, but the two *orderings* are the same: the global
529:   integers assigned to the `IS` on the local list form a strictly increasing sequence.

531:   The `IS` on the input list can belong to subcommunicators of comm, and the subcommunicators
532:   on the input `IS` list are assumed to be in a "deadlock-free" order.

534:   Local lists of `PetscObject`s (or their subcomms) on a comm are "deadlock-free" if subcomm1
535:   preceeds subcomm2 on any local list, then it preceeds subcomm2 on all ranks.
536:   Equivalently, the local numbers of the subcomms on each local list are drawn from some global
537:   numbering. This is ensured, for example, by `ISPairToList()`.

539: .seealso: `IS`, `ISPairToList()`
540: @*/
541: PetscErrorCode ISListToPair(MPI_Comm comm, PetscInt listlen, IS islist[], IS *xis, IS *yis)
542: {
543:   PetscInt        ncolors, *colors, i, leni, len, *xinds, *yinds, k, j;
544:   const PetscInt *indsi;

546:   PetscMalloc1(listlen, &colors);
547:   PetscObjectsListGetGlobalNumbering(comm, listlen, (PetscObject *)islist, &ncolors, colors);
548:   len = 0;
549:   for (i = 0; i < listlen; ++i) {
550:     ISGetLocalSize(islist[i], &leni);
551:     len += leni;
552:   }
553:   PetscMalloc1(len, &xinds);
554:   PetscMalloc1(len, &yinds);
555:   k = 0;
556:   for (i = 0; i < listlen; ++i) {
557:     ISGetLocalSize(islist[i], &leni);
558:     ISGetIndices(islist[i], &indsi);
559:     for (j = 0; j < leni; ++j) {
560:       xinds[k] = indsi[j];
561:       yinds[k] = colors[i];
562:       ++k;
563:     }
564:   }
565:   PetscFree(colors);
566:   ISCreateGeneral(comm, len, xinds, PETSC_OWN_POINTER, xis);
567:   ISCreateGeneral(comm, len, yinds, PETSC_OWN_POINTER, yis);
568:   return 0;
569: }

571: /*@
572:    ISPairToList - Convert an `IS` pair encoding an integer map to a list of `IS`.
573:                   Each `IS` on the output list contains the preimage for each index on the second input `IS`.
574:                   The `IS` on the output list are constructed on the subcommunicators of the input `IS` pair.
575:                   Each subcommunicator corresponds to the preimage of some index j -- this subcomm contains
576:                   exactly the ranks that assign some indices i to j.  This is essentially the inverse of
577:                   `ISListToPair()`.

579:   Collective

581:   Input Parameters:
582: + xis -  domain IS
583: - yis -  range IS

585:   Output Parameters:
586: + listlen -  length of islist
587: - islist  -  list of ISs breaking up indis by color

589:   Level: developer

591:   Notes:
592:   xis and yis must be of the same length and have congruent communicators.

594:   The resulting `IS` have subcommunicators in a "deadlock-free" order (see `ISListToPair()`).

596: .seealso: `IS`, `ISListToPair()`
597:  @*/
598: PetscErrorCode ISPairToList(IS xis, IS yis, PetscInt *listlen, IS **islist)
599: {
600:   IS              indis = xis, coloris = yis;
601:   PetscInt       *inds, *colors, llen, ilen, lstart, lend, lcount, l;
602:   PetscMPIInt     rank, size, llow, lhigh, low, high, color, subsize;
603:   const PetscInt *ccolors, *cinds;
604:   MPI_Comm        comm, subcomm;

611:   PetscObjectGetComm((PetscObject)xis, &comm);
612:   MPI_Comm_rank(comm, &rank);
613:   MPI_Comm_rank(comm, &size);
614:   /* Extract, copy and sort the local indices and colors on the color. */
615:   ISGetLocalSize(coloris, &llen);
616:   ISGetLocalSize(indis, &ilen);
618:   ISGetIndices(coloris, &ccolors);
619:   ISGetIndices(indis, &cinds);
620:   PetscMalloc2(ilen, &inds, llen, &colors);
621:   PetscArraycpy(inds, cinds, ilen);
622:   PetscArraycpy(colors, ccolors, llen);
623:   PetscSortIntWithArray(llen, colors, inds);
624:   /* Determine the global extent of colors. */
625:   llow   = 0;
626:   lhigh  = -1;
627:   lstart = 0;
628:   lcount = 0;
629:   while (lstart < llen) {
630:     lend = lstart + 1;
631:     while (lend < llen && colors[lend] == colors[lstart]) ++lend;
632:     llow  = PetscMin(llow, colors[lstart]);
633:     lhigh = PetscMax(lhigh, colors[lstart]);
634:     ++lcount;
635:   }
636:   MPIU_Allreduce(&llow, &low, 1, MPI_INT, MPI_MIN, comm);
637:   MPIU_Allreduce(&lhigh, &high, 1, MPI_INT, MPI_MAX, comm);
638:   *listlen = 0;
639:   if (low <= high) {
640:     if (lcount > 0) {
641:       *listlen = lcount;
642:       if (!*islist) PetscMalloc1(lcount, islist);
643:     }
644:     /*
645:      Traverse all possible global colors, and participate in the subcommunicators
646:      for the locally-supported colors.
647:      */
648:     lcount = 0;
649:     lstart = 0;
650:     lend   = 0;
651:     for (l = low; l <= high; ++l) {
652:       /*
653:        Find the range of indices with the same color, which is not smaller than l.
654:        Observe that, since colors is sorted, and is a subsequence of [low,high],
655:        as soon as we find a new color, it is >= l.
656:        */
657:       if (lstart < llen) {
658:         /* The start of the next locally-owned color is identified.  Now look for the end. */
659:         if (lstart == lend) {
660:           lend = lstart + 1;
661:           while (lend < llen && colors[lend] == colors[lstart]) ++lend;
662:         }
663:         /* Now check whether the identified color segment matches l. */
665:       }
666:       color = (PetscMPIInt)(colors[lstart] == l);
667:       /* Check whether a proper subcommunicator exists. */
668:       MPIU_Allreduce(&color, &subsize, 1, MPI_INT, MPI_SUM, comm);

670:       if (subsize == 1) subcomm = PETSC_COMM_SELF;
671:       else if (subsize == size) subcomm = comm;
672:       else {
673:         /* a proper communicator is necessary, so we create it. */
674:         MPI_Comm_split(comm, color, rank, &subcomm);
675:       }
676:       if (colors[lstart] == l) {
677:         /* If we have l among the local colors, we create an IS to hold the corresponding indices. */
678:         ISCreateGeneral(subcomm, lend - lstart, inds + lstart, PETSC_COPY_VALUES, *islist + lcount);
679:         /* Position lstart at the beginning of the next local color. */
680:         lstart = lend;
681:         /* Increment the counter of the local colors split off into an IS. */
682:         ++lcount;
683:       }
684:       if (subsize > 0 && subsize < size) {
685:         /*
686:          Irrespective of color, destroy the split off subcomm:
687:          a subcomm used in the IS creation above is duplicated
688:          into a proper PETSc comm.
689:          */
690:         MPI_Comm_free(&subcomm);
691:       }
692:     } /* for (l = low; l < high; ++l) */
693:   }   /* if (low <= high) */
694:   PetscFree2(inds, colors);
695:   return 0;
696: }

698: /*@
699:    ISEmbed - Embed `IS` a into `IS` b by finding the locations in b that have the same indices as in a.
700:              If c is the `IS` of these locations, we have a = b*c, regarded as a composition of the
701:              corresponding `ISLocalToGlobalMap`s.

703:   Not collective.

705:   Input Parameters:
706: + a    -  `IS` to embed
707: . b    -  `IS` to embed into
708: - drop -  flag indicating whether to drop a's indices that are not in b.

710:   Output Parameters:
711: . c    -  local embedding indices

713:   Level: developer

715:   Notes:
716:   If some of a's global indices are not among b's indices the embedding is impossible.  The local indices of a
717:   corresponding to these global indices are either mapped to -1 (if !drop) or are omitted (if drop).  In the former
718:   case the size of c is that same as that of a, in the latter case c's size may be smaller.

720:   The resulting `IS` is sequential, since the index substitution it encodes is purely local.

722: .seealso: `IS`, `ISLocalToGlobalMapping`
723:  @*/
724: PetscErrorCode ISEmbed(IS a, IS b, PetscBool drop, IS *c)
725: {
726:   ISLocalToGlobalMapping     ltog;
727:   ISGlobalToLocalMappingMode gtoltype = IS_GTOLM_DROP;
728:   PetscInt                   alen, clen, *cindices, *cindices2;
729:   const PetscInt            *aindices;

734:   ISLocalToGlobalMappingCreateIS(b, &ltog);
735:   ISGetLocalSize(a, &alen);
736:   ISGetIndices(a, &aindices);
737:   PetscMalloc1(alen, &cindices);
738:   if (!drop) gtoltype = IS_GTOLM_MASK;
739:   ISGlobalToLocalMappingApply(ltog, gtoltype, alen, aindices, &clen, cindices);
740:   ISRestoreIndices(a, &aindices);
741:   ISLocalToGlobalMappingDestroy(&ltog);
742:   if (clen != alen) {
743:     cindices2 = cindices;
744:     PetscMalloc1(clen, &cindices);
745:     PetscArraycpy(cindices, cindices2, clen);
746:     PetscFree(cindices2);
747:   }
748:   ISCreateGeneral(PETSC_COMM_SELF, clen, cindices, PETSC_OWN_POINTER, c);
749:   return 0;
750: }

752: /*@
753:   ISSortPermutation  -  calculate the permutation of the indices into a nondecreasing order.

755:   Not collective.

757:   Input Parameters:
758: + f      -  `IS` to sort
759: - always -  build the permutation even when f's indices are nondecreasing.

761:   Output Parameter:
762: . h    -  permutation or NULL, if f is nondecreasing and always == `PETSC_FALSE`.

764:   Level: advanced

766:   Notes:
767:   Indices in f are unchanged. f[h[i]] is the i-th smallest f index.

769:   If always == `PETSC_FALSE`, an extra check is performed to see whether
770:   the f indices are nondecreasing. h is built on `PETSC_COMM_SELF`, since
771:   the permutation has a local meaning only.

773: .seealso: `IS`, `ISLocalToGlobalMapping`, `ISSort()`
774:  @*/
775: PetscErrorCode ISSortPermutation(IS f, PetscBool always, IS *h)
776: {
777:   const PetscInt *findices;
778:   PetscInt        fsize, *hindices, i;
779:   PetscBool       isincreasing;

783:   ISGetLocalSize(f, &fsize);
784:   ISGetIndices(f, &findices);
785:   *h = NULL;
786:   if (!always) {
787:     isincreasing = PETSC_TRUE;
788:     for (i = 1; i < fsize; ++i) {
789:       if (findices[i] <= findices[i - 1]) {
790:         isincreasing = PETSC_FALSE;
791:         break;
792:       }
793:     }
794:     if (isincreasing) {
795:       ISRestoreIndices(f, &findices);
796:       return 0;
797:     }
798:   }
799:   PetscMalloc1(fsize, &hindices);
800:   for (i = 0; i < fsize; ++i) hindices[i] = i;
801:   PetscSortIntWithPermutation(fsize, findices, hindices);
802:   ISRestoreIndices(f, &findices);
803:   ISCreateGeneral(PETSC_COMM_SELF, fsize, hindices, PETSC_OWN_POINTER, h);
804:   ISSetInfo(*h, IS_PERMUTATION, IS_LOCAL, PETSC_FALSE, PETSC_TRUE);
805:   return 0;
806: }