Actual source code: vi.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
5: /*@C
6: SNESVISetComputeVariableBounds - Sets a function that is called to compute the variable bounds
8: Input parameter
9: + snes - the SNES context
10: - compute - computes the bounds
12: Level: advanced
14: @*/
15: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES,Vec,Vec))
16: {
17: PetscErrorCode ierr,(*f)(SNES,PetscErrorCode(*)(SNES,Vec,Vec));
21: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",(PetscVoidStarFunction)&f);
22: if (!f) {SNESSetType(snes,SNESVIRS);}
23: PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode(*)(SNES,Vec,Vec)),(snes,compute));
24: return(0);
25: }
27: EXTERN_C_BEGIN
30: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
31: {
33: snes->ops->computevariablebounds = compute;
34: return(0);
35: }
36: EXTERN_C_END
40: /*
41: SNESVIComputeInactiveSetIS - Gets the global indices for the bogus inactive set variables
43: Input parameter
44: . snes - the SNES context
45: . X - the snes solution vector
47: Output parameter
48: . ISact - active set index set
50: */
51: PetscErrorCode SNESVIComputeInactiveSetIS(Vec upper,Vec lower,Vec X,Vec F,IS* inact)
52: {
53: PetscErrorCode ierr;
54: const PetscScalar *x,*xl,*xu,*f;
55: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
56:
58: VecGetLocalSize(X,&nlocal);
59: VecGetOwnershipRange(X,&ilow,&ihigh);
60: VecGetArrayRead(X,&x);
61: VecGetArrayRead(lower,&xl);
62: VecGetArrayRead(upper,&xu);
63: VecGetArrayRead(F,&f);
64: /* Compute inactive set size */
65: for (i=0; i < nlocal;i++) {
66: if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++;
67: }
69: PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);
71: /* Set inactive set indices */
72: for(i=0; i < nlocal; i++) {
73: if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i;
74: }
76: /* Create inactive set IS */
77: ISCreateGeneral(((PetscObject)upper)->comm,nloc_isact,idx_act,PETSC_OWN_POINTER,inact);
79: VecRestoreArrayRead(X,&x);
80: VecRestoreArrayRead(lower,&xl);
81: VecRestoreArrayRead(upper,&xu);
82: VecRestoreArrayRead(F,&f);
83: return(0);
84: }
86: /* --------------------------------------------------------------------------------------------------------*/
90: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
91: {
92: PetscErrorCode ierr;
93: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);
94: const PetscScalar *x,*xl,*xu,*f;
95: PetscInt i,n,act[2] = {0,0},fact[2],N;
96: /* Number of components that actually hit the bounds (c.f. active variables) */
97: PetscInt act_bound[2] = {0,0},fact_bound[2];
98: PetscReal rnorm,fnorm;
99: double tmp;
102: VecGetLocalSize(snes->vec_sol,&n);
103: VecGetSize(snes->vec_sol,&N);
104: VecGetArrayRead(snes->xl,&xl);
105: VecGetArrayRead(snes->xu,&xu);
106: VecGetArrayRead(snes->vec_sol,&x);
107: VecGetArrayRead(snes->vec_func,&f);
108:
109: rnorm = 0.0;
110: for (i=0; i<n; i++) {
111: if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]);
112: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8 && PetscRealPart(f[i]) >= 0.0) act[0]++;
113: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8 && PetscRealPart(f[i]) <= 0.0) act[1]++;
114: else SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_PLIB,"Can never get here");
115: }
117: for (i=0; i<n; i++) {
118: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8) act_bound[0]++;
119: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8) act_bound[1]++;
120: }
121: VecRestoreArrayRead(snes->vec_func,&f);
122: VecRestoreArrayRead(snes->xl,&xl);
123: VecRestoreArrayRead(snes->xu,&xu);
124: VecRestoreArrayRead(snes->vec_sol,&x);
125: MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);
126: MPI_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,((PetscObject)snes)->comm);
127: MPI_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,((PetscObject)snes)->comm);
128: fnorm = PetscSqrtReal(fnorm);
129:
130: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
131: if(snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
132: else tmp = 0.0;
133: 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);
134:
135: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
136: return(0);
137: }
139: /*
140: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
141: || 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
142: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
143: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
144: */
147: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
148: {
149: PetscReal a1;
151: PetscBool hastranspose;
154: *ismin = PETSC_FALSE;
155: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
156: if (hastranspose) {
157: /* Compute || J^T F|| */
158: MatMultTranspose(A,F,W);
159: VecNorm(W,NORM_2,&a1);
160: PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
161: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
162: } else {
163: Vec work;
164: PetscScalar result;
165: PetscReal wnorm;
167: VecSetRandom(W,PETSC_NULL);
168: VecNorm(W,NORM_2,&wnorm);
169: VecDuplicate(W,&work);
170: MatMult(A,W,work);
171: VecDot(F,work,&result);
172: VecDestroy(&work);
173: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
174: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
175: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
176: }
177: return(0);
178: }
180: /*
181: Checks if J^T(F - J*X) = 0
182: */
185: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
186: {
187: PetscReal a1,a2;
189: PetscBool hastranspose;
192: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
193: if (hastranspose) {
194: MatMult(A,X,W1);
195: VecAXPY(W1,-1.0,F);
197: /* Compute || J^T W|| */
198: MatMultTranspose(A,W1,W2);
199: VecNorm(W1,NORM_2,&a1);
200: VecNorm(W2,NORM_2,&a2);
201: if (a1 != 0.0) {
202: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
203: }
204: }
205: return(0);
206: }
208: /*
209: SNESDefaultConverged_VI - Checks the convergence of the semismooth newton algorithm.
211: Notes:
212: The convergence criterion currently implemented is
213: merit < abstol
214: merit < rtol*merit_initial
215: */
218: PetscErrorCode SNESDefaultConverged_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
219: {
225:
226: *reason = SNES_CONVERGED_ITERATING;
228: if (!it) {
229: /* set parameter for default relative tolerance convergence test */
230: snes->ttol = fnorm*snes->rtol;
231: }
232: if (fnorm != fnorm) {
233: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
234: *reason = SNES_DIVERGED_FNORM_NAN;
235: } else if (fnorm < snes->abstol) {
236: PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
237: *reason = SNES_CONVERGED_FNORM_ABS;
238: } else if (snes->nfuncs >= snes->max_funcs) {
239: PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
240: *reason = SNES_DIVERGED_FUNCTION_COUNT;
241: }
243: if (it && !*reason) {
244: if (fnorm < snes->ttol) {
245: PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
246: *reason = SNES_CONVERGED_FNORM_RELATIVE;
247: }
248: }
249: return(0);
250: }
253: /* -------------------------------------------------------------------------- */
254: /*
255: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
257: Input Parameters:
258: . SNES - nonlinear solver context
260: Output Parameters:
261: . X - Bound projected X
263: */
267: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
268: {
269: PetscErrorCode ierr;
270: const PetscScalar *xl,*xu;
271: PetscScalar *x;
272: PetscInt i,n;
275: VecGetLocalSize(X,&n);
276: VecGetArray(X,&x);
277: VecGetArrayRead(snes->xl,&xl);
278: VecGetArrayRead(snes->xu,&xu);
280: for(i = 0;i<n;i++) {
281: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
282: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
283: }
284: VecRestoreArray(X,&x);
285: VecRestoreArrayRead(snes->xl,&xl);
286: VecRestoreArrayRead(snes->xu,&xu);
287: return(0);
288: }
293: /*
294: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
296: Input parameter
297: . snes - the SNES context
298: . X - the snes solution vector
299: . F - the nonlinear function vector
301: Output parameter
302: . ISact - active set index set
303: */
304: PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS* ISact)
305: {
306: PetscErrorCode ierr;
307: Vec Xl=snes->xl,Xu=snes->xu;
308: const PetscScalar *x,*f,*xl,*xu;
309: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
310:
312: VecGetLocalSize(X,&nlocal);
313: VecGetOwnershipRange(X,&ilow,&ihigh);
314: VecGetArrayRead(X,&x);
315: VecGetArrayRead(Xl,&xl);
316: VecGetArrayRead(Xu,&xu);
317: VecGetArrayRead(F,&f);
318: /* Compute active set size */
319: for (i=0; i < nlocal;i++) {
320: if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) nloc_isact++;
321: }
323: PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);
325: /* Set active set indices */
326: for(i=0; i < nlocal; i++) {
327: if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) idx_act[i1++] = ilow+i;
328: }
330: /* Create active set IS */
331: ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
333: VecRestoreArrayRead(X,&x);
334: VecRestoreArrayRead(Xl,&xl);
335: VecRestoreArrayRead(Xu,&xu);
336: VecRestoreArrayRead(F,&f);
337: return(0);
338: }
342: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS* ISact,IS* ISinact)
343: {
344: PetscErrorCode ierr;
347: SNESVIGetActiveSetIS(snes,X,F,ISact);
348: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
349: return(0);
350: }
354: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
355: {
356: PetscErrorCode ierr;
357: const PetscScalar *x,*xl,*xu,*f;
358: PetscInt i,n;
359: PetscReal rnorm;
362: VecGetLocalSize(X,&n);
363: VecGetArrayRead(snes->xl,&xl);
364: VecGetArrayRead(snes->xu,&xu);
365: VecGetArrayRead(X,&x);
366: VecGetArrayRead(F,&f);
367: rnorm = 0.0;
368: for (i=0; i<n; i++) {
369: if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + 1.e-8 || (PetscRealPart(f[i]) < 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - 1.e-8) || PetscRealPart(f[i]) > 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]);
370: }
371: VecRestoreArrayRead(F,&f);
372: VecRestoreArrayRead(snes->xl,&xl);
373: VecRestoreArrayRead(snes->xu,&xu);
374: VecRestoreArrayRead(X,&x);
375: MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,((PetscObject)snes)->comm);
376: *fnorm = PetscSqrtReal(*fnorm);
377: return(0);
378: }
382: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
383: {
384: PetscErrorCode ierr;
386: DMComputeVariableBounds(snes->dm, xl, xu);
387: return(0);
388: }
391: /* -------------------------------------------------------------------------- */
392: /*
393: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
394: of the SNESVI nonlinear solver.
396: Input Parameter:
397: . snes - the SNES context
398: . x - the solution vector
400: Application Interface Routine: SNESSetUp()
402: Notes:
403: For basic use of the SNES solvers, the user need not explicitly call
404: SNESSetUp(), since these actions will automatically occur during
405: the call to SNESSolve().
406: */
409: PetscErrorCode SNESSetUp_VI(SNES snes)
410: {
412: PetscInt i_start[3],i_end[3];
416: SNESDefaultGetWork(snes,1);
417: SNESSetUpMatrices(snes);
419: if(!snes->ops->computevariablebounds && snes->dm) {
420: PetscBool flag;
421: DMHasVariableBounds(snes->dm, &flag);
422: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
423: }
424: if(!snes->usersetbounds) {
425: if (snes->ops->computevariablebounds) {
426: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
427: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
428: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
429: }
430: else if(!snes->xl && !snes->xu) {
431: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
432: VecDuplicate(snes->vec_sol, &snes->xl);
433: VecSet(snes->xl,SNES_VI_NINF);
434: VecDuplicate(snes->vec_sol, &snes->xu);
435: VecSet(snes->xu,SNES_VI_INF);
436: } else {
437: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
438: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
439: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
440: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
441: 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]))
442: 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.");
443: }
444: }
446: return(0);
447: }
448: /* -------------------------------------------------------------------------- */
451: PetscErrorCode SNESReset_VI(SNES snes)
452: {
456: VecDestroy(&snes->xl);
457: VecDestroy(&snes->xu);
458: snes->usersetbounds = PETSC_FALSE;
459: return(0);
460: }
462: /*
463: SNESDestroy_VI - Destroys the private SNES_VI context that was created
464: with SNESCreate_VI().
466: Input Parameter:
467: . snes - the SNES context
469: Application Interface Routine: SNESDestroy()
470: */
473: PetscErrorCode SNESDestroy_VI(SNES snes)
474: {
478: PetscFree(snes->data);
480: /* clear composed functions */
481: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);
482: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);
483: return(0);
484: }
488: /*@
489: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
491: Input Parameters:
492: . snes - the SNES context.
493: . xl - lower bound.
494: . xu - upper bound.
496: Notes:
497: If this routine is not called then the lower and upper bounds are set to
498: SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp().
500: Level: advanced
502: @*/
503: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
504: {
505: PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
511: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",(PetscVoidStarFunction)&f);
512: if (!f) {SNESSetType(snes,SNESVIRS);}
513: PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
514: snes->usersetbounds = PETSC_TRUE;
515: return(0);
516: }
518: EXTERN_C_BEGIN
521: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
522: {
523: PetscErrorCode ierr;
524: const PetscScalar *xxl,*xxu;
525: PetscInt i,n, cnt = 0;
528: SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);
529: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
530: if (xl->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xl->map->N,snes->vec_func->map->N);
531: if (xu->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xu->map->N,snes->vec_func->map->N);
532: SNESSetType(snes,SNESVIRS);
533: PetscObjectReference((PetscObject)xl);
534: PetscObjectReference((PetscObject)xu);
535: VecDestroy(&snes->xl);
536: VecDestroy(&snes->xu);
537: snes->xl = xl;
538: snes->xu = xu;
539: VecGetLocalSize(xl,&n);
540: VecGetArrayRead(xl,&xxl);
541: VecGetArrayRead(xu,&xxu);
542: for (i=0; i<n; i++) {
543: cnt += ((xxl[i] != SNES_VI_NINF) || (xxu[i] != SNES_VI_INF));
544: }
545: MPI_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,((PetscObject)snes)->comm);
546: VecRestoreArrayRead(xl,&xxl);
547: VecRestoreArrayRead(xu,&xxu);
548: return(0);
549: }
550: EXTERN_C_END
554: PetscErrorCode SNESSetFromOptions_VI(SNES snes)
555: {
556: PetscErrorCode ierr;
557: PetscBool flg;
558: SNESLineSearch linesearch;
561: PetscOptionsHead("SNES VI options");
562: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);
563: if (flg) {
564: SNESMonitorSet(snes,SNESMonitorVI,0,0);
565: }
566: if (!snes->linesearch) {
567: SNESGetSNESLineSearch(snes, &linesearch);
568: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
569: SNESLineSearchBTSetAlpha(linesearch, 0.0);
570: }
571: PetscOptionsTail();
572: return(0);
573: }