Actual source code: vi.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2: #include <petscdm.h>
6: /*@C
7: SNESVISetComputeVariableBounds - Sets a function that is called to compute the variable bounds
9: Input parameter
10: + snes - the SNES context
11: - compute - computes the bounds
13: Level: advanced
15: .seealso: SNESVISetVariableBounds()
17: @*/
18: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES,Vec,Vec))
19: {
20: PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec));
24: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",&f);
25: if (!f) {
26: SNESVISetComputeVariableBounds_VI(snes,compute);
27: } else {
28: PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));
29: }
30: return(0);
31: }
35: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
36: {
38: snes->ops->computevariablebounds = compute;
39: return(0);
40: }
42: /* --------------------------------------------------------------------------------------------------------*/
46: PetscErrorCode SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
47: {
49: Vec X, F, Finactive;
50: IS isactive;
51: PetscViewer viewer = (PetscViewer) dummy;
54: SNESGetFunction(snes,&F,0,0);
55: SNESGetSolution(snes,&X);
56: SNESVIGetActiveSetIS(snes,X,F,&isactive);
57: VecDuplicate(F,&Finactive);
58: VecCopy(F,Finactive);
59: VecISSet(Finactive,isactive,0.0);
60: ISDestroy(&isactive);
61: VecView(Finactive,viewer);
62: VecDestroy(&Finactive);
63: return(0);
64: }
68: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
69: {
70: PetscErrorCode ierr;
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;
81: VecGetLocalSize(snes->vec_sol,&n);
82: VecGetSize(snes->vec_sol,&N);
83: VecGetArrayRead(snes->xl,&xl);
84: VecGetArrayRead(snes->xu,&xu);
85: VecGetArrayRead(snes->vec_sol,&x);
86: VecGetArrayRead(snes->vec_func,&f);
88: rnorm = 0.0;
89: for (i=0; i<n; i++) {
90: 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]);
91: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
92: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
93: else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
94: }
96: for (i=0; i<n; i++) {
97: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
98: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
99: }
100: VecRestoreArrayRead(snes->vec_func,&f);
101: VecRestoreArrayRead(snes->xl,&xl);
102: VecRestoreArrayRead(snes->xu,&xu);
103: VecRestoreArrayRead(snes->vec_sol,&x);
104: MPIU_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
105: MPIU_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
106: MPIU_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
107: fnorm = PetscSqrtReal(fnorm);
109: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
110: if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
111: else tmp = 0.0;
112: 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);
114: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
115: return(0);
116: }
118: /*
119: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
120: || 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
121: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
122: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
123: */
126: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
127: {
128: PetscReal a1;
130: PetscBool hastranspose;
133: *ismin = PETSC_FALSE;
134: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
135: if (hastranspose) {
136: /* Compute || J^T F|| */
137: MatMultTranspose(A,F,W);
138: VecNorm(W,NORM_2,&a1);
139: PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
140: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
141: } else {
142: Vec work;
143: PetscScalar result;
144: PetscReal wnorm;
146: VecSetRandom(W,NULL);
147: VecNorm(W,NORM_2,&wnorm);
148: VecDuplicate(W,&work);
149: MatMult(A,W,work);
150: VecDot(F,work,&result);
151: VecDestroy(&work);
152: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
153: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
154: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
155: }
156: return(0);
157: }
159: /*
160: Checks if J^T(F - J*X) = 0
161: */
164: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
165: {
166: PetscReal a1,a2;
168: PetscBool hastranspose;
171: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
172: if (hastranspose) {
173: MatMult(A,X,W1);
174: VecAXPY(W1,-1.0,F);
176: /* Compute || J^T W|| */
177: MatMultTranspose(A,W1,W2);
178: VecNorm(W1,NORM_2,&a1);
179: VecNorm(W2,NORM_2,&a2);
180: if (a1 != 0.0) {
181: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
182: }
183: }
184: return(0);
185: }
187: /*
188: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
190: Notes:
191: The convergence criterion currently implemented is
192: merit < abstol
193: merit < rtol*merit_initial
194: */
197: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
198: {
205: *reason = SNES_CONVERGED_ITERATING;
207: if (!it) {
208: /* set parameter for default relative tolerance convergence test */
209: snes->ttol = fnorm*snes->rtol;
210: }
211: if (fnorm != fnorm) {
212: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
213: *reason = SNES_DIVERGED_FNORM_NAN;
214: } else if (fnorm < snes->abstol) {
215: PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
216: *reason = SNES_CONVERGED_FNORM_ABS;
217: } else if (snes->nfuncs >= snes->max_funcs) {
218: PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
219: *reason = SNES_DIVERGED_FUNCTION_COUNT;
220: }
222: if (it && !*reason) {
223: if (fnorm < snes->ttol) {
224: PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
225: *reason = SNES_CONVERGED_FNORM_RELATIVE;
226: }
227: }
228: return(0);
229: }
232: /* -------------------------------------------------------------------------- */
233: /*
234: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
236: Input Parameters:
237: . SNES - nonlinear solver context
239: Output Parameters:
240: . X - Bound projected X
242: */
246: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
247: {
248: PetscErrorCode ierr;
249: const PetscScalar *xl,*xu;
250: PetscScalar *x;
251: PetscInt i,n;
254: VecGetLocalSize(X,&n);
255: VecGetArray(X,&x);
256: VecGetArrayRead(snes->xl,&xl);
257: VecGetArrayRead(snes->xu,&xu);
259: for (i = 0; i<n; i++) {
260: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
261: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
262: }
263: VecRestoreArray(X,&x);
264: VecRestoreArrayRead(snes->xl,&xl);
265: VecRestoreArrayRead(snes->xu,&xu);
266: return(0);
267: }
272: /*
273: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
275: Input parameter
276: . snes - the SNES context
277: . X - the snes solution vector
278: . F - the nonlinear function vector
280: Output parameter
281: . ISact - active set index set
282: */
283: PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact)
284: {
285: PetscErrorCode ierr;
286: Vec Xl=snes->xl,Xu=snes->xu;
287: const PetscScalar *x,*f,*xl,*xu;
288: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
289: PetscReal zerotolerance = snes->vizerotolerance;
290:
292: VecGetLocalSize(X,&nlocal);
293: VecGetOwnershipRange(X,&ilow,&ihigh);
294: VecGetArrayRead(X,&x);
295: VecGetArrayRead(Xl,&xl);
296: VecGetArrayRead(Xu,&xu);
297: VecGetArrayRead(F,&f);
298: /* Compute active set size */
299: for (i=0; i < nlocal;i++) {
300: 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++;
301: }
303: PetscMalloc1(nloc_isact,&idx_act);
305: /* Set active set indices */
306: for (i=0; i < nlocal; i++) {
307: 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;
308: }
310: /* Create active set IS */
311: ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
313: VecRestoreArrayRead(X,&x);
314: VecRestoreArrayRead(Xl,&xl);
315: VecRestoreArrayRead(Xu,&xu);
316: VecRestoreArrayRead(F,&f);
317: return(0);
318: }
322: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
323: {
325: PetscInt rstart,rend;
328: SNESVIGetActiveSetIS(snes,X,F,ISact);
329: VecGetOwnershipRange(X,&rstart,&rend);
330: ISComplement(*ISact,rstart,rend,ISinact);
331: return(0);
332: }
336: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
337: {
338: PetscErrorCode ierr;
339: const PetscScalar *x,*xl,*xu,*f;
340: PetscInt i,n;
341: PetscReal rnorm,zerotolerance = snes->vizerotolerance;;
344: VecGetLocalSize(X,&n);
345: VecGetArrayRead(snes->xl,&xl);
346: VecGetArrayRead(snes->xu,&xu);
347: VecGetArrayRead(X,&x);
348: VecGetArrayRead(F,&f);
349: rnorm = 0.0;
350: for (i=0; i<n; i++) {
351: 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]);
352: }
353: VecRestoreArrayRead(F,&f);
354: VecRestoreArrayRead(snes->xl,&xl);
355: VecRestoreArrayRead(snes->xu,&xu);
356: VecRestoreArrayRead(X,&x);
357: MPIU_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
358: *fnorm = PetscSqrtReal(*fnorm);
359: return(0);
360: }
364: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
365: {
369: DMComputeVariableBounds(snes->dm, xl, xu);
370: return(0);
371: }
374: /* -------------------------------------------------------------------------- */
375: /*
376: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
377: of the SNESVI nonlinear solver.
379: Input Parameter:
380: . snes - the SNES context
382: Application Interface Routine: SNESSetUp()
384: Notes:
385: For basic use of the SNES solvers, the user need not explicitly call
386: SNESSetUp(), since these actions will automatically occur during
387: the call to SNESSolve().
388: */
391: PetscErrorCode SNESSetUp_VI(SNES snes)
392: {
394: PetscInt i_start[3],i_end[3];
397: SNESSetWorkVecs(snes,1);
398: SNESSetUpMatrices(snes);
400: if (!snes->ops->computevariablebounds && snes->dm) {
401: PetscBool flag;
402: DMHasVariableBounds(snes->dm, &flag);
403: if (flag) {
404: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
405: }
406: }
407: if (!snes->usersetbounds) {
408: if (snes->ops->computevariablebounds) {
409: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
410: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
411: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
412: } else if (!snes->xl && !snes->xu) {
413: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
414: VecDuplicate(snes->vec_sol, &snes->xl);
415: VecSet(snes->xl,PETSC_NINFINITY);
416: VecDuplicate(snes->vec_sol, &snes->xu);
417: VecSet(snes->xu,PETSC_INFINITY);
418: } else {
419: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
420: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
421: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
422: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
423: 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]))
424: 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.");
425: }
426: }
427: return(0);
428: }
429: /* -------------------------------------------------------------------------- */
432: PetscErrorCode SNESReset_VI(SNES snes)
433: {
437: VecDestroy(&snes->xl);
438: VecDestroy(&snes->xu);
439: snes->usersetbounds = PETSC_FALSE;
440: return(0);
441: }
443: /*
444: SNESDestroy_VI - Destroys the private SNES_VI context that was created
445: with SNESCreate_VI().
447: Input Parameter:
448: . snes - the SNES context
450: Application Interface Routine: SNESDestroy()
451: */
454: PetscErrorCode SNESDestroy_VI(SNES snes)
455: {
459: PetscFree(snes->data);
461: /* clear composed functions */
462: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
463: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetDefaultMonitor_C",NULL);
464: return(0);
465: }
469: /*@
470: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
472: Input Parameters:
473: . snes - the SNES context.
474: . xl - lower bound.
475: . xu - upper bound.
477: Notes:
478: If this routine is not called then the lower and upper bounds are set to
479: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
481: Level: advanced
483: @*/
484: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
485: {
486: PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
492: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
493: if (!f) {
494: SNESVISetVariableBounds_VI(snes, xl, xu);
495: } else {
496: PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
497: }
498: snes->usersetbounds = PETSC_TRUE;
499: return(0);
500: }
504: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
505: {
506: PetscErrorCode ierr;
507: const PetscScalar *xxl,*xxu;
508: PetscInt i,n, cnt = 0;
511: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
512: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
513: {
514: PetscInt xlN,xuN,N;
515: VecGetSize(xl,&xlN);
516: VecGetSize(xu,&xuN);
517: VecGetSize(snes->vec_func,&N);
518: if (xlN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xlN,N);
519: if (xuN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xuN,N);
520: }
521: PetscObjectReference((PetscObject)xl);
522: PetscObjectReference((PetscObject)xu);
523: VecDestroy(&snes->xl);
524: VecDestroy(&snes->xu);
525: snes->xl = xl;
526: snes->xu = xu;
527: VecGetLocalSize(xl,&n);
528: VecGetArrayRead(xl,&xxl);
529: VecGetArrayRead(xu,&xxu);
530: for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
532: MPIU_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
533: VecRestoreArrayRead(xl,&xxl);
534: VecRestoreArrayRead(xu,&xxu);
535: return(0);
536: }
540: PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems *PetscOptionsObject,SNES snes)
541: {
543: PetscBool flg = PETSC_FALSE;
544: SNESLineSearch linesearch;
547: PetscOptionsHead(PetscOptionsObject,"SNES VI options");
548: PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
549: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
550: if (flg) {
551: SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);
552: }
553: flg = PETSC_FALSE;
554: PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
555: if (flg) {
556: SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
557: }
558: if (!snes->linesearch) {
559: SNESGetLineSearch(snes, &linesearch);
560: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
561: SNESLineSearchBTSetAlpha(linesearch, 0.0);
562: }
563: PetscOptionsTail();
564: return(0);
565: }