Actual source code: vi.c
petsc-3.6.1 2015-08-06
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: if (!viewer) {
62: viewer = PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes));
63: }
64: VecView(Finactive,viewer);
65: VecDestroy(&Finactive);
66: return(0);
67: }
71: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
72: {
73: PetscErrorCode ierr;
74: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
75: const PetscScalar *x,*xl,*xu,*f;
76: PetscInt i,n,act[2] = {0,0},fact[2],N;
77: /* Number of components that actually hit the bounds (c.f. active variables) */
78: PetscInt act_bound[2] = {0,0},fact_bound[2];
79: PetscReal rnorm,fnorm,zerotolerance = snes->vizerotolerance;
80: double tmp;
83: VecGetLocalSize(snes->vec_sol,&n);
84: VecGetSize(snes->vec_sol,&N);
85: VecGetArrayRead(snes->xl,&xl);
86: VecGetArrayRead(snes->xu,&xu);
87: VecGetArrayRead(snes->vec_sol,&x);
88: 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: VecRestoreArrayRead(snes->vec_func,&f);
103: VecRestoreArrayRead(snes->xl,&xl);
104: VecRestoreArrayRead(snes->xu,&xu);
105: VecRestoreArrayRead(snes->vec_sol,&x);
106: MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
107: MPI_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
108: MPI_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
109: fnorm = PetscSqrtReal(fnorm);
111: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
112: if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
113: else tmp = 0.0;
114: PetscViewerASCIIPrintf(viewer,"%3D SNES VI Function norm %14.12e 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);
116: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
117: return(0);
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: */
128: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
129: {
130: PetscReal a1;
132: PetscBool hastranspose;
135: *ismin = PETSC_FALSE;
136: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
137: if (hastranspose) {
138: /* Compute || J^T F|| */
139: MatMultTranspose(A,F,W);
140: VecNorm(W,NORM_2,&a1);
141: PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
142: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
143: } else {
144: Vec work;
145: PetscScalar result;
146: PetscReal wnorm;
148: VecSetRandom(W,NULL);
149: VecNorm(W,NORM_2,&wnorm);
150: VecDuplicate(W,&work);
151: MatMult(A,W,work);
152: VecDot(F,work,&result);
153: VecDestroy(&work);
154: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
155: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
156: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
157: }
158: return(0);
159: }
161: /*
162: Checks if J^T(F - J*X) = 0
163: */
166: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
167: {
168: PetscReal a1,a2;
170: PetscBool hastranspose;
173: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
174: if (hastranspose) {
175: MatMult(A,X,W1);
176: VecAXPY(W1,-1.0,F);
178: /* Compute || J^T W|| */
179: MatMultTranspose(A,W1,W2);
180: VecNorm(W1,NORM_2,&a1);
181: VecNorm(W2,NORM_2,&a2);
182: if (a1 != 0.0) {
183: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
184: }
185: }
186: return(0);
187: }
189: /*
190: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
192: Notes:
193: The convergence criterion currently implemented is
194: merit < abstol
195: merit < rtol*merit_initial
196: */
199: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
200: {
207: *reason = SNES_CONVERGED_ITERATING;
209: if (!it) {
210: /* set parameter for default relative tolerance convergence test */
211: snes->ttol = fnorm*snes->rtol;
212: }
213: if (fnorm != fnorm) {
214: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
215: *reason = SNES_DIVERGED_FNORM_NAN;
216: } else if (fnorm < snes->abstol) {
217: PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
218: *reason = SNES_CONVERGED_FNORM_ABS;
219: } else if (snes->nfuncs >= snes->max_funcs) {
220: PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
221: *reason = SNES_DIVERGED_FUNCTION_COUNT;
222: }
224: if (it && !*reason) {
225: if (fnorm < snes->ttol) {
226: PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
227: *reason = SNES_CONVERGED_FNORM_RELATIVE;
228: }
229: }
230: return(0);
231: }
234: /* -------------------------------------------------------------------------- */
235: /*
236: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
238: Input Parameters:
239: . SNES - nonlinear solver context
241: Output Parameters:
242: . X - Bound projected X
244: */
248: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
249: {
250: PetscErrorCode ierr;
251: const PetscScalar *xl,*xu;
252: PetscScalar *x;
253: PetscInt i,n;
256: VecGetLocalSize(X,&n);
257: VecGetArray(X,&x);
258: VecGetArrayRead(snes->xl,&xl);
259: VecGetArrayRead(snes->xu,&xu);
261: for (i = 0; i<n; i++) {
262: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
263: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
264: }
265: VecRestoreArray(X,&x);
266: VecRestoreArrayRead(snes->xl,&xl);
267: VecRestoreArrayRead(snes->xu,&xu);
268: return(0);
269: }
274: /*
275: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
277: Input parameter
278: . snes - the SNES context
279: . X - the snes solution vector
280: . F - the nonlinear function vector
282: Output parameter
283: . ISact - active set index set
284: */
285: PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact)
286: {
287: PetscErrorCode ierr;
288: Vec Xl=snes->xl,Xu=snes->xu;
289: const PetscScalar *x,*f,*xl,*xu;
290: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
291: PetscReal zerotolerance = snes->vizerotolerance;
292:
294: VecGetLocalSize(X,&nlocal);
295: VecGetOwnershipRange(X,&ilow,&ihigh);
296: VecGetArrayRead(X,&x);
297: VecGetArrayRead(Xl,&xl);
298: VecGetArrayRead(Xu,&xu);
299: VecGetArrayRead(F,&f);
300: /* Compute active set size */
301: for (i=0; i < nlocal;i++) {
302: 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++;
303: }
305: PetscMalloc1(nloc_isact,&idx_act);
307: /* Set active set indices */
308: for (i=0; i < nlocal; i++) {
309: 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;
310: }
312: /* Create active set IS */
313: ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
315: VecRestoreArrayRead(X,&x);
316: VecRestoreArrayRead(Xl,&xl);
317: VecRestoreArrayRead(Xu,&xu);
318: VecRestoreArrayRead(F,&f);
319: return(0);
320: }
324: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
325: {
327: PetscInt rstart,rend;
330: SNESVIGetActiveSetIS(snes,X,F,ISact);
331: VecGetOwnershipRange(X,&rstart,&rend);
332: ISComplement(*ISact,rstart,rend,ISinact);
333: return(0);
334: }
338: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
339: {
340: PetscErrorCode ierr;
341: const PetscScalar *x,*xl,*xu,*f;
342: PetscInt i,n;
343: PetscReal rnorm,zerotolerance = snes->vizerotolerance;;
346: VecGetLocalSize(X,&n);
347: VecGetArrayRead(snes->xl,&xl);
348: VecGetArrayRead(snes->xu,&xu);
349: VecGetArrayRead(X,&x);
350: VecGetArrayRead(F,&f);
351: rnorm = 0.0;
352: for (i=0; i<n; i++) {
353: 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]);
354: }
355: VecRestoreArrayRead(F,&f);
356: VecRestoreArrayRead(snes->xl,&xl);
357: VecRestoreArrayRead(snes->xu,&xu);
358: VecRestoreArrayRead(X,&x);
359: MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
360: *fnorm = PetscSqrtReal(*fnorm);
361: return(0);
362: }
366: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
367: {
371: DMComputeVariableBounds(snes->dm, xl, xu);
372: return(0);
373: }
376: /* -------------------------------------------------------------------------- */
377: /*
378: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
379: of the SNESVI nonlinear solver.
381: Input Parameter:
382: . snes - the SNES context
384: Application Interface Routine: SNESSetUp()
386: Notes:
387: For basic use of the SNES solvers, the user need not explicitly call
388: SNESSetUp(), since these actions will automatically occur during
389: the call to SNESSolve().
390: */
393: PetscErrorCode SNESSetUp_VI(SNES snes)
394: {
396: PetscInt i_start[3],i_end[3];
399: SNESSetWorkVecs(snes,1);
400: SNESSetUpMatrices(snes);
402: if (!snes->ops->computevariablebounds && snes->dm) {
403: PetscBool flag;
404: DMHasVariableBounds(snes->dm, &flag);
406: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
407: }
408: if (!snes->usersetbounds) {
409: if (snes->ops->computevariablebounds) {
410: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
411: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
412: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
413: } else if (!snes->xl && !snes->xu) {
414: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
415: VecDuplicate(snes->vec_sol, &snes->xl);
416: VecSet(snes->xl,PETSC_NINFINITY);
417: VecDuplicate(snes->vec_sol, &snes->xu);
418: VecSet(snes->xu,PETSC_INFINITY);
419: } else {
420: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
421: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
422: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
423: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
424: 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]))
425: 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.");
426: }
427: }
428: return(0);
429: }
430: /* -------------------------------------------------------------------------- */
433: PetscErrorCode SNESReset_VI(SNES snes)
434: {
438: VecDestroy(&snes->xl);
439: VecDestroy(&snes->xu);
440: snes->usersetbounds = PETSC_FALSE;
441: return(0);
442: }
444: /*
445: SNESDestroy_VI - Destroys the private SNES_VI context that was created
446: with SNESCreate_VI().
448: Input Parameter:
449: . snes - the SNES context
451: Application Interface Routine: SNESDestroy()
452: */
455: PetscErrorCode SNESDestroy_VI(SNES snes)
456: {
460: PetscFree(snes->data);
462: /* clear composed functions */
463: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
464: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetMonitor_C",NULL);
465: return(0);
466: }
470: /*@
471: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
473: Input Parameters:
474: . snes - the SNES context.
475: . xl - lower bound.
476: . xu - upper bound.
478: Notes:
479: If this routine is not called then the lower and upper bounds are set to
480: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
482: Level: advanced
484: @*/
485: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
486: {
487: PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
493: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
494: if (!f) {
495: SNESVISetVariableBounds_VI(snes, xl, xu);
496: } else {
497: PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
498: }
499: snes->usersetbounds = PETSC_TRUE;
500: return(0);
501: }
505: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
506: {
507: PetscErrorCode ierr;
508: const PetscScalar *xxl,*xxu;
509: PetscInt i,n, cnt = 0;
512: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
513: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
514: {
515: PetscInt xlN,xuN,N;
516: VecGetSize(xl,&xlN);
517: VecGetSize(xu,&xuN);
518: VecGetSize(snes->vec_func,&N);
519: if (xlN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xlN,N);
520: if (xuN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xuN,N);
521: }
522: PetscObjectReference((PetscObject)xl);
523: PetscObjectReference((PetscObject)xu);
524: VecDestroy(&snes->xl);
525: VecDestroy(&snes->xu);
526: snes->xl = xl;
527: snes->xu = xu;
528: VecGetLocalSize(xl,&n);
529: VecGetArrayRead(xl,&xxl);
530: VecGetArrayRead(xu,&xxu);
531: for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
533: MPI_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
534: VecRestoreArrayRead(xl,&xxl);
535: VecRestoreArrayRead(xu,&xxu);
536: return(0);
537: }
541: PetscErrorCode SNESSetFromOptions_VI(PetscOptions *PetscOptionsObject,SNES snes)
542: {
544: PetscBool flg = PETSC_FALSE;
545: SNESLineSearch linesearch;
548: PetscOptionsHead(PetscOptionsObject,"SNES VI options");
549: PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
550: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
551: if (flg) {
552: SNESMonitorSet(snes,SNESMonitorVI,0,0);
553: }
554: flg = PETSC_FALSE;
555: PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
556: if (flg) {
557: SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
558: }
559: if (!snes->linesearch) {
560: SNESGetLineSearch(snes, &linesearch);
561: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
562: SNESLineSearchBTSetAlpha(linesearch, 0.0);
563: }
564: PetscOptionsTail();
565: return(0);
566: }