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 parameter:
9: + snes - the `SNES` context
10: - compute - function that computes the bounds
12: Calling Sequence of `compute`:
13: $ PetscErrorCode compute(SNES snes, Vec lower, Vec higher)
14: + snes - the `SNES` context
15: . lower - vector to hold lower bounds
16: - higher - vector to hold upper bounds
18: Level: advanced
20: Notes:
21: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.
23: For entries with no bounds you can set `PETSC_NINFINITY` or `PETSC_INFINITY`
25: You may use `SNESVISetVariableBounds()` to provide the bounds once if they will never change
27: If you have associated a `DM` with the `SNES` and provided a function to the `DM` via `DMSetVariableBounds()` that will be used automatically
28: to provide the bounds and you need not use this function.
30: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `DMSetVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`,
31: 'SNESSetType()`
32: @*/
33: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES, Vec, Vec))
34: {
35: PetscErrorCode (*f)(SNES, PetscErrorCode (*)(SNES, Vec, Vec));
37: PetscFunctionBegin;
39: PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", &f));
40: if (f) PetscUseMethod(snes, "SNESVISetComputeVariableBounds_C", (SNES, PetscErrorCode(*)(SNES, Vec, Vec)), (snes, compute));
41: else PetscCall(SNESVISetComputeVariableBounds_VI(snes, compute));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes, SNESVIComputeVariableBoundsFunction compute)
46: {
47: PetscFunctionBegin;
48: snes->ops->computevariablebounds = compute;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: PetscErrorCode SNESVIMonitorResidual(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
53: {
54: Vec X, F, Finactive;
55: IS isactive;
56: PetscViewer viewer = (PetscViewer)dummy;
58: PetscFunctionBegin;
59: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
60: PetscCall(SNESGetSolution(snes, &X));
61: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
62: PetscCall(VecDuplicate(F, &Finactive));
63: PetscCall(VecCopy(F, Finactive));
64: PetscCall(VecISSet(Finactive, isactive, 0.0));
65: PetscCall(ISDestroy(&isactive));
66: PetscCall(VecView(Finactive, viewer));
67: PetscCall(VecDestroy(&Finactive));
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: PetscErrorCode SNESMonitorVI(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
72: {
73: PetscViewer viewer = (PetscViewer)dummy;
74: const PetscScalar *x, *xl, *xu, *f;
75: PetscInt i, n, act[2] = {0, 0}, fact[2], N;
76: /* Number of components that actually hit the bounds (c.f. active variables) */
77: PetscInt act_bound[2] = {0, 0}, fact_bound[2];
78: PetscReal rnorm, fnorm, zerotolerance = snes->vizerotolerance;
79: double tmp;
81: PetscFunctionBegin;
83: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
84: PetscCall(VecGetSize(snes->vec_sol, &N));
85: PetscCall(VecGetArrayRead(snes->xl, &xl));
86: PetscCall(VecGetArrayRead(snes->xu, &xu));
87: PetscCall(VecGetArrayRead(snes->vec_sol, &x));
88: PetscCall(VecGetArrayRead(snes->vec_func, &f));
90: rnorm = 0.0;
91: for (i = 0; i < n; i++) {
92: 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]);
93: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
94: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
95: else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Can never get here");
96: }
98: for (i = 0; i < n; i++) {
99: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
100: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
101: }
102: PetscCall(VecRestoreArrayRead(snes->vec_func, &f));
103: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
104: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
105: PetscCall(VecRestoreArrayRead(snes->vec_sol, &x));
106: PetscCall(MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
107: PetscCall(MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
108: PetscCall(MPIU_Allreduce(act_bound, fact_bound, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
109: fnorm = PetscSqrtReal(fnorm);
111: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
112: if (snes->ntruebounds) tmp = ((double)(fact[0] + fact[1])) / ((double)snes->ntruebounds);
113: else tmp = 0.0;
114: 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));
116: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*
121: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
122: || 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
123: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
124: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
125: */
126: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes, Mat A, Vec F, Vec W, PetscReal fnorm, PetscBool *ismin)
127: {
128: PetscReal a1;
129: PetscBool hastranspose;
131: PetscFunctionBegin;
132: *ismin = PETSC_FALSE;
133: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
134: if (hastranspose) {
135: /* Compute || J^T F|| */
136: PetscCall(MatMultTranspose(A, F, W));
137: PetscCall(VecNorm(W, NORM_2, &a1));
138: PetscCall(PetscInfo(snes, "|| J^T F|| %g near zero implies found a local minimum\n", (double)(a1 / fnorm)));
139: if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
140: } else {
141: Vec work;
142: PetscScalar result;
143: PetscReal wnorm;
145: PetscCall(VecSetRandom(W, NULL));
146: PetscCall(VecNorm(W, NORM_2, &wnorm));
147: PetscCall(VecDuplicate(W, &work));
148: PetscCall(MatMult(A, W, work));
149: PetscCall(VecDot(F, work, &result));
150: PetscCall(VecDestroy(&work));
151: a1 = PetscAbsScalar(result) / (fnorm * wnorm);
152: PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n", (double)a1));
153: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
154: }
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: /*
159: Checks if J^T(F - J*X) = 0
160: */
161: PetscErrorCode SNESVICheckResidual_Private(SNES snes, Mat A, Vec F, Vec X, Vec W1, Vec W2)
162: {
163: PetscReal a1, a2;
164: PetscBool hastranspose;
166: PetscFunctionBegin;
167: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
168: if (hastranspose) {
169: PetscCall(MatMult(A, X, W1));
170: PetscCall(VecAXPY(W1, -1.0, F));
172: /* Compute || J^T W|| */
173: PetscCall(MatMultTranspose(A, W1, W2));
174: PetscCall(VecNorm(W1, NORM_2, &a1));
175: PetscCall(VecNorm(W2, NORM_2, &a2));
176: if (a1 != 0.0) PetscCall(PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n", (double)(a2 / a1)));
177: }
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*
182: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
184: Notes:
185: The convergence criterion currently implemented is
186: merit < abstol
187: merit < rtol*merit_initial
188: */
189: PetscErrorCode SNESConvergedDefault_VI(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gradnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
190: {
191: PetscFunctionBegin;
195: *reason = SNES_CONVERGED_ITERATING;
197: if (!it) {
198: /* set parameter for default relative tolerance convergence test */
199: snes->ttol = fnorm * snes->rtol;
200: }
201: if (fnorm != fnorm) {
202: PetscCall(PetscInfo(snes, "Failed to converged, function norm is NaN\n"));
203: *reason = SNES_DIVERGED_FNORM_NAN;
204: } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
205: PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g\n", (double)fnorm, (double)snes->abstol));
206: *reason = SNES_CONVERGED_FNORM_ABS;
207: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
208: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", snes->nfuncs, snes->max_funcs));
209: *reason = SNES_DIVERGED_FUNCTION_COUNT;
210: }
212: if (it && !*reason) {
213: if (fnorm < snes->ttol) {
214: PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)snes->ttol));
215: *reason = SNES_CONVERGED_FNORM_RELATIVE;
216: }
217: }
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*
222: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
224: Input Parameters:
225: . SNES - nonlinear solver context
227: Output Parameters:
228: . X - Bound projected X
230: */
232: PetscErrorCode SNESVIProjectOntoBounds(SNES snes, Vec X)
233: {
234: const PetscScalar *xl, *xu;
235: PetscScalar *x;
236: PetscInt i, n;
238: PetscFunctionBegin;
239: PetscCall(VecGetLocalSize(X, &n));
240: PetscCall(VecGetArray(X, &x));
241: PetscCall(VecGetArrayRead(snes->xl, &xl));
242: PetscCall(VecGetArrayRead(snes->xu, &xu));
244: for (i = 0; i < n; i++) {
245: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
246: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
247: }
248: PetscCall(VecRestoreArray(X, &x));
249: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
250: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: /*
255: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
257: Input parameter:
258: . snes - the SNES context
259: . X - the snes solution vector
260: . F - the nonlinear function vector
262: Output parameter:
263: . ISact - active set index set
264: */
265: PetscErrorCode SNESVIGetActiveSetIS(SNES snes, Vec X, Vec F, IS *ISact)
266: {
267: Vec Xl = snes->xl, Xu = snes->xu;
268: const PetscScalar *x, *f, *xl, *xu;
269: PetscInt *idx_act, i, nlocal, nloc_isact = 0, ilow, ihigh, i1 = 0;
270: PetscReal zerotolerance = snes->vizerotolerance;
272: PetscFunctionBegin;
273: PetscCall(VecGetLocalSize(X, &nlocal));
274: PetscCall(VecGetOwnershipRange(X, &ilow, &ihigh));
275: PetscCall(VecGetArrayRead(X, &x));
276: PetscCall(VecGetArrayRead(Xl, &xl));
277: PetscCall(VecGetArrayRead(Xu, &xu));
278: PetscCall(VecGetArrayRead(F, &f));
279: /* Compute active set size */
280: for (i = 0; i < nlocal; i++) {
281: 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++;
282: }
284: PetscCall(PetscMalloc1(nloc_isact, &idx_act));
286: /* Set active set indices */
287: for (i = 0; i < nlocal; i++) {
288: 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;
289: }
291: /* Create active set IS */
292: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), nloc_isact, idx_act, PETSC_OWN_POINTER, ISact));
294: PetscCall(VecRestoreArrayRead(X, &x));
295: PetscCall(VecRestoreArrayRead(Xl, &xl));
296: PetscCall(VecRestoreArrayRead(Xu, &xu));
297: PetscCall(VecRestoreArrayRead(F, &f));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
302: {
303: PetscInt rstart, rend;
305: PetscFunctionBegin;
306: PetscCall(SNESVIGetActiveSetIS(snes, X, F, ISact));
307: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
308: PetscCall(ISComplement(*ISact, rstart, rend, ISinact));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes, Vec F, Vec X, PetscReal *fnorm)
313: {
314: const PetscScalar *x, *xl, *xu, *f;
315: PetscInt i, n;
316: PetscReal rnorm, zerotolerance = snes->vizerotolerance;
318: PetscFunctionBegin;
319: PetscCall(VecGetLocalSize(X, &n));
320: PetscCall(VecGetArrayRead(snes->xl, &xl));
321: PetscCall(VecGetArrayRead(snes->xu, &xu));
322: PetscCall(VecGetArrayRead(X, &x));
323: PetscCall(VecGetArrayRead(F, &f));
324: rnorm = 0.0;
325: for (i = 0; i < n; i++) {
326: 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]);
327: }
328: PetscCall(VecRestoreArrayRead(F, &f));
329: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
330: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
331: PetscCall(VecRestoreArrayRead(X, &x));
332: PetscCall(MPIU_Allreduce(&rnorm, fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
333: *fnorm = PetscSqrtReal(*fnorm);
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes, Vec xl, Vec xu)
338: {
339: PetscFunctionBegin;
340: PetscCall(DMComputeVariableBounds(snes->dm, xl, xu));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /*
345: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
346: of the SNESVI nonlinear solver.
348: Input Parameter:
349: . snes - the SNES context
351: Application Interface Routine: SNESSetUp()
353: Notes:
354: For basic use of the SNES solvers, the user need not explicitly call
355: SNESSetUp(), since these actions will automatically occur during
356: the call to SNESSolve().
357: */
358: PetscErrorCode SNESSetUp_VI(SNES snes)
359: {
360: PetscInt i_start[3], i_end[3];
362: PetscFunctionBegin;
363: PetscCall(SNESSetWorkVecs(snes, 1));
364: PetscCall(SNESSetUpMatrices(snes));
366: if (!snes->ops->computevariablebounds && snes->dm) {
367: PetscBool flag;
368: PetscCall(DMHasVariableBounds(snes->dm, &flag));
369: if (flag) snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
370: }
371: if (!snes->usersetbounds) {
372: if (snes->ops->computevariablebounds) {
373: if (!snes->xl) PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
374: if (!snes->xu) PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
375: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
376: } else if (!snes->xl && !snes->xu) {
377: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
378: PetscCall(VecDuplicate(snes->vec_sol, &snes->xl));
379: PetscCall(VecSet(snes->xl, PETSC_NINFINITY));
380: PetscCall(VecDuplicate(snes->vec_sol, &snes->xu));
381: PetscCall(VecSet(snes->xu, PETSC_INFINITY));
382: } else {
383: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
384: PetscCall(VecGetOwnershipRange(snes->vec_sol, i_start, i_end));
385: PetscCall(VecGetOwnershipRange(snes->xl, i_start + 1, i_end + 1));
386: PetscCall(VecGetOwnershipRange(snes->xu, i_start + 2, i_end + 2));
387: 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]))
388: 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.");
389: }
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
393: PetscErrorCode SNESReset_VI(SNES snes)
394: {
395: PetscFunctionBegin;
396: PetscCall(VecDestroy(&snes->xl));
397: PetscCall(VecDestroy(&snes->xu));
398: snes->usersetbounds = PETSC_FALSE;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*
403: SNESDestroy_VI - Destroys the private SNES_VI context that was created
404: with SNESCreate_VI().
406: Input Parameter:
407: . snes - the SNES context
409: Application Interface Routine: SNESDestroy()
410: */
411: PetscErrorCode SNESDestroy_VI(SNES snes)
412: {
413: PetscFunctionBegin;
414: PetscCall(PetscFree(snes->data));
416: /* clear composed functions */
417: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", NULL));
418: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", NULL));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. This allows solving
424: (differential) variable inequalities.
426: Input Parameters:
427: + snes - the `SNES` context.
428: . xl - lower bound.
429: - xu - upper bound.
431: Level: advanced
433: Notes:
434: If this routine is not called then the lower and upper bounds are set to
435: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
437: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.
439: For particular components that have no bounds you can use `PETSC_NINFINITY` or `PETSC_INFINITY`
441: `SNESVISetComputeVariableBounds()` can be used to provide a function that computes the bounds. This should be used if you are using, for example, grid
442: sequencing and need bounds set for a variety of vectors
444: .seealso: [](sec_vi), `SNES`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, 'SNESSetType()`
445: @*/
446: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
447: {
448: PetscErrorCode (*f)(SNES, Vec, Vec);
450: PetscFunctionBegin;
454: PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetVariableBounds_C", &f));
455: if (f) PetscUseMethod(snes, "SNESVISetVariableBounds_C", (SNES, Vec, Vec), (snes, xl, xu));
456: else PetscCall(SNESVISetVariableBounds_VI(snes, xl, xu));
457: snes->usersetbounds = PETSC_TRUE;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes, Vec xl, Vec xu)
462: {
463: const PetscScalar *xxl, *xxu;
464: PetscInt i, n, cnt = 0;
466: PetscFunctionBegin;
467: PetscCall(SNESGetFunction(snes, &snes->vec_func, NULL, NULL));
468: PetscCheck(snes->vec_func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
469: {
470: PetscInt xlN, xuN, N;
471: PetscCall(VecGetSize(xl, &xlN));
472: PetscCall(VecGetSize(xu, &xuN));
473: PetscCall(VecGetSize(snes->vec_func, &N));
474: PetscCheck(xlN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths lower bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xlN, N);
475: PetscCheck(xuN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths: upper bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xuN, N);
476: }
477: PetscCall(PetscObjectReference((PetscObject)xl));
478: PetscCall(PetscObjectReference((PetscObject)xu));
479: PetscCall(VecDestroy(&snes->xl));
480: PetscCall(VecDestroy(&snes->xu));
481: snes->xl = xl;
482: snes->xu = xu;
483: PetscCall(VecGetLocalSize(xl, &n));
484: PetscCall(VecGetArrayRead(xl, &xxl));
485: PetscCall(VecGetArrayRead(xu, &xxu));
486: for (i = 0; i < n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
488: PetscCall(MPIU_Allreduce(&cnt, &snes->ntruebounds, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
489: PetscCall(VecRestoreArrayRead(xl, &xxl));
490: PetscCall(VecRestoreArrayRead(xu, &xxu));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: SNESVIGetVariableBounds - Gets the lower and upper bounds for the solution vector. xl <= x <= xu. This allows solving
496: (differential) variable inequalities.
498: Input Parameters:
499: + snes - the `SNES` context.
500: . xl - lower bound (may be `NULL`)
501: - xu - upper bound (may be `NULL`)
503: Level: advanced
505: Notes:
506: These vectors are owned by the `SNESVI` and should not be destroyed by the caller
508: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, SNESVINEWTONRSLS, SNESVINEWTONSSLS, 'SNESSetType()`
509: @*/
510: PetscErrorCode SNESVIGetVariableBounds(SNES snes, Vec *xl, Vec *xu)
511: {
512: PetscFunctionBegin;
513: PetscCheck(snes->usersetbounds, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must set SNESVI bounds before calling SNESVIGetVariableBounds()");
514: if (xl) *xl = snes->xl;
515: if (xu) *xu = snes->xu;
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: PetscErrorCode SNESSetFromOptions_VI(SNES snes, PetscOptionItems *PetscOptionsObject)
520: {
521: PetscBool flg = PETSC_FALSE;
523: PetscFunctionBegin;
524: PetscOptionsHeadBegin(PetscOptionsObject, "SNES VI options");
525: PetscCall(PetscOptionsReal("-snes_vi_zero_tolerance", "Tolerance for considering x[] value to be on a bound", "None", snes->vizerotolerance, &snes->vizerotolerance, NULL));
526: PetscCall(PetscOptionsBool("-snes_vi_monitor", "Monitor all non-active variables", "SNESMonitorResidual", flg, &flg, NULL));
527: if (flg) PetscCall(SNESMonitorSet(snes, SNESMonitorVI, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)), NULL));
528: flg = PETSC_FALSE;
529: PetscCall(PetscOptionsBool("-snes_vi_monitor_residual", "Monitor residual all non-active variables; using zero for active constraints", "SNESMonitorVIResidual", flg, &flg, NULL));
530: if (flg) PetscCall(SNESMonitorSet(snes, SNESVIMonitorResidual, PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)), NULL));
531: PetscOptionsHeadEnd();
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }