Actual source code: vi.c

      1: #include <petsc/private/snesimpl.h>
      2: #include <petscdm.h>

      4: /*@C
      5:   SNESVISetComputeVariableBounds - Sets a function that is called to compute the bounds on variable for
      6:   (differential) variable inequalities.

      8:   Input Parameters:
      9: + snes    - the `SNES` context
     10: - compute - function that computes the bounds

     12:   Calling sequence of `compute`:
     13: + snes   - the `SNES` context
     14: . lower  - vector to hold lower bounds
     15: - higher - vector to hold upper bounds

     17:   Level: advanced

     19:   Notes:
     20:   Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.

     22:   For entries with no bounds you can set `PETSC_NINFINITY` or `PETSC_INFINITY`

     24:   You may use `SNESVISetVariableBounds()` to provide the bounds once if they will never change

     26:   If you have associated a `DM` with the `SNES` and provided a function to the `DM` via `DMSetVariableBounds()` that will be used automatically
     27:   to provide the bounds and you need not use this function.

     29: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `DMSetVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`,
     30:           `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
     31: @*/
     32: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES snes, Vec lower, Vec higher))
     33: {
     34:   PetscErrorCode (*f)(SNES, PetscErrorCode (*)(SNES, Vec, Vec));

     36:   PetscFunctionBegin;
     38:   PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", &f));
     39:   if (f) PetscUseMethod(snes, "SNESVISetComputeVariableBounds_C", (SNES, PetscErrorCode (*)(SNES, Vec, Vec)), (snes, compute));
     40:   else PetscCall(SNESVISetComputeVariableBounds_VI(snes, compute));
     41:   PetscFunctionReturn(PETSC_SUCCESS);
     42: }

     44: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes, SNESVIComputeVariableBoundsFn *compute)
     45: {
     46:   PetscFunctionBegin;
     47:   snes->ops->computevariablebounds = compute;
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: static PetscErrorCode SNESVIMonitorResidual(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
     52: {
     53:   Vec         X, F, Finactive;
     54:   IS          isactive;
     55:   PetscViewer viewer = (PetscViewer)dummy;

     57:   PetscFunctionBegin;
     58:   PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
     59:   PetscCall(SNESGetSolution(snes, &X));
     60:   PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
     61:   PetscCall(VecDuplicate(F, &Finactive));
     62:   PetscCall(VecCopy(F, Finactive));
     63:   PetscCall(VecISSet(Finactive, isactive, 0.0));
     64:   PetscCall(ISDestroy(&isactive));
     65:   PetscCall(VecView(Finactive, viewer));
     66:   PetscCall(VecDestroy(&Finactive));
     67:   PetscFunctionReturn(PETSC_SUCCESS);
     68: }

     70: static PetscErrorCode SNESMonitorVI(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
     71: {
     72:   PetscViewer        viewer = (PetscViewer)dummy;
     73:   const PetscScalar *x, *xl, *xu, *f;
     74:   PetscInt           i, n, act[2] = {0, 0}, fact[2], N;
     75:   /* Number of components that actually hit the bounds (c.f. active variables) */
     76:   PetscInt  act_bound[2] = {0, 0}, fact_bound[2];
     77:   PetscReal rnorm, fnorm, zerotolerance = snes->vizerotolerance;
     78:   double    tmp;

     80:   PetscFunctionBegin;
     82:   PetscCall(VecGetLocalSize(snes->vec_sol, &n));
     83:   PetscCall(VecGetSize(snes->vec_sol, &N));
     84:   PetscCall(VecGetArrayRead(snes->xl, &xl));
     85:   PetscCall(VecGetArrayRead(snes->xu, &xu));
     86:   PetscCall(VecGetArrayRead(snes->vec_sol, &x));
     87:   PetscCall(VecGetArrayRead(snes->vec_func, &f));

     89:   rnorm = 0.0;
     90:   for (i = 0; i < n; i++) {
     91:     if ((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0)) rnorm += PetscRealPart(PetscConj(f[i]) * f[i]);
     92:     else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
     93:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
     94:     else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Can never get here");
     95:   }

     97:   for (i = 0; i < n; i++) {
     98:     if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
     99:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
    100:   }
    101:   PetscCall(VecRestoreArrayRead(snes->vec_func, &f));
    102:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    103:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    104:   PetscCall(VecRestoreArrayRead(snes->vec_sol, &x));
    105:   PetscCallMPI(MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
    106:   PetscCallMPI(MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    107:   PetscCallMPI(MPIU_Allreduce(act_bound, fact_bound, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    108:   fnorm = PetscSqrtReal(fnorm);

    110:   PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
    111:   if (snes->ntruebounds) tmp = ((double)(fact[0] + fact[1])) / ((double)snes->ntruebounds);
    112:   else tmp = 0.0;
    113:   PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES VI Function norm %g Active lower constraints %" PetscInt_FMT "/%" PetscInt_FMT " upper constraints %" PetscInt_FMT "/%" PetscInt_FMT " Percent of total %g Percent of bounded %g\n", its, (double)fnorm, fact[0], fact_bound[0], fact[1], fact_bound[1], ((double)(fact[0] + fact[1])) / ((double)N), tmp));

    115:   PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: /*
    120:      Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
    121:     || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that
    122:     0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
    123:     for this trick. One assumes that the probability that W is in the null space of J is very, very small.
    124: */
    125: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes, Mat A, Vec F, Vec W, PetscReal fnorm, PetscBool *ismin)
    126: {
    127:   PetscReal a1;
    128:   PetscBool hastranspose;

    130:   PetscFunctionBegin;
    131:   *ismin = PETSC_FALSE;
    132:   PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
    133:   if (hastranspose) {
    134:     /* Compute || J^T F|| */
    135:     PetscCall(MatMultTranspose(A, F, W));
    136:     PetscCall(VecNorm(W, NORM_2, &a1));
    137:     PetscCall(PetscInfo(snes, "|| J^T F|| %g near zero implies found a local minimum\n", (double)(a1 / fnorm)));
    138:     if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
    139:   } else {
    140:     Vec         work;
    141:     PetscScalar result;
    142:     PetscReal   wnorm;

    144:     PetscCall(VecSetRandom(W, NULL));
    145:     PetscCall(VecNorm(W, NORM_2, &wnorm));
    146:     PetscCall(VecDuplicate(W, &work));
    147:     PetscCall(MatMult(A, W, work));
    148:     PetscCall(VecDot(F, work, &result));
    149:     PetscCall(VecDestroy(&work));
    150:     a1 = PetscAbsScalar(result) / (fnorm * wnorm);
    151:     PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n", (double)a1));
    152:     if (a1 < 1.e-4) *ismin = PETSC_TRUE;
    153:   }
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: /*
    158:   SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.

    160:   Notes:
    161:   The convergence criterion currently implemented is
    162:   merit < abstol
    163:   merit < rtol*merit_initial
    164: */
    165: PetscErrorCode SNESConvergedDefault_VI(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gradnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
    166: {
    167:   PetscFunctionBegin;
    169:   PetscAssertPointer(reason, 6);

    171:   *reason = SNES_CONVERGED_ITERATING;

    173:   if (!it) {
    174:     /* set parameter for default relative tolerance convergence test */
    175:     snes->ttol = fnorm * snes->rtol;
    176:   }
    177:   if (fnorm != fnorm) {
    178:     PetscCall(PetscInfo(snes, "Failed to converged, function norm is NaN\n"));
    179:     *reason = SNES_DIVERGED_FNORM_NAN;
    180:   } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
    181:     PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g\n", (double)fnorm, (double)snes->abstol));
    182:     *reason = SNES_CONVERGED_FNORM_ABS;
    183:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
    184:     PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", snes->nfuncs, snes->max_funcs));
    185:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
    186:   }

    188:   if (it && !*reason) {
    189:     if (fnorm < snes->ttol) {
    190:       PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)snes->ttol));
    191:       *reason = SNES_CONVERGED_FNORM_RELATIVE;
    192:     }
    193:   }
    194:   PetscFunctionReturn(PETSC_SUCCESS);
    195: }

    197: /*
    198:    SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.

    200:    Input Parameters:
    201: .  SNES - nonlinear solver context

    203:    Output Parameters:
    204: .  X - Bound projected X

    206: */

    208: PetscErrorCode SNESVIProjectOntoBounds(SNES snes, Vec X)
    209: {
    210:   const PetscScalar *xl, *xu;
    211:   PetscScalar       *x;
    212:   PetscInt           i, n;

    214:   PetscFunctionBegin;
    215:   PetscCall(VecGetLocalSize(X, &n));
    216:   PetscCall(VecGetArray(X, &x));
    217:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    218:   PetscCall(VecGetArrayRead(snes->xu, &xu));

    220:   for (i = 0; i < n; i++) {
    221:     if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
    222:     else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
    223:   }
    224:   PetscCall(VecRestoreArray(X, &x));
    225:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    226:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: /*@
    231:   SNESVIGetActiveSetIS - Gets the global indices for the active set variables

    233:   Input Parameters:
    234: + snes - the `SNES` context
    235: . X    - the `snes` solution vector
    236: - F    - the nonlinear function vector

    238:   Output Parameter:
    239: . ISact - active set index set

    241:   Level: developer

    243: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`
    244: @*/
    245: PetscErrorCode SNESVIGetActiveSetIS(SNES snes, Vec X, Vec F, IS *ISact)
    246: {
    247:   Vec                Xl = snes->xl, Xu = snes->xu;
    248:   const PetscScalar *x, *f, *xl, *xu;
    249:   PetscInt          *idx_act, i, nlocal, nloc_isact = 0, ilow, ihigh, i1 = 0;
    250:   PetscReal          zerotolerance = snes->vizerotolerance;

    252:   PetscFunctionBegin;
    253:   PetscCall(VecGetLocalSize(X, &nlocal));
    254:   PetscCall(VecGetOwnershipRange(X, &ilow, &ihigh));
    255:   PetscCall(VecGetArrayRead(X, &x));
    256:   PetscCall(VecGetArrayRead(Xl, &xl));
    257:   PetscCall(VecGetArrayRead(Xu, &xu));
    258:   PetscCall(VecGetArrayRead(F, &f));
    259:   /* Compute active set size */
    260:   for (i = 0; i < nlocal; i++) {
    261:     if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) nloc_isact++;
    262:   }

    264:   PetscCall(PetscMalloc1(nloc_isact, &idx_act));

    266:   /* Set active set indices */
    267:   for (i = 0; i < nlocal; i++) {
    268:     if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) idx_act[i1++] = ilow + i;
    269:   }

    271:   /* Create active set IS */
    272:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), nloc_isact, idx_act, PETSC_OWN_POINTER, ISact));

    274:   PetscCall(VecRestoreArrayRead(X, &x));
    275:   PetscCall(VecRestoreArrayRead(Xl, &xl));
    276:   PetscCall(VecRestoreArrayRead(Xu, &xu));
    277:   PetscCall(VecRestoreArrayRead(F, &f));
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }

    281: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes, Vec F, Vec X, PetscReal *fnorm)
    282: {
    283:   const PetscScalar *x, *xl, *xu, *f;
    284:   PetscInt           i, n;
    285:   PetscReal          rnorm, zerotolerance = snes->vizerotolerance;

    287:   PetscFunctionBegin;
    288:   PetscCall(VecGetLocalSize(X, &n));
    289:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    290:   PetscCall(VecGetArrayRead(snes->xu, &xu));
    291:   PetscCall(VecGetArrayRead(X, &x));
    292:   PetscCall(VecGetArrayRead(F, &f));
    293:   rnorm = 0.0;
    294:   for (i = 0; i < n; i++) {
    295:     if ((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0)) rnorm += PetscRealPart(PetscConj(f[i]) * f[i]);
    296:   }
    297:   PetscCall(VecRestoreArrayRead(F, &f));
    298:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    299:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    300:   PetscCall(VecRestoreArrayRead(X, &x));
    301:   PetscCallMPI(MPIU_Allreduce(&rnorm, fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
    302:   *fnorm = PetscSqrtReal(*fnorm);
    303:   PetscFunctionReturn(PETSC_SUCCESS);
    304: }

    306: static PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes, Vec xl, Vec xu)
    307: {
    308:   PetscFunctionBegin;
    309:   PetscCall(DMComputeVariableBounds(snes->dm, xl, xu));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: /*
    314:    SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
    315:    of the SNESVI nonlinear solver.

    317:    Input Parameter:
    318: .  snes - the SNES context

    320:    Application Interface Routine: SNESSetUp()

    322:    Notes:
    323:    For basic use of the SNES solvers, the user need not explicitly call
    324:    SNESSetUp(), since these actions will automatically occur during
    325:    the call to SNESSolve().
    326:  */
    327: PetscErrorCode SNESSetUp_VI(SNES snes)
    328: {
    329:   PetscInt i_start[3], i_end[3];

    331:   PetscFunctionBegin;
    332:   PetscCall(SNESSetWorkVecs(snes, 1));
    333:   PetscCall(SNESSetUpMatrices(snes));

    335:   if (!snes->ops->computevariablebounds && snes->dm) {
    336:     PetscBool flag;
    337:     PetscCall(DMHasVariableBounds(snes->dm, &flag));
    338:     if (flag) snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
    339:   }
    340:   if (!snes->usersetbounds) {
    341:     if (snes->ops->computevariablebounds) {
    342:       if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
    343:       if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
    344:       PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
    345:     } else if (!snes->xl && !snes->xu) {
    346:       /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
    347:       PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
    348:       PetscCall(VecSet(snes->xl, PETSC_NINFINITY));
    349:       PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
    350:       PetscCall(VecSet(snes->xu, PETSC_INFINITY));
    351:     } else {
    352:       /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
    353:       PetscCall(VecGetOwnershipRange(snes->vec_sol, i_start, i_end));
    354:       PetscCall(VecGetOwnershipRange(snes->xl, i_start + 1, i_end + 1));
    355:       PetscCall(VecGetOwnershipRange(snes->xu, i_start + 2, i_end + 2));
    356:       if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2]))
    357:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Distribution of lower bound, upper bound and the solution vector should be identical across all the processors.");
    358:     }
    359:   }
    360:   PetscFunctionReturn(PETSC_SUCCESS);
    361: }
    362: PetscErrorCode SNESReset_VI(SNES snes)
    363: {
    364:   PetscFunctionBegin;
    365:   PetscCall(VecDestroy(&snes->xl));
    366:   PetscCall(VecDestroy(&snes->xu));
    367:   snes->usersetbounds = PETSC_FALSE;
    368:   PetscFunctionReturn(PETSC_SUCCESS);
    369: }

    371: /*
    372:    SNESDestroy_VI - Destroys the private SNES_VI context that was created
    373:    with SNESCreate_VI().

    375:    Input Parameter:
    376: .  snes - the SNES context

    378:    Application Interface Routine: SNESDestroy()
    379:  */
    380: PetscErrorCode SNESDestroy_VI(SNES snes)
    381: {
    382:   PetscFunctionBegin;
    383:   PetscCall(PetscFree(snes->data));

    385:   /* clear composed functions */
    386:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", NULL));
    387:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", NULL));
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

    391: /*@
    392:   SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. `xl` <= x <= `xu`. This allows solving
    393:   (differential) variable inequalities.

    395:   Input Parameters:
    396: + snes - the `SNES` context.
    397: . xl   - lower bound.
    398: - xu   - upper bound.

    400:   Level: advanced

    402:   Notes:
    403:   If this routine is not called then the lower and upper bounds are set to
    404:   `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.

    406:   Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS` or semi-smooth `SNESVINEWTONSSLS` solvers.

    408:   For particular components that have no bounds you can use `PETSC_NINFINITY` or `PETSC_INFINITY`

    410:   `SNESVISetComputeVariableBounds()` can be used to provide a function that computes the bounds. This should be used if you are using, for example, grid
    411:   sequencing and need bounds set for a variety of vectors

    413: .seealso: [](sec_vi), `SNES`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
    414: @*/
    415: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
    416: {
    417:   PetscErrorCode (*f)(SNES, Vec, Vec);

    419:   PetscFunctionBegin;
    423:   PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetVariableBounds_C", &f));
    424:   if (f) PetscUseMethod(snes, "SNESVISetVariableBounds_C", (SNES, Vec, Vec), (snes, xl, xu));
    425:   else PetscCall(SNESVISetVariableBounds_VI(snes, xl, xu));
    426:   snes->usersetbounds = PETSC_TRUE;
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes, Vec xl, Vec xu)
    431: {
    432:   const PetscScalar *xxl, *xxu;
    433:   PetscInt           i, n, cnt = 0;

    435:   PetscFunctionBegin;
    436:   PetscCall(SNESGetFunction(snes, &snes->vec_func, NULL, NULL));
    437:   PetscCheck(snes->vec_func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
    438:   {
    439:     PetscInt xlN, xuN, N;
    440:     PetscCall(VecGetSize(xl, &xlN));
    441:     PetscCall(VecGetSize(xu, &xuN));
    442:     PetscCall(VecGetSize(snes->vec_func, &N));
    443:     PetscCheck(xlN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths lower bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xlN, N);
    444:     PetscCheck(xuN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths: upper bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xuN, N);
    445:   }
    446:   PetscCall(PetscObjectReference((PetscObject)xl));
    447:   PetscCall(PetscObjectReference((PetscObject)xu));
    448:   PetscCall(VecDestroy(&snes->xl));
    449:   PetscCall(VecDestroy(&snes->xu));
    450:   snes->xl = xl;
    451:   snes->xu = xu;
    452:   PetscCall(VecGetLocalSize(xl, &n));
    453:   PetscCall(VecGetArrayRead(xl, &xxl));
    454:   PetscCall(VecGetArrayRead(xu, &xxu));
    455:   for (i = 0; i < n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));

    457:   PetscCallMPI(MPIU_Allreduce(&cnt, &snes->ntruebounds, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    458:   PetscCall(VecRestoreArrayRead(xl, &xxl));
    459:   PetscCall(VecRestoreArrayRead(xu, &xxu));
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }

    463: /*@
    464:   SNESVIGetVariableBounds - Gets the lower and upper bounds for the solution vector. `xl` <= x <= `xu`. These are used in solving
    465:   (differential) variable inequalities.

    467:   Input Parameters:
    468: + snes - the `SNES` context.
    469: . xl   - lower bound (may be `NULL`)
    470: - xu   - upper bound (may be `NULL`)

    472:   Level: advanced

    474:   Note:
    475:   These vectors are owned by the `SNESVI` and should not be destroyed by the caller

    477: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
    478: @*/
    479: PetscErrorCode SNESVIGetVariableBounds(SNES snes, Vec *xl, Vec *xu)
    480: {
    481:   PetscFunctionBegin;
    482:   PetscCheck(snes->usersetbounds, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must set SNESVI bounds before calling SNESVIGetVariableBounds()");
    483:   if (xl) *xl = snes->xl;
    484:   if (xu) *xu = snes->xu;
    485:   PetscFunctionReturn(PETSC_SUCCESS);
    486: }

    488: PetscErrorCode SNESSetFromOptions_VI(SNES snes, PetscOptionItems *PetscOptionsObject)
    489: {
    490:   PetscBool flg = PETSC_FALSE;

    492:   PetscFunctionBegin;
    493:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES VI options");
    494:   PetscCall(PetscOptionsReal("-snes_vi_zero_tolerance", "Tolerance for considering x[] value to be on a bound", "None", snes->vizerotolerance, &snes->vizerotolerance, NULL));
    495:   PetscCall(PetscOptionsBool("-snes_vi_monitor", "Monitor all non-active variables", "SNESMonitorResidual", flg, &flg, NULL));
    496:   if (flg) PetscCall(SNESMonitorSet(snes, SNESMonitorVI, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)), NULL));
    497:   flg = PETSC_FALSE;
    498:   PetscCall(PetscOptionsBool("-snes_vi_monitor_residual", "Monitor residual all non-active variables; using zero for active constraints", "SNESMonitorVIResidual", flg, &flg, NULL));
    499:   if (flg) PetscCall(SNESMonitorSet(snes, SNESVIMonitorResidual, PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)), NULL));
    500:   PetscOptionsHeadEnd();
    501:   PetscFunctionReturn(PETSC_SUCCESS);
    502: }