Actual source code: tao_util.c

      1: #include <petsc/private/petscimpl.h>
      2: #include <petsctao.h>
      3: #include <petscsys.h>

      5: static inline PetscReal Fischer(PetscReal a, PetscReal b)
      6: {
      7:   /* Method suggested by Bob Vanderbei */
      8:   if (a + b <= 0) return PetscSqrtReal(a * a + b * b) - (a + b);
      9:   return -2.0 * a * b / (PetscSqrtReal(a * a + b * b) + (a + b));
     10: }

     12: /*@
     13:   VecFischer - Evaluates the Fischer-Burmeister function for complementarity
     14:   problems.

     16:   Logically Collective

     18:   Input Parameters:
     19: + X - current point
     20: . F - function evaluated at x
     21: . L - lower bounds
     22: - U - upper bounds

     24:   Output Parameter:
     25: . FB - The Fischer-Burmeister function vector

     27:   Level: developer

     29:   Notes:
     30:   The Fischer-Burmeister function is defined as
     31: $        phi(a,b) := sqrt(a*a + b*b) - a - b
     32:   and is used reformulate a complementarity problem as a semismooth
     33:   system of equations.

     35:   The result of this function is done by cases\:
     36: +  l[i] == -infinity, u[i] == infinity  -- fb[i] = -f[i]
     37: .  l[i] == -infinity, u[i] finite       -- fb[i] = phi(u[i]-x[i], -f[i])
     38: .  l[i] finite,       u[i] == infinity  -- fb[i] = phi(x[i]-l[i],  f[i])
     39: .  l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
     40: -  otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]

     42: .seealso: `Vec`, `VecSFischer()`, `MatDFischer()`, `MatDSFischer()`
     43: @*/
     44: PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB)
     45: {
     46:   const PetscScalar *x, *f, *l, *u;
     47:   PetscScalar       *fb;
     48:   PetscReal          xval, fval, lval, uval;
     49:   PetscInt           low[5], high[5], n, i;

     51:   PetscFunctionBegin;

     58:   if (!L && !U) {
     59:     PetscCall(VecAXPBY(FB, -1.0, 0.0, F));
     60:     PetscFunctionReturn(PETSC_SUCCESS);
     61:   }

     63:   PetscCall(VecGetOwnershipRange(X, low, high));
     64:   PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
     65:   PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
     66:   PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
     67:   PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));

     69:   for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors");

     71:   PetscCall(VecGetArrayRead(X, &x));
     72:   PetscCall(VecGetArrayRead(F, &f));
     73:   PetscCall(VecGetArrayRead(L, &l));
     74:   PetscCall(VecGetArrayRead(U, &u));
     75:   PetscCall(VecGetArray(FB, &fb));

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

     79:   for (i = 0; i < n; ++i) {
     80:     xval = PetscRealPart(x[i]);
     81:     fval = PetscRealPart(f[i]);
     82:     lval = PetscRealPart(l[i]);
     83:     uval = PetscRealPart(u[i]);

     85:     if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) {
     86:       fb[i] = -fval;
     87:     } else if (lval <= -PETSC_INFINITY) {
     88:       fb[i] = -Fischer(uval - xval, -fval);
     89:     } else if (uval >= PETSC_INFINITY) {
     90:       fb[i] = Fischer(xval - lval, fval);
     91:     } else if (lval == uval) {
     92:       fb[i] = lval - xval;
     93:     } else {
     94:       fval  = Fischer(uval - xval, -fval);
     95:       fb[i] = Fischer(xval - lval, fval);
     96:     }
     97:   }

     99:   PetscCall(VecRestoreArrayRead(X, &x));
    100:   PetscCall(VecRestoreArrayRead(F, &f));
    101:   PetscCall(VecRestoreArrayRead(L, &l));
    102:   PetscCall(VecRestoreArrayRead(U, &u));
    103:   PetscCall(VecRestoreArray(FB, &fb));
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: static inline PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c)
    108: {
    109:   /* Method suggested by Bob Vanderbei */
    110:   if (a + b <= 0) return PetscSqrtReal(a * a + b * b + 2.0 * c * c) - (a + b);
    111:   return 2.0 * (c * c - a * b) / (PetscSqrtReal(a * a + b * b + 2.0 * c * c) + (a + b));
    112: }

    114: /*@
    115:   VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
    116:   complementarity problems.

    118:   Logically Collective

    120:   Input Parameters:
    121: + X  - current point
    122: . F  - function evaluated at x
    123: . L  - lower bounds
    124: . U  - upper bounds
    125: - mu - smoothing parameter

    127:   Output Parameter:
    128: . FB - The Smoothed Fischer-Burmeister function vector

    130:   Notes:
    131:   The Smoothed Fischer-Burmeister function is defined as
    132: $        phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b
    133:   and is used reformulate a complementarity problem as a semismooth
    134:   system of equations.

    136:   The result of this function is done by cases\:
    137: +  l[i] == -infinity, u[i] == infinity  -- fb[i] = -f[i] - 2*mu*x[i]
    138: .  l[i] == -infinity, u[i] finite       -- fb[i] = phi(u[i]-x[i], -f[i], mu)
    139: .  l[i] finite,       u[i] == infinity  -- fb[i] = phi(x[i]-l[i],  f[i], mu)
    140: .  l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
    141: -  otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]

    143:   Level: developer

    145: .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()`
    146: @*/
    147: PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB)
    148: {
    149:   const PetscScalar *x, *f, *l, *u;
    150:   PetscScalar       *fb;
    151:   PetscReal          xval, fval, lval, uval;
    152:   PetscInt           low[5], high[5], n, i;

    154:   PetscFunctionBegin;

    161:   PetscCall(VecGetOwnershipRange(X, low, high));
    162:   PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
    163:   PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
    164:   PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
    165:   PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));

    167:   for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors");

    169:   PetscCall(VecGetArrayRead(X, &x));
    170:   PetscCall(VecGetArrayRead(F, &f));
    171:   PetscCall(VecGetArrayRead(L, &l));
    172:   PetscCall(VecGetArrayRead(U, &u));
    173:   PetscCall(VecGetArray(FB, &fb));

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

    177:   for (i = 0; i < n; ++i) {
    178:     xval = PetscRealPart(*x++);
    179:     fval = PetscRealPart(*f++);
    180:     lval = PetscRealPart(*l++);
    181:     uval = PetscRealPart(*u++);

    183:     if ((lval <= -PETSC_INFINITY) && (uval >= PETSC_INFINITY)) {
    184:       (*fb++) = -fval - mu * xval;
    185:     } else if (lval <= -PETSC_INFINITY) {
    186:       (*fb++) = -SFischer(uval - xval, -fval, mu);
    187:     } else if (uval >= PETSC_INFINITY) {
    188:       (*fb++) = SFischer(xval - lval, fval, mu);
    189:     } else if (lval == uval) {
    190:       (*fb++) = lval - xval;
    191:     } else {
    192:       fval    = SFischer(uval - xval, -fval, mu);
    193:       (*fb++) = SFischer(xval - lval, fval, mu);
    194:     }
    195:   }
    196:   x -= n;
    197:   f -= n;
    198:   l -= n;
    199:   u -= n;
    200:   fb -= n;

    202:   PetscCall(VecRestoreArrayRead(X, &x));
    203:   PetscCall(VecRestoreArrayRead(F, &f));
    204:   PetscCall(VecRestoreArrayRead(L, &l));
    205:   PetscCall(VecRestoreArrayRead(U, &u));
    206:   PetscCall(VecRestoreArray(FB, &fb));
    207:   PetscFunctionReturn(PETSC_SUCCESS);
    208: }

    210: static inline PetscReal fischnorm(PetscReal a, PetscReal b)
    211: {
    212:   return PetscSqrtReal(a * a + b * b);
    213: }

    215: static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
    216: {
    217:   return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
    218: }

    220: /*@
    221:   MatDFischer - Calculates an element of the B-subdifferential of the
    222:   Fischer-Burmeister function for complementarity problems.

    224:   Collective

    226:   Input Parameters:
    227: + jac - the jacobian of `f` at `X`
    228: . X   - current point
    229: . Con - constraints function evaluated at `X`
    230: . XL  - lower bounds
    231: . XU  - upper bounds
    232: . T1  - work vector
    233: - T2  - work vector

    235:   Output Parameters:
    236: + Da - diagonal perturbation component of the result
    237: - Db - row scaling component of the result

    239:   Level: developer

    241: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()`
    242: @*/
    243: PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db)
    244: {
    245:   PetscInt           i, nn;
    246:   const PetscScalar *x, *f, *l, *u, *t2;
    247:   PetscScalar       *da, *db, *t1;
    248:   PetscReal          ai, bi, ci, di, ei;

    250:   PetscFunctionBegin;
    251:   PetscCall(VecGetLocalSize(X, &nn));
    252:   PetscCall(VecGetArrayRead(X, &x));
    253:   PetscCall(VecGetArrayRead(Con, &f));
    254:   PetscCall(VecGetArrayRead(XL, &l));
    255:   PetscCall(VecGetArrayRead(XU, &u));
    256:   PetscCall(VecGetArray(Da, &da));
    257:   PetscCall(VecGetArray(Db, &db));
    258:   PetscCall(VecGetArray(T1, &t1));
    259:   PetscCall(VecGetArrayRead(T2, &t2));

    261:   for (i = 0; i < nn; i++) {
    262:     da[i] = 0.0;
    263:     db[i] = 0.0;
    264:     t1[i] = 0.0;

    266:     if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) {
    267:       if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) {
    268:         t1[i] = 1.0;
    269:         da[i] = 1.0;
    270:       }

    272:       if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) {
    273:         t1[i] = 1.0;
    274:         db[i] = 1.0;
    275:       }
    276:     }
    277:   }

    279:   PetscCall(VecRestoreArray(T1, &t1));
    280:   PetscCall(VecRestoreArrayRead(T2, &t2));
    281:   PetscCall(MatMult(jac, T1, T2));
    282:   PetscCall(VecGetArrayRead(T2, &t2));

    284:   for (i = 0; i < nn; i++) {
    285:     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
    286:       da[i] = 0.0;
    287:       db[i] = -1.0;
    288:     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
    289:       if (PetscRealPart(db[i]) >= 1) {
    290:         ai = fischnorm(1.0, PetscRealPart(t2[i]));

    292:         da[i] = -1.0 / ai - 1.0;
    293:         db[i] = -t2[i] / ai - 1.0;
    294:       } else {
    295:         bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
    296:         ai = fischnorm(bi, PetscRealPart(f[i]));
    297:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    299:         da[i] = bi / ai - 1.0;
    300:         db[i] = -f[i] / ai - 1.0;
    301:       }
    302:     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
    303:       if (PetscRealPart(da[i]) >= 1) {
    304:         ai = fischnorm(1.0, PetscRealPart(t2[i]));

    306:         da[i] = 1.0 / ai - 1.0;
    307:         db[i] = t2[i] / ai - 1.0;
    308:       } else {
    309:         bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
    310:         ai = fischnorm(bi, PetscRealPart(f[i]));
    311:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    313:         da[i] = bi / ai - 1.0;
    314:         db[i] = f[i] / ai - 1.0;
    315:       }
    316:     } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
    317:       da[i] = -1.0;
    318:       db[i] = 0.0;
    319:     } else {
    320:       if (PetscRealPart(db[i]) >= 1) {
    321:         ai = fischnorm(1.0, PetscRealPart(t2[i]));

    323:         ci = 1.0 / ai + 1.0;
    324:         di = PetscRealPart(t2[i]) / ai + 1.0;
    325:       } else {
    326:         bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
    327:         ai = fischnorm(bi, PetscRealPart(f[i]));
    328:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    330:         ci = bi / ai + 1.0;
    331:         di = PetscRealPart(f[i]) / ai + 1.0;
    332:       }

    334:       if (PetscRealPart(da[i]) >= 1) {
    335:         bi = ci + di * PetscRealPart(t2[i]);
    336:         ai = fischnorm(1.0, bi);

    338:         bi = bi / ai - 1.0;
    339:         ai = 1.0 / ai - 1.0;
    340:       } else {
    341:         ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]));
    342:         ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei);
    343:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    345:         bi = ei / ai - 1.0;
    346:         ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
    347:       }

    349:       da[i] = ai + bi * ci;
    350:       db[i] = bi * di;
    351:     }
    352:   }

    354:   PetscCall(VecRestoreArray(Da, &da));
    355:   PetscCall(VecRestoreArray(Db, &db));
    356:   PetscCall(VecRestoreArrayRead(X, &x));
    357:   PetscCall(VecRestoreArrayRead(Con, &f));
    358:   PetscCall(VecRestoreArrayRead(XL, &l));
    359:   PetscCall(VecRestoreArrayRead(XU, &u));
    360:   PetscCall(VecRestoreArrayRead(T2, &t2));
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }

    364: /*@
    365:   MatDSFischer - Calculates an element of the B-subdifferential of the
    366:   smoothed Fischer-Burmeister function for complementarity problems.

    368:   Collective

    370:   Input Parameters:
    371: + jac - the jacobian of f at X
    372: . X   - current point
    373: . Con - constraint function evaluated at X
    374: . XL  - lower bounds
    375: . XU  - upper bounds
    376: . mu  - smoothing parameter
    377: . T1  - work vector
    378: - T2  - work vector

    380:   Output Parameters:
    381: + Da - diagonal perturbation component of the result
    382: . Db - row scaling component of the result
    383: - Dm - derivative with respect to scaling parameter

    385:   Level: developer

    387: .seealso: `Mat`, `VecFischer()`, `VecDFischer()`, `MatDFischer()`
    388: @*/
    389: PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, PetscReal mu, Vec T1, Vec T2, Vec Da, Vec Db, Vec Dm)
    390: {
    391:   PetscInt           i, nn;
    392:   const PetscScalar *x, *f, *l, *u;
    393:   PetscScalar       *da, *db, *dm;
    394:   PetscReal          ai, bi, ci, di, ei, fi;

    396:   PetscFunctionBegin;
    397:   if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) {
    398:     PetscCall(VecZeroEntries(Dm));
    399:     PetscCall(MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db));
    400:   } else {
    401:     PetscCall(VecGetLocalSize(X, &nn));
    402:     PetscCall(VecGetArrayRead(X, &x));
    403:     PetscCall(VecGetArrayRead(Con, &f));
    404:     PetscCall(VecGetArrayRead(XL, &l));
    405:     PetscCall(VecGetArrayRead(XU, &u));
    406:     PetscCall(VecGetArray(Da, &da));
    407:     PetscCall(VecGetArray(Db, &db));
    408:     PetscCall(VecGetArray(Dm, &dm));

    410:     for (i = 0; i < nn; ++i) {
    411:       if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
    412:         da[i] = -mu;
    413:         db[i] = -1.0;
    414:         dm[i] = -x[i];
    415:       } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
    416:         bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
    417:         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
    418:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    420:         da[i] = bi / ai - 1.0;
    421:         db[i] = -PetscRealPart(f[i]) / ai - 1.0;
    422:         dm[i] = 2.0 * mu / ai;
    423:       } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
    424:         bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
    425:         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
    426:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    428:         da[i] = bi / ai - 1.0;
    429:         db[i] = PetscRealPart(f[i]) / ai - 1.0;
    430:         dm[i] = 2.0 * mu / ai;
    431:       } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
    432:         da[i] = -1.0;
    433:         db[i] = 0.0;
    434:         dm[i] = 0.0;
    435:       } else {
    436:         bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
    437:         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
    438:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    440:         ci = bi / ai + 1.0;
    441:         di = PetscRealPart(f[i]) / ai + 1.0;
    442:         fi = 2.0 * mu / ai;

    444:         ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu);
    445:         ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu);
    446:         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);

    448:         bi = ei / ai - 1.0;
    449:         ei = 2.0 * mu / ei;
    450:         ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;

    452:         da[i] = ai + bi * ci;
    453:         db[i] = bi * di;
    454:         dm[i] = ei + bi * fi;
    455:       }
    456:     }

    458:     PetscCall(VecRestoreArrayRead(X, &x));
    459:     PetscCall(VecRestoreArrayRead(Con, &f));
    460:     PetscCall(VecRestoreArrayRead(XL, &l));
    461:     PetscCall(VecRestoreArrayRead(XU, &u));
    462:     PetscCall(VecRestoreArray(Da, &da));
    463:     PetscCall(VecRestoreArray(Db, &db));
    464:     PetscCall(VecRestoreArray(Dm, &dm));
    465:   }
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: static inline PetscReal ST_InternalPN(PetscScalar in, PetscReal lb, PetscReal ub)
    470: {
    471:   return PetscMax(0, (PetscReal)PetscRealPart(in) - ub) - PetscMax(0, -(PetscReal)PetscRealPart(in) - PetscAbsReal(lb));
    472: }

    474: static inline PetscReal ST_InternalNN(PetscScalar in, PetscReal lb, PetscReal ub)
    475: {
    476:   return PetscMax(0, (PetscReal)PetscRealPart(in) + PetscAbsReal(ub)) - PetscMax(0, -(PetscReal)PetscRealPart(in) - PetscAbsReal(lb));
    477: }

    479: static inline PetscReal ST_InternalPP(PetscScalar in, PetscReal lb, PetscReal ub)
    480: {
    481:   return PetscMax(0, (PetscReal)PetscRealPart(in) - ub) + PetscMin(0, (PetscReal)PetscRealPart(in) - lb);
    482: }

    484: /*@
    485:   TaoSoftThreshold - Calculates soft thresholding routine with input vector
    486:   and given lower and upper bound and returns it to output vector.

    488:   Input Parameters:
    489: + in - input vector to be thresholded
    490: . lb - lower bound
    491: - ub - upper bound

    493:   Output Parameter:
    494: . out - Soft thresholded output vector

    496:   Notes:
    497:   Soft thresholding is defined as
    498:   \[ S(input,lb,ub) =
    499:   \begin{cases}
    500:   input - ub  \text{input > ub} \\
    501:   0           \text{lb =< input <= ub} \\
    502:   input + lb  \text{input < lb} \\
    503:   \]

    505:   Level: developer

    507: .seealso: `Tao`, `Vec`
    508: @*/
    509: PetscErrorCode TaoSoftThreshold(Vec in, PetscReal lb, PetscReal ub, Vec out)
    510: {
    511:   PetscInt     i, nlocal, mlocal;
    512:   PetscScalar *inarray, *outarray;

    514:   PetscFunctionBegin;
    515:   PetscCall(VecGetArrayPair(in, out, &inarray, &outarray));
    516:   PetscCall(VecGetLocalSize(in, &nlocal));
    517:   PetscCall(VecGetLocalSize(out, &mlocal));

    519:   PetscCheck(nlocal == mlocal, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Input and output vectors need to be of same size.");
    520:   if (lb == ub) {
    521:     PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
    522:     PetscFunctionReturn(PETSC_SUCCESS);
    523:   }
    524:   PetscCheck(lb < ub, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Lower bound needs to be lower than upper bound.");

    526:   if (ub >= 0 && lb < 0) {
    527:     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPN(inarray[i], lb, ub);
    528:   } else if (ub < 0 && lb < 0) {
    529:     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalNN(inarray[i], lb, ub);
    530:   } else {
    531:     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPP(inarray[i], lb, ub);
    532:   }

    534:   PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
    535:   PetscFunctionReturn(PETSC_SUCCESS);
    536: }