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, SNESVIComputeVariableBoundsFunction 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:   PetscCall(MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
106:   PetscCall(MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
107:   PetscCall(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:   PetscCall(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:   PetscCall(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: }