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: }