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 function:
13: PetscErrorCode compute(SNES snes,Vec lower,Vec higher, void *ctx)
15: + snes - the `SNES` context
16: . lower - vector to hold lower bounds
17: - higher - vector to hold upper bounds
19: Level: advanced
21: Notes:
22: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.
24: For entries with no bounds you can set `PETSC_NINFINITY` or `PETSC_INFINITY`
26: You may use `SNESVISetVariableBounds()` to provide the bounds once if they will never change
28: If you have associated a `DM` with the `SNES` and provided a function to the `DM` via `DMSetVariableBounds()` that will be used automatically
29: to provide the bounds and you need not use this function.
31: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `DMSetVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`,
32: 'SNESSetType()`
33: @*/
34: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES, Vec, Vec))
35: {
36: PetscErrorCode (*f)(SNES, PetscErrorCode(*)(SNES, Vec, Vec));
39: PetscObjectQueryFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", &f);
40: if (f) PetscUseMethod(snes, "SNESVISetComputeVariableBounds_C", (SNES, PetscErrorCode(*)(SNES, Vec, Vec)), (snes, compute));
41: else SNESVISetComputeVariableBounds_VI(snes, compute);
42: return 0;
43: }
45: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes, SNESVIComputeVariableBoundsFunction compute)
46: {
47: snes->ops->computevariablebounds = compute;
48: return 0;
49: }
51: 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: SNESGetFunction(snes, &F, NULL, NULL);
58: SNESGetSolution(snes, &X);
59: SNESVIGetActiveSetIS(snes, X, F, &isactive);
60: VecDuplicate(F, &Finactive);
61: VecCopy(F, Finactive);
62: VecISSet(Finactive, isactive, 0.0);
63: ISDestroy(&isactive);
64: VecView(Finactive, viewer);
65: VecDestroy(&Finactive);
66: return 0;
67: }
69: PetscErrorCode SNESMonitorVI(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
70: {
71: PetscViewer viewer = (PetscViewer)dummy;
72: const PetscScalar *x, *xl, *xu, *f;
73: PetscInt i, n, act[2] = {0, 0}, fact[2], N;
74: /* Number of components that actually hit the bounds (c.f. active variables) */
75: PetscInt act_bound[2] = {0, 0}, fact_bound[2];
76: PetscReal rnorm, fnorm, zerotolerance = snes->vizerotolerance;
77: double tmp;
80: VecGetLocalSize(snes->vec_sol, &n);
81: VecGetSize(snes->vec_sol, &N);
82: VecGetArrayRead(snes->xl, &xl);
83: VecGetArrayRead(snes->xu, &xu);
84: VecGetArrayRead(snes->vec_sol, &x);
85: VecGetArrayRead(snes->vec_func, &f);
87: rnorm = 0.0;
88: for (i = 0; i < n; i++) {
89: 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]);
90: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
91: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
92: else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Can never get here");
93: }
95: for (i = 0; i < n; i++) {
96: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
97: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
98: }
99: VecRestoreArrayRead(snes->vec_func, &f);
100: VecRestoreArrayRead(snes->xl, &xl);
101: VecRestoreArrayRead(snes->xu, &xu);
102: VecRestoreArrayRead(snes->vec_sol, &x);
103: MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes));
104: MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes));
105: MPIU_Allreduce(act_bound, fact_bound, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes));
106: fnorm = PetscSqrtReal(fnorm);
108: PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
109: if (snes->ntruebounds) tmp = ((double)(fact[0] + fact[1])) / ((double)snes->ntruebounds);
110: else tmp = 0.0;
111: 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);
113: PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
114: return 0;
115: }
117: /*
118: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
119: || 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
120: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
121: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
122: */
123: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes, Mat A, Vec F, Vec W, PetscReal fnorm, PetscBool *ismin)
124: {
125: PetscReal a1;
126: PetscBool hastranspose;
128: *ismin = PETSC_FALSE;
129: MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose);
130: if (hastranspose) {
131: /* Compute || J^T F|| */
132: MatMultTranspose(A, F, W);
133: VecNorm(W, NORM_2, &a1);
134: PetscInfo(snes, "|| J^T F|| %g near zero implies found a local minimum\n", (double)(a1 / fnorm));
135: if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
136: } else {
137: Vec work;
138: PetscScalar result;
139: PetscReal wnorm;
141: VecSetRandom(W, NULL);
142: VecNorm(W, NORM_2, &wnorm);
143: VecDuplicate(W, &work);
144: MatMult(A, W, work);
145: VecDot(F, work, &result);
146: VecDestroy(&work);
147: a1 = PetscAbsScalar(result) / (fnorm * wnorm);
148: PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n", (double)a1);
149: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
150: }
151: return 0;
152: }
154: /*
155: Checks if J^T(F - J*X) = 0
156: */
157: PetscErrorCode SNESVICheckResidual_Private(SNES snes, Mat A, Vec F, Vec X, Vec W1, Vec W2)
158: {
159: PetscReal a1, a2;
160: PetscBool hastranspose;
162: MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose);
163: if (hastranspose) {
164: MatMult(A, X, W1);
165: VecAXPY(W1, -1.0, F);
167: /* Compute || J^T W|| */
168: MatMultTranspose(A, W1, W2);
169: VecNorm(W1, NORM_2, &a1);
170: VecNorm(W2, NORM_2, &a2);
171: if (a1 != 0.0) PetscInfo(snes, "||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n", (double)(a2 / a1));
172: }
173: return 0;
174: }
176: /*
177: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
179: Notes:
180: The convergence criterion currently implemented is
181: merit < abstol
182: merit < rtol*merit_initial
183: */
184: PetscErrorCode SNESConvergedDefault_VI(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gradnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
185: {
189: *reason = SNES_CONVERGED_ITERATING;
191: if (!it) {
192: /* set parameter for default relative tolerance convergence test */
193: snes->ttol = fnorm * snes->rtol;
194: }
195: if (fnorm != fnorm) {
196: PetscInfo(snes, "Failed to converged, function norm is NaN\n");
197: *reason = SNES_DIVERGED_FNORM_NAN;
198: } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
199: PetscInfo(snes, "Converged due to function norm %g < %g\n", (double)fnorm, (double)snes->abstol);
200: *reason = SNES_CONVERGED_FNORM_ABS;
201: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
202: PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", snes->nfuncs, snes->max_funcs);
203: *reason = SNES_DIVERGED_FUNCTION_COUNT;
204: }
206: if (it && !*reason) {
207: if (fnorm < snes->ttol) {
208: PetscInfo(snes, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)snes->ttol);
209: *reason = SNES_CONVERGED_FNORM_RELATIVE;
210: }
211: }
212: return 0;
213: }
215: /*
216: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
218: Input Parameters:
219: . SNES - nonlinear solver context
221: Output Parameters:
222: . X - Bound projected X
224: */
226: PetscErrorCode SNESVIProjectOntoBounds(SNES snes, Vec X)
227: {
228: const PetscScalar *xl, *xu;
229: PetscScalar *x;
230: PetscInt i, n;
232: VecGetLocalSize(X, &n);
233: VecGetArray(X, &x);
234: VecGetArrayRead(snes->xl, &xl);
235: VecGetArrayRead(snes->xu, &xu);
237: for (i = 0; i < n; i++) {
238: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
239: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
240: }
241: VecRestoreArray(X, &x);
242: VecRestoreArrayRead(snes->xl, &xl);
243: VecRestoreArrayRead(snes->xu, &xu);
244: return 0;
245: }
247: /*
248: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
250: Input parameter:
251: . snes - the SNES context
252: . X - the snes solution vector
253: . F - the nonlinear function vector
255: Output parameter:
256: . ISact - active set index set
257: */
258: PetscErrorCode SNESVIGetActiveSetIS(SNES snes, Vec X, Vec F, IS *ISact)
259: {
260: Vec Xl = snes->xl, Xu = snes->xu;
261: const PetscScalar *x, *f, *xl, *xu;
262: PetscInt *idx_act, i, nlocal, nloc_isact = 0, ilow, ihigh, i1 = 0;
263: PetscReal zerotolerance = snes->vizerotolerance;
265: VecGetLocalSize(X, &nlocal);
266: VecGetOwnershipRange(X, &ilow, &ihigh);
267: VecGetArrayRead(X, &x);
268: VecGetArrayRead(Xl, &xl);
269: VecGetArrayRead(Xu, &xu);
270: VecGetArrayRead(F, &f);
271: /* Compute active set size */
272: for (i = 0; i < nlocal; i++) {
273: 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++;
274: }
276: PetscMalloc1(nloc_isact, &idx_act);
278: /* Set active set indices */
279: for (i = 0; i < nlocal; i++) {
280: 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;
281: }
283: /* Create active set IS */
284: ISCreateGeneral(PetscObjectComm((PetscObject)snes), nloc_isact, idx_act, PETSC_OWN_POINTER, ISact);
286: VecRestoreArrayRead(X, &x);
287: VecRestoreArrayRead(Xl, &xl);
288: VecRestoreArrayRead(Xu, &xu);
289: VecRestoreArrayRead(F, &f);
290: return 0;
291: }
293: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
294: {
295: PetscInt rstart, rend;
297: SNESVIGetActiveSetIS(snes, X, F, ISact);
298: VecGetOwnershipRange(X, &rstart, &rend);
299: ISComplement(*ISact, rstart, rend, ISinact);
300: return 0;
301: }
303: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes, Vec F, Vec X, PetscReal *fnorm)
304: {
305: const PetscScalar *x, *xl, *xu, *f;
306: PetscInt i, n;
307: PetscReal rnorm, zerotolerance = snes->vizerotolerance;
309: VecGetLocalSize(X, &n);
310: VecGetArrayRead(snes->xl, &xl);
311: VecGetArrayRead(snes->xu, &xu);
312: VecGetArrayRead(X, &x);
313: VecGetArrayRead(F, &f);
314: rnorm = 0.0;
315: for (i = 0; i < n; i++) {
316: 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]);
317: }
318: VecRestoreArrayRead(F, &f);
319: VecRestoreArrayRead(snes->xl, &xl);
320: VecRestoreArrayRead(snes->xu, &xu);
321: VecRestoreArrayRead(X, &x);
322: MPIU_Allreduce(&rnorm, fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes));
323: *fnorm = PetscSqrtReal(*fnorm);
324: return 0;
325: }
327: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes, Vec xl, Vec xu)
328: {
329: DMComputeVariableBounds(snes->dm, xl, xu);
330: return 0;
331: }
333: /*
334: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
335: of the SNESVI nonlinear solver.
337: Input Parameter:
338: . snes - the SNES context
340: Application Interface Routine: SNESSetUp()
342: Notes:
343: For basic use of the SNES solvers, the user need not explicitly call
344: SNESSetUp(), since these actions will automatically occur during
345: the call to SNESSolve().
346: */
347: PetscErrorCode SNESSetUp_VI(SNES snes)
348: {
349: PetscInt i_start[3], i_end[3];
351: SNESSetWorkVecs(snes, 1);
352: SNESSetUpMatrices(snes);
354: if (!snes->ops->computevariablebounds && snes->dm) {
355: PetscBool flag;
356: DMHasVariableBounds(snes->dm, &flag);
357: if (flag) snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
358: }
359: if (!snes->usersetbounds) {
360: if (snes->ops->computevariablebounds) {
361: if (!snes->xl) VecDuplicate(snes->vec_sol, &snes->xl);
362: if (!snes->xu) VecDuplicate(snes->vec_sol, &snes->xu);
363: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
364: } else if (!snes->xl && !snes->xu) {
365: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
366: VecDuplicate(snes->vec_sol, &snes->xl);
367: VecSet(snes->xl, PETSC_NINFINITY);
368: VecDuplicate(snes->vec_sol, &snes->xu);
369: VecSet(snes->xu, PETSC_INFINITY);
370: } else {
371: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
372: VecGetOwnershipRange(snes->vec_sol, i_start, i_end);
373: VecGetOwnershipRange(snes->xl, i_start + 1, i_end + 1);
374: VecGetOwnershipRange(snes->xu, i_start + 2, i_end + 2);
375: 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]))
376: 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.");
377: }
378: }
379: return 0;
380: }
381: PetscErrorCode SNESReset_VI(SNES snes)
382: {
383: VecDestroy(&snes->xl);
384: VecDestroy(&snes->xu);
385: snes->usersetbounds = PETSC_FALSE;
386: return 0;
387: }
389: /*
390: SNESDestroy_VI - Destroys the private SNES_VI context that was created
391: with SNESCreate_VI().
393: Input Parameter:
394: . snes - the SNES context
396: Application Interface Routine: SNESDestroy()
397: */
398: PetscErrorCode SNESDestroy_VI(SNES snes)
399: {
400: PetscFree(snes->data);
402: /* clear composed functions */
403: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", NULL);
404: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", NULL);
405: return 0;
406: }
408: /*@
409: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. This allows solving
410: (differential) variable inequalities.
412: Input Parameters:
413: + snes - the `SNES` context.
414: . xl - lower bound.
415: - xu - upper bound.
417: Notes:
418: If this routine is not called then the lower and upper bounds are set to
419: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
421: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.
423: For particular components that have no bounds you can use `PETSC_NINFINITY` or `PETSC_INFINITY`
425: `SNESVISetComputeVariableBounds()` can be used to provide a function that computes the bounds. This should be used if you are using, for example, grid
426: sequencing and need bounds set for a variety of vectors
428: Level: advanced
430: .seealso: [](sec_vi), `SNES`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, SNESVINEWTONRSLS, SNESVINEWTONSSLS, 'SNESSetType()`
431: @*/
432: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
433: {
434: PetscErrorCode (*f)(SNES, Vec, Vec);
439: PetscObjectQueryFunction((PetscObject)snes, "SNESVISetVariableBounds_C", &f);
440: if (f) PetscUseMethod(snes, "SNESVISetVariableBounds_C", (SNES, Vec, Vec), (snes, xl, xu));
441: else SNESVISetVariableBounds_VI(snes, xl, xu);
442: snes->usersetbounds = PETSC_TRUE;
443: return 0;
444: }
446: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes, Vec xl, Vec xu)
447: {
448: const PetscScalar *xxl, *xxu;
449: PetscInt i, n, cnt = 0;
451: SNESGetFunction(snes, &snes->vec_func, NULL, NULL);
453: {
454: PetscInt xlN, xuN, N;
455: VecGetSize(xl, &xlN);
456: VecGetSize(xu, &xuN);
457: VecGetSize(snes->vec_func, &N);
460: }
461: PetscObjectReference((PetscObject)xl);
462: PetscObjectReference((PetscObject)xu);
463: VecDestroy(&snes->xl);
464: VecDestroy(&snes->xu);
465: snes->xl = xl;
466: snes->xu = xu;
467: VecGetLocalSize(xl, &n);
468: VecGetArrayRead(xl, &xxl);
469: VecGetArrayRead(xu, &xxu);
470: for (i = 0; i < n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
472: MPIU_Allreduce(&cnt, &snes->ntruebounds, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes));
473: VecRestoreArrayRead(xl, &xxl);
474: VecRestoreArrayRead(xu, &xxu);
475: return 0;
476: }
478: PetscErrorCode SNESSetFromOptions_VI(SNES snes, PetscOptionItems *PetscOptionsObject)
479: {
480: PetscBool flg = PETSC_FALSE;
482: PetscOptionsHeadBegin(PetscOptionsObject, "SNES VI options");
483: PetscOptionsReal("-snes_vi_zero_tolerance", "Tolerance for considering x[] value to be on a bound", "None", snes->vizerotolerance, &snes->vizerotolerance, NULL);
484: PetscOptionsBool("-snes_vi_monitor", "Monitor all non-active variables", "SNESMonitorResidual", flg, &flg, NULL);
485: if (flg) SNESMonitorSet(snes, SNESMonitorVI, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)), NULL);
486: flg = PETSC_FALSE;
487: PetscOptionsBool("-snes_vi_monitor_residual", "Monitor residual all non-active variables; using zero for active constraints", "SNESMonitorVIResidual", flg, &flg, NULL);
488: if (flg) SNESMonitorSet(snes, SNESVIMonitorResidual, PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)), NULL);
489: PetscOptionsHeadEnd();
490: return 0;
491: }