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 variable bounds
7: Input parameter:
8: + snes - the SNES context
9: - compute - computes the bounds
11: Level: advanced
13: .seealso: SNESVISetVariableBounds()
15: @*/
16: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES,Vec,Vec))
17: {
18: PetscErrorCode (*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec));
21: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",&f);
22: if (f) PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));
23: else SNESVISetComputeVariableBounds_VI(snes,compute);
24: return 0;
25: }
27: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
28: {
29: snes->ops->computevariablebounds = compute;
30: return 0;
31: }
33: /* --------------------------------------------------------------------------------------------------------*/
35: PetscErrorCode SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
36: {
37: Vec X, F, Finactive;
38: IS isactive;
39: PetscViewer viewer = (PetscViewer) dummy;
41: SNESGetFunction(snes,&F,NULL,NULL);
42: SNESGetSolution(snes,&X);
43: SNESVIGetActiveSetIS(snes,X,F,&isactive);
44: VecDuplicate(F,&Finactive);
45: VecCopy(F,Finactive);
46: VecISSet(Finactive,isactive,0.0);
47: ISDestroy(&isactive);
48: VecView(Finactive,viewer);
49: VecDestroy(&Finactive);
50: return 0;
51: }
53: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
54: {
55: PetscViewer viewer = (PetscViewer) dummy;
56: const PetscScalar *x,*xl,*xu,*f;
57: PetscInt i,n,act[2] = {0,0},fact[2],N;
58: /* Number of components that actually hit the bounds (c.f. active variables) */
59: PetscInt act_bound[2] = {0,0},fact_bound[2];
60: PetscReal rnorm,fnorm,zerotolerance = snes->vizerotolerance;
61: double tmp;
64: VecGetLocalSize(snes->vec_sol,&n);
65: VecGetSize(snes->vec_sol,&N);
66: VecGetArrayRead(snes->xl,&xl);
67: VecGetArrayRead(snes->xu,&xu);
68: VecGetArrayRead(snes->vec_sol,&x);
69: VecGetArrayRead(snes->vec_func,&f);
71: rnorm = 0.0;
72: for (i=0; i<n; i++) {
73: 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]);
74: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
75: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
76: else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
77: }
79: for (i=0; i<n; i++) {
80: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
81: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
82: }
83: VecRestoreArrayRead(snes->vec_func,&f);
84: VecRestoreArrayRead(snes->xl,&xl);
85: VecRestoreArrayRead(snes->xu,&xu);
86: VecRestoreArrayRead(snes->vec_sol,&x);
87: MPIU_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
88: MPIU_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
89: MPIU_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
90: fnorm = PetscSqrtReal(fnorm);
92: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
93: if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
94: else tmp = 0.0;
95: PetscViewerASCIIPrintf(viewer,"%3D SNES VI Function norm %g Active lower constraints %D/%D upper constraints %D/%D 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);
97: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
98: return 0;
99: }
101: /*
102: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
103: || 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
104: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
105: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
106: */
107: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
108: {
109: PetscReal a1;
110: PetscBool hastranspose;
112: *ismin = PETSC_FALSE;
113: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
114: if (hastranspose) {
115: /* Compute || J^T F|| */
116: MatMultTranspose(A,F,W);
117: VecNorm(W,NORM_2,&a1);
118: PetscInfo(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
119: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
120: } else {
121: Vec work;
122: PetscScalar result;
123: PetscReal wnorm;
125: VecSetRandom(W,NULL);
126: VecNorm(W,NORM_2,&wnorm);
127: VecDuplicate(W,&work);
128: MatMult(A,W,work);
129: VecDot(F,work,&result);
130: VecDestroy(&work);
131: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
132: PetscInfo(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
133: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
134: }
135: return 0;
136: }
138: /*
139: Checks if J^T(F - J*X) = 0
140: */
141: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
142: {
143: PetscReal a1,a2;
144: PetscBool hastranspose;
146: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
147: if (hastranspose) {
148: MatMult(A,X,W1);
149: VecAXPY(W1,-1.0,F);
151: /* Compute || J^T W|| */
152: MatMultTranspose(A,W1,W2);
153: VecNorm(W1,NORM_2,&a1);
154: VecNorm(W2,NORM_2,&a2);
155: if (a1 != 0.0) {
156: PetscInfo(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
157: }
158: }
159: return 0;
160: }
162: /*
163: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
165: Notes:
166: The convergence criterion currently implemented is
167: merit < abstol
168: merit < rtol*merit_initial
169: */
170: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
171: {
175: *reason = SNES_CONVERGED_ITERATING;
177: if (!it) {
178: /* set parameter for default relative tolerance convergence test */
179: snes->ttol = fnorm*snes->rtol;
180: }
181: if (fnorm != fnorm) {
182: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
183: *reason = SNES_DIVERGED_FNORM_NAN;
184: } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
185: PetscInfo(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
186: *reason = SNES_CONVERGED_FNORM_ABS;
187: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
188: PetscInfo(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
189: *reason = SNES_DIVERGED_FUNCTION_COUNT;
190: }
192: if (it && !*reason) {
193: if (fnorm < snes->ttol) {
194: PetscInfo(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
195: *reason = SNES_CONVERGED_FNORM_RELATIVE;
196: }
197: }
198: return 0;
199: }
201: /* -------------------------------------------------------------------------- */
202: /*
203: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
205: Input Parameters:
206: . SNES - nonlinear solver context
208: Output Parameters:
209: . X - Bound projected X
211: */
213: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
214: {
215: const PetscScalar *xl,*xu;
216: PetscScalar *x;
217: PetscInt i,n;
219: VecGetLocalSize(X,&n);
220: VecGetArray(X,&x);
221: VecGetArrayRead(snes->xl,&xl);
222: VecGetArrayRead(snes->xu,&xu);
224: for (i = 0; i<n; i++) {
225: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
226: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
227: }
228: VecRestoreArray(X,&x);
229: VecRestoreArrayRead(snes->xl,&xl);
230: VecRestoreArrayRead(snes->xu,&xu);
231: return 0;
232: }
234: /*
235: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
237: Input parameter:
238: . snes - the SNES context
239: . X - the snes solution vector
240: . F - the nonlinear function vector
242: Output parameter:
243: . ISact - active set index set
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: VecGetLocalSize(X,&nlocal);
253: VecGetOwnershipRange(X,&ilow,&ihigh);
254: VecGetArrayRead(X,&x);
255: VecGetArrayRead(Xl,&xl);
256: VecGetArrayRead(Xu,&xu);
257: VecGetArrayRead(F,&f);
258: /* Compute active set size */
259: for (i=0; i < nlocal;i++) {
260: 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++;
261: }
263: PetscMalloc1(nloc_isact,&idx_act);
265: /* Set active set indices */
266: for (i=0; i < nlocal; i++) {
267: 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;
268: }
270: /* Create active set IS */
271: ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
273: VecRestoreArrayRead(X,&x);
274: VecRestoreArrayRead(Xl,&xl);
275: VecRestoreArrayRead(Xu,&xu);
276: VecRestoreArrayRead(F,&f);
277: return 0;
278: }
280: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
281: {
282: PetscInt rstart,rend;
284: SNESVIGetActiveSetIS(snes,X,F,ISact);
285: VecGetOwnershipRange(X,&rstart,&rend);
286: ISComplement(*ISact,rstart,rend,ISinact);
287: return 0;
288: }
290: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
291: {
292: const PetscScalar *x,*xl,*xu,*f;
293: PetscInt i,n;
294: PetscReal rnorm,zerotolerance = snes->vizerotolerance;
296: VecGetLocalSize(X,&n);
297: VecGetArrayRead(snes->xl,&xl);
298: VecGetArrayRead(snes->xu,&xu);
299: VecGetArrayRead(X,&x);
300: VecGetArrayRead(F,&f);
301: rnorm = 0.0;
302: for (i=0; i<n; i++) {
303: 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]);
304: }
305: VecRestoreArrayRead(F,&f);
306: VecRestoreArrayRead(snes->xl,&xl);
307: VecRestoreArrayRead(snes->xu,&xu);
308: VecRestoreArrayRead(X,&x);
309: MPIU_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
310: *fnorm = PetscSqrtReal(*fnorm);
311: return 0;
312: }
314: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
315: {
316: DMComputeVariableBounds(snes->dm, xl, xu);
317: return 0;
318: }
320: /* -------------------------------------------------------------------------- */
321: /*
322: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
323: of the SNESVI nonlinear solver.
325: Input Parameter:
326: . snes - the SNES context
328: Application Interface Routine: SNESSetUp()
330: Notes:
331: For basic use of the SNES solvers, the user need not explicitly call
332: SNESSetUp(), since these actions will automatically occur during
333: the call to SNESSolve().
334: */
335: PetscErrorCode SNESSetUp_VI(SNES snes)
336: {
337: PetscInt i_start[3],i_end[3];
339: SNESSetWorkVecs(snes,1);
340: SNESSetUpMatrices(snes);
342: if (!snes->ops->computevariablebounds && snes->dm) {
343: PetscBool flag;
344: DMHasVariableBounds(snes->dm, &flag);
345: if (flag) {
346: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
347: }
348: }
349: if (!snes->usersetbounds) {
350: if (snes->ops->computevariablebounds) {
351: if (!snes->xl) VecDuplicate(snes->vec_sol,&snes->xl);
352: if (!snes->xu) VecDuplicate(snes->vec_sol,&snes->xu);
353: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
354: } else if (!snes->xl && !snes->xu) {
355: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
356: VecDuplicate(snes->vec_sol, &snes->xl);
357: VecSet(snes->xl,PETSC_NINFINITY);
358: VecDuplicate(snes->vec_sol, &snes->xu);
359: VecSet(snes->xu,PETSC_INFINITY);
360: } else {
361: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
362: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
363: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
364: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
365: 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]))
366: 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.");
367: }
368: }
369: return 0;
370: }
371: /* -------------------------------------------------------------------------- */
372: PetscErrorCode SNESReset_VI(SNES snes)
373: {
374: VecDestroy(&snes->xl);
375: VecDestroy(&snes->xu);
376: snes->usersetbounds = PETSC_FALSE;
377: return 0;
378: }
380: /*
381: SNESDestroy_VI - Destroys the private SNES_VI context that was created
382: with SNESCreate_VI().
384: Input Parameter:
385: . snes - the SNES context
387: Application Interface Routine: SNESDestroy()
388: */
389: PetscErrorCode SNESDestroy_VI(SNES snes)
390: {
391: PetscFree(snes->data);
393: /* clear composed functions */
394: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
395: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetDefaultMonitor_C",NULL);
396: return 0;
397: }
399: /*@
400: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
402: Input Parameters:
403: + snes - the SNES context.
404: . xl - lower bound.
405: - xu - upper bound.
407: Notes:
408: If this routine is not called then the lower and upper bounds are set to
409: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
411: Level: advanced
413: @*/
414: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
415: {
416: PetscErrorCode (*f)(SNES,Vec,Vec);
421: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
422: if (f) PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
423: else SNESVISetVariableBounds_VI(snes, xl, xu);
424: snes->usersetbounds = PETSC_TRUE;
425: return 0;
426: }
428: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
429: {
430: const PetscScalar *xxl,*xxu;
431: PetscInt i,n, cnt = 0;
433: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
435: {
436: PetscInt xlN,xuN,N;
437: VecGetSize(xl,&xlN);
438: VecGetSize(xu,&xuN);
439: VecGetSize(snes->vec_func,&N);
442: }
443: PetscObjectReference((PetscObject)xl);
444: PetscObjectReference((PetscObject)xu);
445: VecDestroy(&snes->xl);
446: VecDestroy(&snes->xu);
447: snes->xl = xl;
448: snes->xu = xu;
449: VecGetLocalSize(xl,&n);
450: VecGetArrayRead(xl,&xxl);
451: VecGetArrayRead(xu,&xxu);
452: for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
454: MPIU_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
455: VecRestoreArrayRead(xl,&xxl);
456: VecRestoreArrayRead(xu,&xxu);
457: return 0;
458: }
460: PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems *PetscOptionsObject,SNES snes)
461: {
462: PetscBool flg = PETSC_FALSE;
464: PetscOptionsHead(PetscOptionsObject,"SNES VI options");
465: PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
466: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
467: if (flg) {
468: SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);
469: }
470: flg = PETSC_FALSE;
471: PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
472: if (flg) {
473: SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
474: }
475: PetscOptionsTail();
476: return 0;
477: }