Actual source code: vi.c
petsc-3.10.5 2019-03-28
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 ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec));
22: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",&f);
23: if (!f) {
24: SNESVISetComputeVariableBounds_VI(snes,compute);
25: } else {
26: PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));
27: }
28: return(0);
29: }
31: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
32: {
34: snes->ops->computevariablebounds = compute;
35: return(0);
36: }
38: /* --------------------------------------------------------------------------------------------------------*/
40: PetscErrorCode SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
41: {
43: Vec X, F, Finactive;
44: IS isactive;
45: PetscViewer viewer = (PetscViewer) dummy;
48: SNESGetFunction(snes,&F,0,0);
49: SNESGetSolution(snes,&X);
50: SNESVIGetActiveSetIS(snes,X,F,&isactive);
51: VecDuplicate(F,&Finactive);
52: VecCopy(F,Finactive);
53: VecISSet(Finactive,isactive,0.0);
54: ISDestroy(&isactive);
55: VecView(Finactive,viewer);
56: VecDestroy(&Finactive);
57: return(0);
58: }
60: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
61: {
62: PetscErrorCode ierr;
63: PetscViewer viewer = (PetscViewer) dummy;
64: const PetscScalar *x,*xl,*xu,*f;
65: PetscInt i,n,act[2] = {0,0},fact[2],N;
66: /* Number of components that actually hit the bounds (c.f. active variables) */
67: PetscInt act_bound[2] = {0,0},fact_bound[2];
68: PetscReal rnorm,fnorm,zerotolerance = snes->vizerotolerance;
69: double tmp;
73: VecGetLocalSize(snes->vec_sol,&n);
74: VecGetSize(snes->vec_sol,&N);
75: VecGetArrayRead(snes->xl,&xl);
76: VecGetArrayRead(snes->xu,&xu);
77: VecGetArrayRead(snes->vec_sol,&x);
78: VecGetArrayRead(snes->vec_func,&f);
80: rnorm = 0.0;
81: for (i=0; i<n; i++) {
82: 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]);
83: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
84: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
85: else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
86: }
88: for (i=0; i<n; i++) {
89: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
90: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
91: }
92: VecRestoreArrayRead(snes->vec_func,&f);
93: VecRestoreArrayRead(snes->xl,&xl);
94: VecRestoreArrayRead(snes->xu,&xu);
95: VecRestoreArrayRead(snes->vec_sol,&x);
96: MPIU_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
97: MPIU_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
98: MPIU_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
99: fnorm = PetscSqrtReal(fnorm);
101: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
102: if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
103: else tmp = 0.0;
104: 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);
106: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
107: return(0);
108: }
110: /*
111: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
112: || 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
113: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
114: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
115: */
116: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
117: {
118: PetscReal a1;
120: PetscBool hastranspose;
123: *ismin = PETSC_FALSE;
124: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
125: if (hastranspose) {
126: /* Compute || J^T F|| */
127: MatMultTranspose(A,F,W);
128: VecNorm(W,NORM_2,&a1);
129: PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
130: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
131: } else {
132: Vec work;
133: PetscScalar result;
134: PetscReal wnorm;
136: VecSetRandom(W,NULL);
137: VecNorm(W,NORM_2,&wnorm);
138: VecDuplicate(W,&work);
139: MatMult(A,W,work);
140: VecDot(F,work,&result);
141: VecDestroy(&work);
142: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
143: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
144: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
145: }
146: return(0);
147: }
149: /*
150: Checks if J^T(F - J*X) = 0
151: */
152: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
153: {
154: PetscReal a1,a2;
156: PetscBool hastranspose;
159: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
160: if (hastranspose) {
161: MatMult(A,X,W1);
162: VecAXPY(W1,-1.0,F);
164: /* Compute || J^T W|| */
165: MatMultTranspose(A,W1,W2);
166: VecNorm(W1,NORM_2,&a1);
167: VecNorm(W2,NORM_2,&a2);
168: if (a1 != 0.0) {
169: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
170: }
171: }
172: return(0);
173: }
175: /*
176: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
178: Notes:
179: The convergence criterion currently implemented is
180: merit < abstol
181: merit < rtol*merit_initial
182: */
183: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
184: {
191: *reason = SNES_CONVERGED_ITERATING;
193: if (!it) {
194: /* set parameter for default relative tolerance convergence test */
195: snes->ttol = fnorm*snes->rtol;
196: }
197: if (fnorm != fnorm) {
198: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
199: *reason = SNES_DIVERGED_FNORM_NAN;
200: } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
201: PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
202: *reason = SNES_CONVERGED_FNORM_ABS;
203: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
204: PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
205: *reason = SNES_DIVERGED_FUNCTION_COUNT;
206: }
208: if (it && !*reason) {
209: if (fnorm < snes->ttol) {
210: PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
211: *reason = SNES_CONVERGED_FNORM_RELATIVE;
212: }
213: }
214: return(0);
215: }
218: /* -------------------------------------------------------------------------- */
219: /*
220: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
222: Input Parameters:
223: . SNES - nonlinear solver context
225: Output Parameters:
226: . X - Bound projected X
228: */
230: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
231: {
232: PetscErrorCode ierr;
233: const PetscScalar *xl,*xu;
234: PetscScalar *x;
235: PetscInt i,n;
238: VecGetLocalSize(X,&n);
239: VecGetArray(X,&x);
240: VecGetArrayRead(snes->xl,&xl);
241: VecGetArrayRead(snes->xu,&xu);
243: for (i = 0; i<n; i++) {
244: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
245: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
246: }
247: VecRestoreArray(X,&x);
248: VecRestoreArrayRead(snes->xl,&xl);
249: VecRestoreArrayRead(snes->xu,&xu);
250: return(0);
251: }
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: PetscErrorCode ierr;
268: Vec Xl=snes->xl,Xu=snes->xu;
269: const PetscScalar *x,*f,*xl,*xu;
270: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
271: PetscReal zerotolerance = snes->vizerotolerance;
272:
274: VecGetLocalSize(X,&nlocal);
275: VecGetOwnershipRange(X,&ilow,&ihigh);
276: VecGetArrayRead(X,&x);
277: VecGetArrayRead(Xl,&xl);
278: VecGetArrayRead(Xu,&xu);
279: VecGetArrayRead(F,&f);
280: /* Compute active set size */
281: for (i=0; i < nlocal;i++) {
282: 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++;
283: }
285: PetscMalloc1(nloc_isact,&idx_act);
287: /* Set active set indices */
288: for (i=0; i < nlocal; i++) {
289: 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;
290: }
292: /* Create active set IS */
293: ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
295: VecRestoreArrayRead(X,&x);
296: VecRestoreArrayRead(Xl,&xl);
297: VecRestoreArrayRead(Xu,&xu);
298: VecRestoreArrayRead(F,&f);
299: return(0);
300: }
302: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
303: {
305: PetscInt rstart,rend;
308: SNESVIGetActiveSetIS(snes,X,F,ISact);
309: VecGetOwnershipRange(X,&rstart,&rend);
310: ISComplement(*ISact,rstart,rend,ISinact);
311: return(0);
312: }
314: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
315: {
316: PetscErrorCode ierr;
317: const PetscScalar *x,*xl,*xu,*f;
318: PetscInt i,n;
319: PetscReal rnorm,zerotolerance = snes->vizerotolerance;;
322: VecGetLocalSize(X,&n);
323: VecGetArrayRead(snes->xl,&xl);
324: VecGetArrayRead(snes->xu,&xu);
325: VecGetArrayRead(X,&x);
326: VecGetArrayRead(F,&f);
327: rnorm = 0.0;
328: for (i=0; i<n; i++) {
329: 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]);
330: }
331: VecRestoreArrayRead(F,&f);
332: VecRestoreArrayRead(snes->xl,&xl);
333: VecRestoreArrayRead(snes->xu,&xu);
334: VecRestoreArrayRead(X,&x);
335: MPIU_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
336: *fnorm = PetscSqrtReal(*fnorm);
337: return(0);
338: }
340: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
341: {
345: DMComputeVariableBounds(snes->dm, xl, xu);
346: return(0);
347: }
350: /* -------------------------------------------------------------------------- */
351: /*
352: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
353: of the SNESVI nonlinear solver.
355: Input Parameter:
356: . snes - the SNES context
358: Application Interface Routine: SNESSetUp()
360: Notes:
361: For basic use of the SNES solvers, the user need not explicitly call
362: SNESSetUp(), since these actions will automatically occur during
363: the call to SNESSolve().
364: */
365: PetscErrorCode SNESSetUp_VI(SNES snes)
366: {
368: PetscInt i_start[3],i_end[3];
371: SNESSetWorkVecs(snes,1);
372: SNESSetUpMatrices(snes);
374: if (!snes->ops->computevariablebounds && snes->dm) {
375: PetscBool flag;
376: DMHasVariableBounds(snes->dm, &flag);
377: if (flag) {
378: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
379: }
380: }
381: if (!snes->usersetbounds) {
382: if (snes->ops->computevariablebounds) {
383: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
384: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
385: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
386: } else if (!snes->xl && !snes->xu) {
387: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
388: VecDuplicate(snes->vec_sol, &snes->xl);
389: VecSet(snes->xl,PETSC_NINFINITY);
390: VecDuplicate(snes->vec_sol, &snes->xu);
391: VecSet(snes->xu,PETSC_INFINITY);
392: } else {
393: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
394: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
395: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
396: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
397: 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]))
398: 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.");
399: }
400: }
401: return(0);
402: }
403: /* -------------------------------------------------------------------------- */
404: PetscErrorCode SNESReset_VI(SNES snes)
405: {
409: VecDestroy(&snes->xl);
410: VecDestroy(&snes->xu);
411: snes->usersetbounds = PETSC_FALSE;
412: return(0);
413: }
415: /*
416: SNESDestroy_VI - Destroys the private SNES_VI context that was created
417: with SNESCreate_VI().
419: Input Parameter:
420: . snes - the SNES context
422: Application Interface Routine: SNESDestroy()
423: */
424: PetscErrorCode SNESDestroy_VI(SNES snes)
425: {
429: PetscFree(snes->data);
431: /* clear composed functions */
432: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
433: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetDefaultMonitor_C",NULL);
434: return(0);
435: }
437: /*@
438: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
440: Input Parameters:
441: . snes - the SNES context.
442: . xl - lower bound.
443: . xu - upper bound.
445: Notes:
446: If this routine is not called then the lower and upper bounds are set to
447: PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
449: Level: advanced
451: @*/
452: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
453: {
454: PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
460: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
461: if (!f) {
462: SNESVISetVariableBounds_VI(snes, xl, xu);
463: } else {
464: PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
465: }
466: snes->usersetbounds = PETSC_TRUE;
467: return(0);
468: }
470: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
471: {
472: PetscErrorCode ierr;
473: const PetscScalar *xxl,*xxu;
474: PetscInt i,n, cnt = 0;
477: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
478: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
479: {
480: PetscInt xlN,xuN,N;
481: VecGetSize(xl,&xlN);
482: VecGetSize(xu,&xuN);
483: VecGetSize(snes->vec_func,&N);
484: if (xlN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xlN,N);
485: if (xuN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xuN,N);
486: }
487: PetscObjectReference((PetscObject)xl);
488: PetscObjectReference((PetscObject)xu);
489: VecDestroy(&snes->xl);
490: VecDestroy(&snes->xu);
491: snes->xl = xl;
492: snes->xu = xu;
493: VecGetLocalSize(xl,&n);
494: VecGetArrayRead(xl,&xxl);
495: VecGetArrayRead(xu,&xxu);
496: for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
498: MPIU_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
499: VecRestoreArrayRead(xl,&xxl);
500: VecRestoreArrayRead(xu,&xxu);
501: return(0);
502: }
504: PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems *PetscOptionsObject,SNES snes)
505: {
507: PetscBool flg = PETSC_FALSE;
508: SNESLineSearch linesearch;
511: PetscOptionsHead(PetscOptionsObject,"SNES VI options");
512: PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
513: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
514: if (flg) {
515: SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);
516: }
517: flg = PETSC_FALSE;
518: PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
519: if (flg) {
520: SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
521: }
522: if (!snes->linesearch) {
523: SNESGetLineSearch(snes, &linesearch);
524: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
525: SNESLineSearchBTSetAlpha(linesearch, 0.0);
526: }
527: PetscOptionsTail();
528: return(0);
529: }