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:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

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

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

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

     53: static PetscErrorCode SNESVIMonitorResidual(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
     54: {
     55:   Vec X, F, Finactive;
     56:   IS  isactive;

     58:   PetscFunctionBegin;
     60:   PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
     61:   PetscCall(SNESGetSolution(snes, &X));
     62:   PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
     63:   PetscCall(VecDuplicate(F, &Finactive));
     64:   PetscCall(PetscObjectCompose((PetscObject)Finactive, "__Vec_bc_zero__", (PetscObject)snes));
     65:   PetscCall(VecCopy(F, Finactive));
     66:   PetscCall(VecISSet(Finactive, isactive, 0.0));
     67:   PetscCall(ISDestroy(&isactive));
     68:   PetscCall(PetscViewerPushFormat(vf->viewer, vf->format));
     69:   PetscCall(VecView(Finactive, vf->viewer));
     70:   PetscCall(PetscViewerPopFormat(vf->viewer));
     71:   PetscCall(PetscObjectCompose((PetscObject)Finactive, "__Vec_bc_zero__", NULL));
     72:   PetscCall(VecDestroy(&Finactive));
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: static PetscErrorCode SNESVIMonitorActive(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
     77: {
     78:   Vec X, F, A;
     79:   IS  isactive;

     81:   PetscFunctionBegin;
     83:   PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
     84:   PetscCall(SNESGetSolution(snes, &X));
     85:   PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
     86:   PetscCall(VecDuplicate(F, &A));
     87:   PetscCall(PetscObjectCompose((PetscObject)A, "__Vec_bc_zero__", (PetscObject)snes));
     88:   PetscCall(VecSet(A, 0.));
     89:   PetscCall(VecISSet(A, isactive, 1.));
     90:   PetscCall(ISDestroy(&isactive));
     91:   PetscCall(PetscViewerPushFormat(vf->viewer, vf->format));
     92:   PetscCall(VecView(A, vf->viewer));
     93:   PetscCall(PetscViewerPopFormat(vf->viewer));
     94:   PetscCall(PetscObjectCompose((PetscObject)A, "__Vec_bc_zero__", NULL));
     95:   PetscCall(VecDestroy(&A));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode SNESMonitorVI(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
    100: {
    101:   PetscViewer        viewer = (PetscViewer)dummy;
    102:   const PetscScalar *x, *xl, *xu, *f;
    103:   PetscInt           i, n, act[2] = {0, 0}, fact[2], N;
    104:   /* Number of components that actually hit the bounds (c.f. active variables) */
    105:   PetscInt  act_bound[2] = {0, 0}, fact_bound[2];
    106:   PetscReal rnorm, fnorm, zerotolerance = snes->vizerotolerance;
    107:   double    tmp;

    109:   PetscFunctionBegin;
    111:   PetscCall(VecGetLocalSize(snes->vec_sol, &n));
    112:   PetscCall(VecGetSize(snes->vec_sol, &N));
    113:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    114:   PetscCall(VecGetArrayRead(snes->xu, &xu));
    115:   PetscCall(VecGetArrayRead(snes->vec_sol, &x));
    116:   PetscCall(VecGetArrayRead(snes->vec_func, &f));

    118:   rnorm = 0.0;
    119:   for (i = 0; i < n; i++) {
    120:     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]);
    121:     else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
    122:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
    123:     else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Can never get here");
    124:   }

    126:   for (i = 0; i < n; i++) {
    127:     if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
    128:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
    129:   }
    130:   PetscCall(VecRestoreArrayRead(snes->vec_func, &f));
    131:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    132:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    133:   PetscCall(VecRestoreArrayRead(snes->vec_sol, &x));
    134:   PetscCallMPI(MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
    135:   PetscCallMPI(MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    136:   PetscCallMPI(MPIU_Allreduce(act_bound, fact_bound, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    137:   fnorm = PetscSqrtReal(fnorm);

    139:   PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
    140:   if (snes->ntruebounds) tmp = ((double)(fact[0] + fact[1])) / ((double)snes->ntruebounds);
    141:   else tmp = 0.0;
    142:   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));

    144:   PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
    145:   PetscFunctionReturn(PETSC_SUCCESS);
    146: }

    148: /*
    149:      Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
    150:     || 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
    151:     0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
    152:     for this trick. One assumes that the probability that W is in the null space of J is very, very small.
    153: */
    154: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes, Mat A, Vec F, Vec W, PetscReal fnorm, PetscBool *ismin)
    155: {
    156:   PetscReal a1;
    157:   PetscBool hastranspose;

    159:   PetscFunctionBegin;
    160:   *ismin = PETSC_FALSE;
    161:   PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
    162:   if (hastranspose) {
    163:     /* Compute || J^T F|| */
    164:     PetscCall(MatMultTranspose(A, F, W));
    165:     PetscCall(VecNorm(W, NORM_2, &a1));
    166:     PetscCall(PetscInfo(snes, "|| J^T F|| %g near zero implies found a local minimum\n", (double)(a1 / fnorm)));
    167:     if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
    168:   } else {
    169:     Vec         work;
    170:     PetscScalar result;
    171:     PetscReal   wnorm;

    173:     PetscCall(VecSetRandom(W, NULL));
    174:     PetscCall(VecNorm(W, NORM_2, &wnorm));
    175:     PetscCall(VecDuplicate(W, &work));
    176:     PetscCall(MatMult(A, W, work));
    177:     PetscCall(VecDot(F, work, &result));
    178:     PetscCall(VecDestroy(&work));
    179:     a1 = PetscAbsScalar(result) / (fnorm * wnorm);
    180:     PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n", (double)a1));
    181:     if (a1 < 1.e-4) *ismin = PETSC_TRUE;
    182:   }
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: /*
    187:   SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.

    189:   Notes:
    190:   The convergence criterion currently implemented is
    191:   merit < abstol
    192:   merit < rtol*merit_initial
    193: */
    194: PetscErrorCode SNESConvergedDefault_VI(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gradnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
    195: {
    196:   PetscFunctionBegin;
    198:   PetscAssertPointer(reason, 6);

    200:   *reason = SNES_CONVERGED_ITERATING;

    202:   if (!it) {
    203:     /* set parameter for default relative tolerance convergence test */
    204:     snes->ttol = fnorm * snes->rtol;
    205:   }
    206:   if (fnorm != fnorm) {
    207:     PetscCall(PetscInfo(snes, "Failed to converged, function norm is NaN\n"));
    208:     *reason = SNES_DIVERGED_FUNCTION_NANORINF;
    209:   } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
    210:     PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g\n", (double)fnorm, (double)snes->abstol));
    211:     *reason = SNES_CONVERGED_FNORM_ABS;
    212:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
    213:     PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", snes->nfuncs, snes->max_funcs));
    214:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
    215:   }

    217:   if (it && !*reason) {
    218:     if (fnorm < snes->ttol) {
    219:       PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)snes->ttol));
    220:       *reason = SNES_CONVERGED_FNORM_RELATIVE;
    221:     }
    222:   }
    223:   PetscFunctionReturn(PETSC_SUCCESS);
    224: }

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

    229:    Input Parameters:
    230: .  SNES - nonlinear solver context

    232:    Output Parameters:
    233: .  X - Bound projected X

    235: */

    237: PetscErrorCode SNESVIProjectOntoBounds(SNES snes, Vec X)
    238: {
    239:   const PetscScalar *xl, *xu;
    240:   PetscScalar       *x;
    241:   PetscInt           i, n;

    243:   PetscFunctionBegin;
    244:   PetscCall(VecGetLocalSize(X, &n));
    245:   PetscCall(VecGetArray(X, &x));
    246:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    247:   PetscCall(VecGetArrayRead(snes->xu, &xu));

    249:   for (i = 0; i < n; i++) {
    250:     if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
    251:     else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
    252:   }
    253:   PetscCall(VecRestoreArray(X, &x));
    254:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    255:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    256:   PetscFunctionReturn(PETSC_SUCCESS);
    257: }

    259: /*@
    260:   SNESVIGetActiveSetIS - Gets the global indices for the active set variables

    262:   Input Parameters:
    263: + snes - the `SNES` context
    264: . X    - the `snes` solution vector
    265: - F    - the nonlinear function vector

    267:   Output Parameter:
    268: . ISact - active set index set

    270:   Level: developer

    272:   Note:
    273:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

    275: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`
    276: @*/
    277: PetscErrorCode SNESVIGetActiveSetIS(SNES snes, Vec X, Vec F, IS *ISact)
    278: {
    279:   Vec                Xl = snes->xl, Xu = snes->xu;
    280:   const PetscScalar *x, *f, *xl, *xu;
    281:   PetscInt          *idx_act, i, nlocal, nloc_isact = 0, ilow, ihigh, i1 = 0;
    282:   PetscReal          zerotolerance = snes->vizerotolerance;

    284:   PetscFunctionBegin;
    285:   PetscCall(VecGetLocalSize(X, &nlocal));
    286:   PetscCall(VecGetOwnershipRange(X, &ilow, &ihigh));
    287:   PetscCall(VecGetArrayRead(X, &x));
    288:   PetscCall(VecGetArrayRead(Xl, &xl));
    289:   PetscCall(VecGetArrayRead(Xu, &xu));
    290:   PetscCall(VecGetArrayRead(F, &f));
    291:   /* Compute active set size */
    292:   for (i = 0; i < nlocal; i++) {
    293:     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++;
    294:   }

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

    298:   /* Set active set indices */
    299:   for (i = 0; i < nlocal; i++) {
    300:     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;
    301:   }

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

    306:   PetscCall(VecRestoreArrayRead(X, &x));
    307:   PetscCall(VecRestoreArrayRead(Xl, &xl));
    308:   PetscCall(VecRestoreArrayRead(Xu, &xu));
    309:   PetscCall(VecRestoreArrayRead(F, &f));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }

    313: /*@
    314:   SNESVIComputeInactiveSetFnorm - Computes the function norm for variational inequalities on the inactive set

    316:   Input Parameters:
    317: + snes - the `SNES` context
    318: . F    - the nonlinear function vector
    319: - X    - the `SNES` solution vector

    321:   Output Parameter:
    322: . fnorm - the function norm

    324:   Level: developer

    326:   Note:
    327:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

    329: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESLineSearchSetVIFunctions()`
    330: @*/
    331: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes, Vec F, Vec X, PetscReal *fnorm)
    332: {
    333:   const PetscScalar *x, *xl, *xu, *f;
    334:   PetscInt           i, n;
    335:   PetscReal          zerotolerance = snes->vizerotolerance;

    337:   PetscFunctionBegin;
    339:   PetscAssertPointer(fnorm, 4);
    340:   PetscCall(VecGetLocalSize(X, &n));
    341:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    342:   PetscCall(VecGetArrayRead(snes->xu, &xu));
    343:   PetscCall(VecGetArrayRead(X, &x));
    344:   PetscCall(VecGetArrayRead(F, &f));
    345:   *fnorm = 0.0;
    346:   for (i = 0; i < n; i++) {
    347:     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)) *fnorm += PetscRealPart(PetscConj(f[i]) * f[i]);
    348:   }
    349:   PetscCall(VecRestoreArrayRead(F, &f));
    350:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    351:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    352:   PetscCall(VecRestoreArrayRead(X, &x));
    353:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
    354:   *fnorm = PetscSqrtReal(*fnorm);
    355:   PetscFunctionReturn(PETSC_SUCCESS);
    356: }

    358: /*@
    359:   SNESVIComputeInactiveSetFtY - Computes the directional derivative for variational inequalities on the inactive set,
    360:   assuming that there exists some $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$ (relevant for some line search algorithms)

    362:   Input Parameters:
    363: + snes - the `SNES` context
    364: . F    - the nonlinear function vector
    365: . X    - the `SNES` solution vector
    366: - Y    - the direction vector

    368:   Output Parameter:
    369: . fty - the directional derivative

    371:   Level: developer

    373:   Note:
    374:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

    376: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`
    377: @*/
    378: PetscErrorCode SNESVIComputeInactiveSetFtY(SNES snes, Vec F, Vec X, Vec Y, PetscScalar *fty)
    379: {
    380:   const PetscScalar *x, *xl, *xu, *y, *f;
    381:   PetscInt           i, n;
    382:   PetscReal          zerotolerance = snes->vizerotolerance;

    384:   PetscFunctionBegin;
    386:   PetscAssertPointer(fty, 5);
    387:   PetscCall(VecGetLocalSize(X, &n));
    388:   PetscCall(VecGetArrayRead(F, &f));
    389:   PetscCall(VecGetArrayRead(X, &x));
    390:   PetscCall(VecGetArrayRead(snes->xl, &xl));
    391:   PetscCall(VecGetArrayRead(snes->xu, &xu));
    392:   PetscCall(VecGetArrayRead(Y, &y));
    393:   *fty = 0.0;
    394:   for (i = 0; i < n; i++) {
    395:     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)) *fty += f[i] * PetscConj(y[i]);
    396:   }
    397:   PetscCall(VecRestoreArrayRead(F, &f));
    398:   PetscCall(VecRestoreArrayRead(X, &x));
    399:   PetscCall(VecRestoreArrayRead(snes->xl, &xl));
    400:   PetscCall(VecRestoreArrayRead(snes->xu, &xu));
    401:   PetscCall(VecRestoreArrayRead(Y, &y));
    402:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fty, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
    403:   PetscFunctionReturn(PETSC_SUCCESS);
    404: }

    406: static PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes, Vec xl, Vec xu)
    407: {
    408:   PetscFunctionBegin;
    409:   PetscCall(DMComputeVariableBounds(snes->dm, xl, xu));
    410:   PetscFunctionReturn(PETSC_SUCCESS);
    411: }

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

    417:    Input Parameter:
    418: .  snes - the SNES context

    420:    Application Interface Routine: SNESSetUp()

    422:    Notes:
    423:    For basic use of the SNES solvers, the user need not explicitly call
    424:    SNESSetUp(), since these actions will automatically occur during
    425:    the call to SNESSolve().
    426:  */
    427: PetscErrorCode SNESSetUp_VI(SNES snes)
    428: {
    429:   PetscInt i_start[3], i_end[3];

    431:   PetscFunctionBegin;
    432:   PetscCall(SNESSetWorkVecs(snes, 1));
    433:   PetscCall(SNESSetUpMatrices(snes));

    435:   if (!snes->ops->computevariablebounds && snes->dm) {
    436:     PetscBool flag;
    437:     PetscCall(DMHasVariableBounds(snes->dm, &flag));
    438:     if (flag) snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
    439:   }
    440:   if (!snes->usersetbounds) {
    441:     if (snes->ops->computevariablebounds) {
    442:       if (!snes->xl) PetscCall(VecDuplicate(snes->work[0], &snes->xl));
    443:       if (!snes->xu) PetscCall(VecDuplicate(snes->work[0], &snes->xu));
    444:       PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
    445:     } else if (!snes->xl && !snes->xu) {
    446:       /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
    447:       PetscCall(VecDuplicate(snes->work[0], &snes->xl));
    448:       PetscCall(VecSet(snes->xl, PETSC_NINFINITY));
    449:       PetscCall(VecDuplicate(snes->work[0], &snes->xu));
    450:       PetscCall(VecSet(snes->xu, PETSC_INFINITY));
    451:     } else {
    452:       /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
    453:       PetscCall(VecGetOwnershipRange(snes->work[0], i_start, i_end));
    454:       PetscCall(VecGetOwnershipRange(snes->xl, i_start + 1, i_end + 1));
    455:       PetscCall(VecGetOwnershipRange(snes->xu, i_start + 2, i_end + 2));
    456:       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]))
    457:         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.");
    458:     }
    459:   }
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }
    462: PetscErrorCode SNESReset_VI(SNES snes)
    463: {
    464:   PetscFunctionBegin;
    465:   PetscCall(VecDestroy(&snes->xl));
    466:   PetscCall(VecDestroy(&snes->xu));
    467:   snes->usersetbounds = PETSC_FALSE;
    468:   PetscFunctionReturn(PETSC_SUCCESS);
    469: }

    471: /*
    472:    SNESDestroy_VI - Destroys the private SNES_VI context that was created
    473:    with SNESCreate_VI().

    475:    Input Parameter:
    476: .  snes - the SNES context

    478:    Application Interface Routine: SNESDestroy()
    479:  */
    480: PetscErrorCode SNESDestroy_VI(SNES snes)
    481: {
    482:   PetscFunctionBegin;
    483:   PetscCall(PetscFree(snes->data));

    485:   /* clear composed functions */
    486:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", NULL));
    487:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", NULL));
    488:   PetscFunctionReturn(PETSC_SUCCESS);
    489: }

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

    495:   Input Parameters:
    496: + snes - the `SNES` context.
    497: . xl   - lower bound.
    498: - xu   - upper bound.

    500:   Level: advanced

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

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

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

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

    513:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

    515: .seealso: [](sec_vi), `SNES`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
    516: @*/
    517: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
    518: {
    519:   PetscErrorCode (*f)(SNES, Vec, Vec);

    521:   PetscFunctionBegin;
    525:   PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetVariableBounds_C", &f));
    526:   if (f) PetscUseMethod(snes, "SNESVISetVariableBounds_C", (SNES, Vec, Vec), (snes, xl, xu));
    527:   else PetscCall(SNESVISetVariableBounds_VI(snes, xl, xu));
    528:   snes->usersetbounds = PETSC_TRUE;
    529:   PetscFunctionReturn(PETSC_SUCCESS);
    530: }

    532: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes, Vec xl, Vec xu)
    533: {
    534:   const PetscScalar *xxl, *xxu;
    535:   PetscInt           i, n, cnt = 0;

    537:   PetscFunctionBegin;
    538:   PetscCall(SNESGetFunction(snes, &snes->vec_func, NULL, NULL));
    539:   PetscCheck(snes->vec_func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
    540:   {
    541:     PetscInt xlN, xuN, N;
    542:     PetscCall(VecGetSize(xl, &xlN));
    543:     PetscCall(VecGetSize(xu, &xuN));
    544:     PetscCall(VecGetSize(snes->vec_func, &N));
    545:     PetscCheck(xlN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths lower bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xlN, N);
    546:     PetscCheck(xuN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths: upper bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xuN, N);
    547:   }
    548:   PetscCall(PetscObjectReference((PetscObject)xl));
    549:   PetscCall(PetscObjectReference((PetscObject)xu));
    550:   PetscCall(VecDestroy(&snes->xl));
    551:   PetscCall(VecDestroy(&snes->xu));
    552:   snes->xl = xl;
    553:   snes->xu = xu;
    554:   PetscCall(VecGetLocalSize(xl, &n));
    555:   PetscCall(VecGetArrayRead(xl, &xxl));
    556:   PetscCall(VecGetArrayRead(xu, &xxu));
    557:   for (i = 0; i < n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));

    559:   PetscCallMPI(MPIU_Allreduce(&cnt, &snes->ntruebounds, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
    560:   PetscCall(VecRestoreArrayRead(xl, &xxl));
    561:   PetscCall(VecRestoreArrayRead(xu, &xxu));
    562:   PetscFunctionReturn(PETSC_SUCCESS);
    563: }

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

    569:   Input Parameters:
    570: + snes - the `SNES` context.
    571: . xl   - lower bound (may be `NULL`)
    572: - xu   - upper bound (may be `NULL`)

    574:   Level: advanced

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

    579: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
    580: @*/
    581: PetscErrorCode SNESVIGetVariableBounds(SNES snes, Vec *xl, Vec *xu)
    582: {
    583:   PetscFunctionBegin;
    584:   PetscCheck(snes->usersetbounds, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must set SNESVI bounds before calling SNESVIGetVariableBounds()");
    585:   if (xl) *xl = snes->xl;
    586:   if (xu) *xu = snes->xu;
    587:   PetscFunctionReturn(PETSC_SUCCESS);
    588: }

    590: PetscErrorCode SNESSetFromOptions_VI(SNES snes, PetscOptionItems PetscOptionsObject)
    591: {
    592:   PetscBool flg = PETSC_FALSE;

    594:   PetscFunctionBegin;
    595:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES VI options");
    596:   PetscCall(PetscOptionsReal("-snes_vi_zero_tolerance", "Tolerance for considering x[] value to be on a bound", "None", snes->vizerotolerance, &snes->vizerotolerance, NULL));
    597:   PetscCall(PetscOptionsBool("-snes_vi_monitor", "Monitor all non-active variables", "SNESMonitorResidual", flg, &flg, NULL));
    598:   if (flg) PetscCall(SNESMonitorSet(snes, SNESMonitorVI, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)), NULL));
    599:   flg = PETSC_FALSE;
    600:   PetscCall(SNESMonitorSetFromOptions(snes, "-snes_vi_monitor_residual", "View residual at each iteration, using zero for active constraints", "SNESVIMonitorResidual", SNESVIMonitorResidual, NULL));
    601:   PetscCall(SNESMonitorSetFromOptions(snes, "-snes_vi_monitor_active", "View active set at each iteration, using zero for inactive dofs", "SNESVIMonitorActive", SNESVIMonitorActive, NULL));
    602:   PetscOptionsHeadEnd();
    603:   PetscFunctionReturn(PETSC_SUCCESS);
    604: }