Actual source code: projection.c

      1: #include <petsc/private/vecimpl.h>

      3: /*@
      4:   VecWhichEqual - Creates an index set containing the indices
      5:   where the vectors `Vec1` and `Vec2` have identical elements.

      7:   Collective

      9:   Input Parameters:
     10: + Vec1 - the first vector to compare
     11: - Vec2 - the second two vector to compare

     13:   Output Parameter:
     14: . S - The index set containing the indices i where vec1[i] == vec2[i]

     16:   Level: advanced

     18:   Note:
     19:   The two vectors must have the same parallel layout

     21: .seealso: `Vec`
     22: @*/
     23: PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
     24: {
     25:   PetscInt           i, n_same = 0;
     26:   PetscInt           n, low, high;
     27:   PetscInt          *same = NULL;
     28:   const PetscScalar *v1, *v2;

     30:   PetscFunctionBegin;
     33:   PetscCheckSameComm(Vec1, 1, Vec2, 2);
     34:   VecCheckSameSize(Vec1, 1, Vec2, 2);

     36:   PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
     37:   PetscCall(VecGetLocalSize(Vec1, &n));
     38:   if (n > 0) {
     39:     if (Vec1 == Vec2) {
     40:       PetscCall(VecGetArrayRead(Vec1, &v1));
     41:       v2 = v1;
     42:     } else {
     43:       PetscCall(VecGetArrayRead(Vec1, &v1));
     44:       PetscCall(VecGetArrayRead(Vec2, &v2));
     45:     }

     47:     PetscCall(PetscMalloc1(n, &same));

     49:     for (i = 0; i < n; ++i) {
     50:       if (v1[i] == v2[i]) {
     51:         same[n_same] = low + i;
     52:         ++n_same;
     53:       }
     54:     }

     56:     if (Vec1 == Vec2) {
     57:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
     58:     } else {
     59:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
     60:       PetscCall(VecRestoreArrayRead(Vec2, &v2));
     61:     }
     62:   }
     63:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_same, same, PETSC_OWN_POINTER, S));
     64:   PetscFunctionReturn(PETSC_SUCCESS);
     65: }

     67: /*@
     68:   VecWhichLessThan - Creates an index set containing the indices
     69:   where the vectors `Vec1` < `Vec2`

     71:   Collective

     73:   Input Parameters:
     74: + Vec1 - the first vector to compare
     75: - Vec2 - the second vector to compare

     77:   Output Parameter:
     78: . S - The index set containing the indices i where vec1[i] < vec2[i]

     80:   Level: advanced

     82:   Notes:
     83:   The two vectors must have the same parallel layout

     85:   For complex numbers this only compares the real part

     87: .seealso: `Vec`
     88: @*/
     89: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
     90: {
     91:   PetscInt           i, n_lt = 0;
     92:   PetscInt           n, low, high;
     93:   PetscInt          *lt = NULL;
     94:   const PetscScalar *v1, *v2;

     96:   PetscFunctionBegin;
     99:   PetscCheckSameComm(Vec1, 1, Vec2, 2);
    100:   VecCheckSameSize(Vec1, 1, Vec2, 2);

    102:   PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
    103:   PetscCall(VecGetLocalSize(Vec1, &n));
    104:   if (n > 0) {
    105:     if (Vec1 == Vec2) {
    106:       PetscCall(VecGetArrayRead(Vec1, &v1));
    107:       v2 = v1;
    108:     } else {
    109:       PetscCall(VecGetArrayRead(Vec1, &v1));
    110:       PetscCall(VecGetArrayRead(Vec2, &v2));
    111:     }

    113:     PetscCall(PetscMalloc1(n, &lt));

    115:     for (i = 0; i < n; ++i) {
    116:       if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {
    117:         lt[n_lt] = low + i;
    118:         ++n_lt;
    119:       }
    120:     }

    122:     if (Vec1 == Vec2) {
    123:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
    124:     } else {
    125:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
    126:       PetscCall(VecRestoreArrayRead(Vec2, &v2));
    127:     }
    128:   }
    129:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_lt, lt, PETSC_OWN_POINTER, S));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*@
    134:   VecWhichGreaterThan - Creates an index set containing the indices
    135:   where the vectors `Vec1` > `Vec2`

    137:   Collective

    139:   Input Parameters:
    140: + Vec1 - the first vector to compare
    141: - Vec2 - the second vector to compare

    143:   Output Parameter:
    144: . S - The index set containing the indices i where vec1[i] > vec2[i]

    146:   Level: advanced

    148:   Notes:
    149:   The two vectors must have the same parallel layout

    151:   For complex numbers this only compares the real part

    153: .seealso: `Vec`
    154: @*/
    155: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
    156: {
    157:   PetscInt           i, n_gt = 0;
    158:   PetscInt           n, low, high;
    159:   PetscInt          *gt = NULL;
    160:   const PetscScalar *v1, *v2;

    162:   PetscFunctionBegin;
    165:   PetscCheckSameComm(Vec1, 1, Vec2, 2);
    166:   VecCheckSameSize(Vec1, 1, Vec2, 2);

    168:   PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
    169:   PetscCall(VecGetLocalSize(Vec1, &n));
    170:   if (n > 0) {
    171:     if (Vec1 == Vec2) {
    172:       PetscCall(VecGetArrayRead(Vec1, &v1));
    173:       v2 = v1;
    174:     } else {
    175:       PetscCall(VecGetArrayRead(Vec1, &v1));
    176:       PetscCall(VecGetArrayRead(Vec2, &v2));
    177:     }

    179:     PetscCall(PetscMalloc1(n, &gt));

    181:     for (i = 0; i < n; ++i) {
    182:       if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {
    183:         gt[n_gt] = low + i;
    184:         ++n_gt;
    185:       }
    186:     }

    188:     if (Vec1 == Vec2) {
    189:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
    190:     } else {
    191:       PetscCall(VecRestoreArrayRead(Vec1, &v1));
    192:       PetscCall(VecRestoreArrayRead(Vec2, &v2));
    193:     }
    194:   }
    195:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_gt, gt, PETSC_OWN_POINTER, S));
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: /*@
    200:   VecWhichBetween - Creates an index set containing the indices
    201:   where  `VecLow` < `V` < `VecHigh`

    203:   Collective

    205:   Input Parameters:
    206: + VecLow  - lower bound
    207: . V       - Vector to compare
    208: - VecHigh - higher bound

    210:   Output Parameter:
    211: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]

    213:   Level: advanced

    215:   Notes:
    216:   The vectors must have the same parallel layout

    218:   For complex numbers this only compares the real part

    220: .seealso: `Vec`
    221: @*/
    222: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
    223: {
    224:   PetscInt           i, n_vm = 0;
    225:   PetscInt           n, low, high;
    226:   PetscInt          *vm = NULL;
    227:   const PetscScalar *v1, *v2, *vmiddle;

    229:   PetscFunctionBegin;
    233:   PetscCheckSameComm(V, 2, VecLow, 1);
    234:   PetscCheckSameComm(V, 2, VecHigh, 3);
    235:   VecCheckSameSize(V, 2, VecLow, 1);
    236:   VecCheckSameSize(V, 2, VecHigh, 3);

    238:   PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
    239:   PetscCall(VecGetLocalSize(VecLow, &n));
    240:   if (n > 0) {
    241:     PetscCall(VecGetArrayRead(VecLow, &v1));
    242:     if (VecLow != VecHigh) {
    243:       PetscCall(VecGetArrayRead(VecHigh, &v2));
    244:     } else {
    245:       v2 = v1;
    246:     }
    247:     if (V != VecLow && V != VecHigh) {
    248:       PetscCall(VecGetArrayRead(V, &vmiddle));
    249:     } else if (V == VecLow) {
    250:       vmiddle = v1;
    251:     } else {
    252:       vmiddle = v2;
    253:     }

    255:     PetscCall(PetscMalloc1(n, &vm));

    257:     for (i = 0; i < n; ++i) {
    258:       if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {
    259:         vm[n_vm] = low + i;
    260:         ++n_vm;
    261:       }
    262:     }

    264:     PetscCall(VecRestoreArrayRead(VecLow, &v1));
    265:     if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
    266:     if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
    267:   }
    268:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: /*@
    273:   VecWhichBetweenOrEqual - Creates an index set containing the indices
    274:   where  `VecLow` <= `V` <= `VecHigh`

    276:   Collective

    278:   Input Parameters:
    279: + VecLow  - lower bound
    280: . V       - Vector to compare
    281: - VecHigh - higher bound

    283:   Output Parameter:
    284: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]

    286:   Level: advanced

    288: .seealso: `Vec`
    289: @*/
    290: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
    291: {
    292:   PetscInt           i, n_vm = 0;
    293:   PetscInt           n, low, high;
    294:   PetscInt          *vm = NULL;
    295:   const PetscScalar *v1, *v2, *vmiddle;

    297:   PetscFunctionBegin;
    301:   PetscCheckSameComm(V, 2, VecLow, 1);
    302:   PetscCheckSameComm(V, 2, VecHigh, 3);
    303:   VecCheckSameSize(V, 2, VecLow, 1);
    304:   VecCheckSameSize(V, 2, VecHigh, 3);

    306:   PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
    307:   PetscCall(VecGetLocalSize(VecLow, &n));
    308:   if (n > 0) {
    309:     PetscCall(VecGetArrayRead(VecLow, &v1));
    310:     if (VecLow != VecHigh) {
    311:       PetscCall(VecGetArrayRead(VecHigh, &v2));
    312:     } else {
    313:       v2 = v1;
    314:     }
    315:     if (V != VecLow && V != VecHigh) {
    316:       PetscCall(VecGetArrayRead(V, &vmiddle));
    317:     } else if (V == VecLow) {
    318:       vmiddle = v1;
    319:     } else {
    320:       vmiddle = v2;
    321:     }

    323:     PetscCall(PetscMalloc1(n, &vm));

    325:     for (i = 0; i < n; ++i) {
    326:       if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
    327:         vm[n_vm] = low + i;
    328:         ++n_vm;
    329:       }
    330:     }

    332:     PetscCall(VecRestoreArrayRead(VecLow, &v1));
    333:     if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
    334:     if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
    335:   }
    336:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
    337:   PetscFunctionReturn(PETSC_SUCCESS);
    338: }

    340: /*@
    341:   VecWhichInactive - Creates an `IS` based on a set of vectors

    343:   Collective

    345:   Input Parameters:
    346: + VecLow  - lower bound
    347: . V       - Vector to compare
    348: . D       - Direction to compare
    349: . VecHigh - higher bound
    350: - Strong  - indicator for applying strongly inactive test

    352:   Output Parameter:
    353: . S - The index set containing the indices i where the bound is inactive

    355:   Level: advanced

    357:   Notes:
    358:   Creates an index set containing the indices where one of the following holds\:
    359: .vb
    360:   - VecLow(i)  < V(i) < VecHigh(i)
    361:   - VecLow(i)  = V(i) and D(i) <= 0 (< 0 when Strong is true)
    362:   - VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
    363: .ve

    365: .seealso: `Vec`
    366: @*/
    367: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
    368: {
    369:   PetscInt           i, n_vm = 0;
    370:   PetscInt           n, low, high;
    371:   PetscInt          *vm = NULL;
    372:   const PetscScalar *v1, *v2, *v, *d;

    374:   PetscFunctionBegin;
    375:   if (!VecLow && !VecHigh) {
    376:     PetscCall(VecGetOwnershipRange(V, &low, &high));
    377:     PetscCall(ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S));
    378:     PetscFunctionReturn(PETSC_SUCCESS);
    379:   }
    384:   PetscCheckSameComm(V, 2, VecLow, 1);
    385:   PetscCheckSameComm(V, 2, D, 3);
    386:   PetscCheckSameComm(V, 2, VecHigh, 4);
    387:   VecCheckSameSize(V, 2, VecLow, 1);
    388:   VecCheckSameSize(V, 2, D, 3);
    389:   VecCheckSameSize(V, 2, VecHigh, 4);

    391:   PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
    392:   PetscCall(VecGetLocalSize(VecLow, &n));
    393:   if (n > 0) {
    394:     PetscCall(VecGetArrayRead(VecLow, &v1));
    395:     if (VecLow != VecHigh) {
    396:       PetscCall(VecGetArrayRead(VecHigh, &v2));
    397:     } else {
    398:       v2 = v1;
    399:     }
    400:     if (V != VecLow && V != VecHigh) {
    401:       PetscCall(VecGetArrayRead(V, &v));
    402:     } else if (V == VecLow) {
    403:       v = v1;
    404:     } else {
    405:       v = v2;
    406:     }
    407:     if (D != VecLow && D != VecHigh && D != V) {
    408:       PetscCall(VecGetArrayRead(D, &d));
    409:     } else if (D == VecLow) {
    410:       d = v1;
    411:     } else if (D == VecHigh) {
    412:       d = v2;
    413:     } else {
    414:       d = v;
    415:     }

    417:     PetscCall(PetscMalloc1(n, &vm));

    419:     if (Strong) {
    420:       for (i = 0; i < n; ++i) {
    421:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
    422:           vm[n_vm] = low + i;
    423:           ++n_vm;
    424:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
    425:           vm[n_vm] = low + i;
    426:           ++n_vm;
    427:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
    428:           vm[n_vm] = low + i;
    429:           ++n_vm;
    430:         }
    431:       }
    432:     } else {
    433:       for (i = 0; i < n; ++i) {
    434:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
    435:           vm[n_vm] = low + i;
    436:           ++n_vm;
    437:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
    438:           vm[n_vm] = low + i;
    439:           ++n_vm;
    440:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
    441:           vm[n_vm] = low + i;
    442:           ++n_vm;
    443:         }
    444:       }
    445:     }

    447:     PetscCall(VecRestoreArrayRead(VecLow, &v1));
    448:     if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
    449:     if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &v));
    450:     if (D != VecLow && D != VecHigh && D != V) PetscCall(VecRestoreArrayRead(D, &d));
    451:   }
    452:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
    453:   PetscFunctionReturn(PETSC_SUCCESS);
    454: }

    456: /*@
    457:   VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
    458:   vfull[is[i]] += alpha*vreduced[i]

    460:   Logically Collective

    462:   Input Parameters:
    463: + vfull    - the full-space vector
    464: . is       - the index set for the reduced space
    465: . alpha    - the scalar coefficient
    466: - vreduced - the reduced-space vector

    468:   Output Parameter:
    469: . vfull - the sum of the full-space vector and reduced-space vector

    471:   Level: advanced

    473:   Notes:
    474:   The index set identifies entries in the global vector.
    475:   Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.

    477: .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
    478: @*/
    479: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
    480: {
    481:   PetscInt  nfull, nreduced;
    482:   PetscBool sorted = PETSC_FALSE;

    484:   PetscFunctionBegin;
    487:   PetscCheckSameComm(vfull, 1, is, 2);
    490:   PetscCall(ISGetSize(is, &nfull));
    491:   if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
    492:   PetscCall(VecGetSize(vfull, &nfull));
    493:   PetscCall(VecGetSize(vreduced, &nreduced));
    494:   if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
    495:   if (sorted) PetscCall(VecAXPY(vfull, alpha, vreduced));
    496:   else {
    497:     PetscScalar       *y;
    498:     const PetscScalar *x;
    499:     PetscInt           i, n, m, rstart, rend;
    500:     const PetscInt    *id;

    502:     PetscCall(VecGetArray(vfull, &y));
    503:     PetscCall(VecGetArrayRead(vreduced, &x));
    504:     PetscCall(ISGetIndices(is, &id));
    505:     PetscCall(ISGetLocalSize(is, &n));
    506:     PetscCall(VecGetLocalSize(vreduced, &m));
    507:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
    508:     PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
    509:     y = PetscSafePointerPlusOffset(y, -rstart);
    510:     for (i = 0; i < n; ++i) {
    511:       if (id[i] < 0) continue;
    512:       PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
    513:       y[id[i]] += alpha * x[i];
    514:     }
    515:     y = PetscSafePointerPlusOffset(y, rstart);
    516:     PetscCall(ISRestoreIndices(is, &id));
    517:     PetscCall(VecRestoreArray(vfull, &y));
    518:     PetscCall(VecRestoreArrayRead(vreduced, &x));
    519:   }
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }

    523: /*@
    524:   VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.

    526:   Logically Collective

    528:   Input Parameters:
    529: + vfull    - the full-space vector
    530: . is       - the index set for the reduced space
    531: . mode     - the direction of copying, `SCATTER_FORWARD` or `SCATTER_REVERSE`
    532: - vreduced - the reduced-space vector

    534:   Output Parameter:
    535: . vfull - the sum of the full-space vector and reduced-space vector

    537:   Level: advanced

    539:   Notes:
    540:   The index set identifies entries in the global vector.
    541:   Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
    542: .vb
    543:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
    544:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
    545: .ve

    547: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
    548: @*/
    549: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
    550: {
    551:   PetscInt  nfull, nreduced;
    552:   PetscBool sorted = PETSC_FALSE;

    554:   PetscFunctionBegin;
    557:   PetscCheckSameComm(vfull, 1, is, 2);
    560:   PetscCall(ISGetSize(is, &nfull));
    561:   if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
    562:   PetscCall(VecGetSize(vfull, &nfull));
    563:   PetscCall(VecGetSize(vreduced, &nreduced));
    564:   if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
    565:   if (sorted) {
    566:     if (mode == SCATTER_FORWARD) {
    567:       PetscCall(VecCopy(vreduced, vfull));
    568:     } else {
    569:       PetscCall(VecCopy(vfull, vreduced));
    570:     }
    571:   } else {
    572:     const PetscInt *id;
    573:     PetscInt        i, n, m, rstart, rend;

    575:     PetscCall(ISGetIndices(is, &id));
    576:     PetscCall(ISGetLocalSize(is, &n));
    577:     PetscCall(VecGetLocalSize(vreduced, &m));
    578:     PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
    579:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
    580:     if (mode == SCATTER_FORWARD) {
    581:       PetscScalar       *y;
    582:       const PetscScalar *x;

    584:       PetscCall(VecGetArray(vfull, &y));
    585:       PetscCall(VecGetArrayRead(vreduced, &x));
    586:       y = PetscSafePointerPlusOffset(y, -rstart);
    587:       for (i = 0; i < n; ++i) {
    588:         if (id[i] < 0) continue;
    589:         PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
    590:         y[id[i]] = x[i];
    591:       }
    592:       y = PetscSafePointerPlusOffset(y, rstart);
    593:       PetscCall(VecRestoreArrayRead(vreduced, &x));
    594:       PetscCall(VecRestoreArray(vfull, &y));
    595:     } else if (mode == SCATTER_REVERSE) {
    596:       PetscScalar       *x;
    597:       const PetscScalar *y;

    599:       PetscCall(VecGetArrayRead(vfull, &y));
    600:       PetscCall(VecGetArray(vreduced, &x));
    601:       for (i = 0; i < n; ++i) {
    602:         if (id[i] < 0) continue;
    603:         PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
    604:         x[i] = y[id[i] - rstart];
    605:       }
    606:       PetscCall(VecRestoreArray(vreduced, &x));
    607:       PetscCall(VecRestoreArrayRead(vfull, &y));
    608:     } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
    609:     PetscCall(ISRestoreIndices(is, &id));
    610:   }
    611:   PetscFunctionReturn(PETSC_SUCCESS);
    612: }

    614: /*@
    615:   ISComplementVec - Creates the complement of the index set relative to a layout defined by a `Vec`

    617:   Collective

    619:   Input Parameters:
    620: + S - a PETSc `IS`
    621: - V - the reference vector space

    623:   Output Parameter:
    624: . T - the complement of S

    626:   Level: advanced

    628: .seealso: `IS`, `Vec`, `ISCreateGeneral()`
    629: @*/
    630: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
    631: {
    632:   PetscInt start, end;

    634:   PetscFunctionBegin;
    635:   PetscCall(VecGetOwnershipRange(V, &start, &end));
    636:   PetscCall(ISComplement(S, start, end, T));
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

    640: /*@
    641:   VecISSet - Sets the elements of a vector, specified by an index set, to a constant

    643:   Logically Collective

    645:   Input Parameters:
    646: + V - the vector
    647: . S - index set for the locations in the vector
    648: - c - the constant

    650:   Level: advanced

    652:   Notes:
    653:   The index set identifies entries in the global vector.
    654:   Negative indices are skipped; indices outside the ownership range of V will raise an error.

    656: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISShift()`, `VecSet()`
    657: @*/
    658: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
    659: {
    660:   PetscInt        nloc, low, high, i;
    661:   const PetscInt *s;
    662:   PetscScalar    *v;

    664:   PetscFunctionBegin;
    667:   PetscCheckSameComm(V, 1, S, 2);
    668:   PetscCall(ISGetSize(S, &nloc));
    669:   if (nloc) {
    670:     PetscCall(VecGetOwnershipRange(V, &low, &high));
    671:     PetscCall(ISGetLocalSize(S, &nloc));
    672:     PetscCall(ISGetIndices(S, &s));
    673:     PetscCall(VecGetArray(V, &v));
    674:     for (i = 0; i < nloc; ++i) {
    675:       if (s[i] < 0) continue;
    676:       PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
    677:       v[s[i] - low] = c;
    678:     }
    679:     PetscCall(ISRestoreIndices(S, &s));
    680:     PetscCall(VecRestoreArray(V, &v));
    681:   }
    682:   PetscFunctionReturn(PETSC_SUCCESS);
    683: }

    685: /*@
    686:   VecISShift - Shifts the elements of a vector, specified by an index set, by a constant

    688:   Logically Collective

    690:   Input Parameters:
    691: + V - the vector
    692: . S - index set for the locations in the vector
    693: - c - the constant

    695:   Level: advanced

    697:   Notes:
    698:   The index set identifies entries in the global vector.
    699:   Negative indices are skipped; indices outside the ownership range of V will raise an error.

    701: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISSet()`, `VecShift()`
    702: @*/
    703: PetscErrorCode VecISShift(Vec V, IS S, PetscScalar c)
    704: {
    705:   PetscInt        nloc, low, high, i;
    706:   const PetscInt *s;
    707:   PetscScalar    *v;

    709:   PetscFunctionBegin;
    712:   PetscCheckSameComm(V, 1, S, 2);
    713:   PetscCall(ISGetSize(S, &nloc));
    714:   if (nloc) {
    715:     PetscCall(VecGetOwnershipRange(V, &low, &high));
    716:     PetscCall(ISGetLocalSize(S, &nloc));
    717:     PetscCall(ISGetIndices(S, &s));
    718:     PetscCall(VecGetArray(V, &v));
    719:     for (i = 0; i < nloc; ++i) {
    720:       if (s[i] < 0) continue;
    721:       PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
    722:       v[s[i] - low] += c;
    723:     }
    724:     PetscCall(ISRestoreIndices(S, &s));
    725:     PetscCall(VecRestoreArray(V, &v));
    726:   }
    727:   PetscFunctionReturn(PETSC_SUCCESS);
    728: }

    730: /*@
    731:   VecBoundGradientProjection - Projects vector according to this definition.
    732:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
    733:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
    734:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

    736:   Input Parameters:
    737: + G  - current gradient vector
    738: . X  - current solution vector with XL[i] <= X[i] <= XU[i]
    739: . XL - lower bounds
    740: - XU - upper bounds

    742:   Output Parameter:
    743: . GP - gradient projection vector

    745:   Level: advanced

    747:   Note:
    748:   `GP` may be the same vector as `G`

    750:   For complex numbers only the real part is used in the bounds.

    752: .seealso: `Vec`
    753: @*/
    754: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
    755: {
    756:   PetscInt           n, i;
    757:   const PetscScalar *xptr, *xlptr, *xuptr;
    758:   PetscScalar       *gptr, *gpptr;
    759:   PetscScalar        xval, gpval;

    761:   /* Project variables at the lower and upper bound */
    762:   PetscFunctionBegin;
    768:   if (!XL && !XU) {
    769:     PetscCall(VecCopy(G, GP));
    770:     PetscFunctionReturn(PETSC_SUCCESS);
    771:   }

    773:   PetscCall(VecGetLocalSize(X, &n));

    775:   PetscCall(VecGetArrayRead(X, &xptr));
    776:   PetscCall(VecGetArrayRead(XL, &xlptr));
    777:   PetscCall(VecGetArrayRead(XU, &xuptr));
    778:   PetscCall(VecGetArrayPair(G, GP, &gptr, &gpptr));

    780:   for (i = 0; i < n; ++i) {
    781:     gpval = gptr[i];
    782:     xval  = xptr[i];
    783:     if (PetscRealPart(gpval) > 0.0 && PetscRealPart(xval) <= PetscRealPart(xlptr[i])) {
    784:       gpval = 0.0;
    785:     } else if (PetscRealPart(gpval) < 0.0 && PetscRealPart(xval) >= PetscRealPart(xuptr[i])) {
    786:       gpval = 0.0;
    787:     }
    788:     gpptr[i] = gpval;
    789:   }

    791:   PetscCall(VecRestoreArrayRead(X, &xptr));
    792:   PetscCall(VecRestoreArrayRead(XL, &xlptr));
    793:   PetscCall(VecRestoreArrayRead(XU, &xuptr));
    794:   PetscCall(VecRestoreArrayPair(G, GP, &gptr, &gpptr));
    795:   PetscFunctionReturn(PETSC_SUCCESS);
    796: }

    798: /*@
    799:   VecStepMaxBounded - See below

    801:   Collective

    803:   Input Parameters:
    804: + X  - vector with no negative entries
    805: . XL - lower bounds
    806: . XU - upper bounds
    807: - DX - step direction, can have negative, positive or zero entries

    809:   Output Parameter:
    810: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + stepmax*DX[i]

    812:   Level: intermediate

    814: .seealso: `Vec`
    815: @*/
    816: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
    817: {
    818:   PetscInt           i, nn;
    819:   const PetscScalar *xx, *dx, *xl, *xu;
    820:   PetscReal          localmax = 0;

    822:   PetscFunctionBegin;

    828:   PetscCall(VecGetArrayRead(X, &xx));
    829:   PetscCall(VecGetArrayRead(XL, &xl));
    830:   PetscCall(VecGetArrayRead(XU, &xu));
    831:   PetscCall(VecGetArrayRead(DX, &dx));
    832:   PetscCall(VecGetLocalSize(X, &nn));
    833:   for (i = 0; i < nn; i++) {
    834:     if (PetscRealPart(dx[i]) > 0) {
    835:       localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
    836:     } else if (PetscRealPart(dx[i]) < 0) {
    837:       localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
    838:     }
    839:   }
    840:   PetscCall(VecRestoreArrayRead(X, &xx));
    841:   PetscCall(VecRestoreArrayRead(XL, &xl));
    842:   PetscCall(VecRestoreArrayRead(XU, &xu));
    843:   PetscCall(VecRestoreArrayRead(DX, &dx));
    844:   PetscCallMPI(MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X)));
    845:   PetscFunctionReturn(PETSC_SUCCESS);
    846: }

    848: /*@
    849:   VecStepBoundInfo - See below

    851:   Collective

    853:   Input Parameters:
    854: + X  - vector with no negative entries
    855: . XL - lower bounds
    856: . XU - upper bounds
    857: - DX - step direction, can have negative, positive or zero entries

    859:   Output Parameters:
    860: + boundmin - (may be `NULL` this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
    861: . wolfemin - (may be `NULL` this it is not computed)
    862: - boundmax - (may be `NULL` this it is not computed) minimum value so that X[i] + boundmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + boundmax*DX[i]

    864:   Level: advanced

    866:   Note:
    867:   For complex numbers only compares the real part

    869: .seealso: `Vec`
    870: @*/
    871: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
    872: {
    873:   PetscInt           n, i;
    874:   const PetscScalar *x, *xl, *xu, *dx;
    875:   PetscReal          t;
    876:   PetscReal          localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
    877:   MPI_Comm           comm;

    879:   PetscFunctionBegin;

    885:   PetscCall(VecGetArrayRead(X, &x));
    886:   PetscCall(VecGetArrayRead(XL, &xl));
    887:   PetscCall(VecGetArrayRead(XU, &xu));
    888:   PetscCall(VecGetArrayRead(DX, &dx));
    889:   PetscCall(VecGetLocalSize(X, &n));
    890:   for (i = 0; i < n; ++i) {
    891:     if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
    892:       t        = PetscRealPart((xu[i] - x[i]) / dx[i]);
    893:       localmin = PetscMin(t, localmin);
    894:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
    895:       localmax = PetscMax(t, localmax);
    896:     } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
    897:       t        = PetscRealPart((xl[i] - x[i]) / dx[i]);
    898:       localmin = PetscMin(t, localmin);
    899:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
    900:       localmax = PetscMax(t, localmax);
    901:     }
    902:   }

    904:   PetscCall(VecRestoreArrayRead(X, &x));
    905:   PetscCall(VecRestoreArrayRead(XL, &xl));
    906:   PetscCall(VecRestoreArrayRead(XU, &xu));
    907:   PetscCall(VecRestoreArrayRead(DX, &dx));
    908:   PetscCall(PetscObjectGetComm((PetscObject)X, &comm));

    910:   if (boundmin) {
    911:     PetscCallMPI(MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm));
    912:     PetscCall(PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin));
    913:   }
    914:   if (wolfemin) {
    915:     PetscCallMPI(MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm));
    916:     PetscCall(PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin));
    917:   }
    918:   if (boundmax) {
    919:     PetscCallMPI(MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm));
    920:     if (*boundmax < 0) *boundmax = PETSC_INFINITY;
    921:     PetscCall(PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax));
    922:   }
    923:   PetscFunctionReturn(PETSC_SUCCESS);
    924: }

    926: /*@
    927:   VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i

    929:   Collective

    931:   Input Parameters:
    932: + X  - vector with no negative entries
    933: - DX - a step direction, can have negative, positive or zero entries

    935:   Output Parameter:
    936: . step - largest value such that x[i] + step*DX[i] >= 0 for all i

    938:   Level: advanced

    940:   Note:
    941:   For complex numbers only compares the real part

    943: .seealso: `Vec`
    944: @*/
    945: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
    946: {
    947:   PetscInt           i, nn;
    948:   PetscReal          stepmax = PETSC_INFINITY;
    949:   const PetscScalar *xx, *dx;

    951:   PetscFunctionBegin;

    955:   PetscCall(VecGetLocalSize(X, &nn));
    956:   PetscCall(VecGetArrayRead(X, &xx));
    957:   PetscCall(VecGetArrayRead(DX, &dx));
    958:   for (i = 0; i < nn; ++i) {
    959:     PetscCheck(PetscRealPart(xx[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector must be positive");
    960:     if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
    961:   }
    962:   PetscCall(VecRestoreArrayRead(X, &xx));
    963:   PetscCall(VecRestoreArrayRead(DX, &dx));
    964:   PetscCallMPI(MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X)));
    965:   PetscFunctionReturn(PETSC_SUCCESS);
    966: }

    968: /*@
    969:   VecPow - Replaces each component of a vector by $ x_i^p $

    971:   Logically Collective

    973:   Input Parameters:
    974: + v - the vector
    975: - p - the exponent to use on each element

    977:   Level: intermediate

    979:   Note:
    980:   This handles negative values, Inf, and Nan in the expected IEEE floating pointing manner. For example, the square root of a negative real number is Nan
    981:   and 1/0.0 is Inf.

    983: .seealso: `Vec`
    984: @*/
    985: PetscErrorCode VecPow(Vec v, PetscScalar p)
    986: {
    987:   PetscInt     n, i;
    988:   PetscScalar *v1;

    990:   PetscFunctionBegin;

    993:   PetscCall(VecGetArray(v, &v1));
    994:   PetscCall(VecGetLocalSize(v, &n));

    996:   if (1.0 == p) {
    997:   } else if (-1.0 == p) {
    998:     for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
    999:   } else if (0.0 == p) {
    1000:     for (i = 0; i < n; ++i) {
    1001:       /*  Not-a-number left alone
    1002:           Infinity set to one  */
    1003:       if (!PetscIsNanScalar(v1[i])) v1[i] = 1.0;
    1004:     }
    1005:   } else if (0.5 == p) {
    1006:     for (i = 0; i < n; ++i) v1[i] = PetscSqrtScalar(v1[i]);
    1007:   } else if (-0.5 == p) {
    1008:     for (i = 0; i < n; ++i) v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
    1009:   } else if (2.0 == p) {
    1010:     for (i = 0; i < n; ++i) v1[i] *= v1[i];
    1011:   } else if (-2.0 == p) {
    1012:     for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
    1013:   } else {
    1014:     for (i = 0; i < n; ++i) v1[i] = PetscPowScalar(v1[i], p);
    1015:   }
    1016:   PetscCall(VecRestoreArray(v, &v1));
    1017:   PetscFunctionReturn(PETSC_SUCCESS);
    1018: }

    1020: /*@
    1021:   VecMedian - Computes the componentwise median of three vectors
    1022:   and stores the result in this vector.  Used primarily for projecting
    1023:   a vector within upper and lower bounds.

    1025:   Logically Collective

    1027:   Input Parameters:
    1028: + Vec1 - The first vector
    1029: . Vec2 - The second vector
    1030: - Vec3 - The third vector

    1032:   Output Parameter:
    1033: . VMedian - The median vector (this can be any one of the input vectors)

    1035:   Level: advanced

    1037: .seealso: `Vec`
    1038: @*/
    1039: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
    1040: {
    1041:   PetscInt           i, n, low1, high1;
    1042:   const PetscScalar *v1, *v2, *v3;
    1043:   PetscScalar       *vmed;

    1045:   PetscFunctionBegin;

    1051:   if (!Vec1 && !Vec3) {
    1052:     PetscCall(VecCopy(Vec2, VMedian));
    1053:     PetscFunctionReturn(PETSC_SUCCESS);
    1054:   }
    1055:   if (Vec1 == Vec2 || Vec1 == Vec3) {
    1056:     PetscCall(VecCopy(Vec1, VMedian));
    1057:     PetscFunctionReturn(PETSC_SUCCESS);
    1058:   }
    1059:   if (Vec2 == Vec3) {
    1060:     PetscCall(VecCopy(Vec2, VMedian));
    1061:     PetscFunctionReturn(PETSC_SUCCESS);
    1062:   }

    1064:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
    1069:   PetscCheckSameType(Vec1, 1, Vec2, 2);
    1070:   PetscCheckSameType(Vec1, 1, Vec3, 3);
    1071:   PetscCheckSameType(Vec1, 1, VMedian, 4);
    1072:   PetscCheckSameComm(Vec1, 1, Vec2, 2);
    1073:   PetscCheckSameComm(Vec1, 1, Vec3, 3);
    1074:   PetscCheckSameComm(Vec1, 1, VMedian, 4);
    1075:   VecCheckSameSize(Vec1, 1, Vec2, 2);
    1076:   VecCheckSameSize(Vec1, 1, Vec3, 3);
    1077:   VecCheckSameSize(Vec1, 1, VMedian, 4);

    1079:   PetscCall(VecGetOwnershipRange(Vec1, &low1, &high1));
    1080:   PetscCall(VecGetLocalSize(Vec1, &n));
    1081:   if (n > 0) {
    1082:     PetscCall(VecGetArray(VMedian, &vmed));
    1083:     if (Vec1 != VMedian) {
    1084:       PetscCall(VecGetArrayRead(Vec1, &v1));
    1085:     } else {
    1086:       v1 = vmed;
    1087:     }
    1088:     if (Vec2 != VMedian) {
    1089:       PetscCall(VecGetArrayRead(Vec2, &v2));
    1090:     } else {
    1091:       v2 = vmed;
    1092:     }
    1093:     if (Vec3 != VMedian) {
    1094:       PetscCall(VecGetArrayRead(Vec3, &v3));
    1095:     } else {
    1096:       v3 = vmed;
    1097:     }

    1099:     for (i = 0; i < n; ++i) vmed[i] = PetscMax(PetscMax(PetscMin(PetscRealPart(v1[i]), PetscRealPart(v2[i])), PetscMin(PetscRealPart(v1[i]), PetscRealPart(v3[i]))), PetscMin(PetscRealPart(v2[i]), PetscRealPart(v3[i])));

    1101:     PetscCall(VecRestoreArray(VMedian, &vmed));
    1102:     if (VMedian != Vec1) PetscCall(VecRestoreArrayRead(Vec1, &v1));
    1103:     if (VMedian != Vec2) PetscCall(VecRestoreArrayRead(Vec2, &v2));
    1104:     if (VMedian != Vec3) PetscCall(VecRestoreArrayRead(Vec3, &v3));
    1105:   }
    1106:   PetscFunctionReturn(PETSC_SUCCESS);
    1107: }