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