Actual source code: vi.c
petsc-3.4.5 2014-06-29
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) {SNESSetType(snes,SNESVINEWTONRSLS);}
26: PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));
27: return(0);
28: }
32: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
33: {
35: snes->ops->computevariablebounds = compute;
36: return(0);
37: }
41: /*
42: SNESVIComputeInactiveSetIS - Gets the global indices for the bogus inactive set variables
44: Input parameter
45: . snes - the SNES context
46: . X - the snes solution vector
48: Output parameter
49: . ISact - active set index set
51: */
52: PetscErrorCode SNESVIComputeInactiveSetIS(Vec upper,Vec lower,Vec X,Vec F,IS *inact)
53: {
54: PetscErrorCode ierr;
55: const PetscScalar *x,*xl,*xu,*f;
56: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
59: VecGetLocalSize(X,&nlocal);
60: VecGetOwnershipRange(X,&ilow,&ihigh);
61: VecGetArrayRead(X,&x);
62: VecGetArrayRead(lower,&xl);
63: VecGetArrayRead(upper,&xu);
64: VecGetArrayRead(F,&f);
65: /* Compute inactive set size */
66: for (i=0; i < nlocal; i++) {
67: 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++;
68: }
70: PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);
72: /* Set inactive set indices */
73: for (i=0; i < nlocal; i++) {
74: 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;
75: }
77: /* Create inactive set IS */
78: ISCreateGeneral(PetscObjectComm((PetscObject)upper),nloc_isact,idx_act,PETSC_OWN_POINTER,inact);
80: VecRestoreArrayRead(X,&x);
81: VecRestoreArrayRead(lower,&xl);
82: VecRestoreArrayRead(upper,&xu);
83: VecRestoreArrayRead(F,&f);
84: return(0);
85: }
87: /* --------------------------------------------------------------------------------------------------------*/
91: PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
92: {
93: PetscErrorCode ierr;
94: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
95: const PetscScalar *x,*xl,*xu,*f;
96: PetscInt i,n,act[2] = {0,0},fact[2],N;
97: /* Number of components that actually hit the bounds (c.f. active variables) */
98: PetscInt act_bound[2] = {0,0},fact_bound[2];
99: PetscReal rnorm,fnorm;
100: double tmp;
103: VecGetLocalSize(snes->vec_sol,&n);
104: VecGetSize(snes->vec_sol,&N);
105: VecGetArrayRead(snes->xl,&xl);
106: VecGetArrayRead(snes->xu,&xu);
107: VecGetArrayRead(snes->vec_sol,&x);
108: VecGetArrayRead(snes->vec_func,&f);
110: rnorm = 0.0;
111: for (i=0; i<n; i++) {
112: 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]);
113: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8 && PetscRealPart(f[i]) >= 0.0) act[0]++;
114: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8 && PetscRealPart(f[i]) <= 0.0) act[1]++;
115: else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
116: }
118: for (i=0; i<n; i++) {
119: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + 1.e-8) act_bound[0]++;
120: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - 1.e-8) act_bound[1]++;
121: }
122: VecRestoreArrayRead(snes->vec_func,&f);
123: VecRestoreArrayRead(snes->xl,&xl);
124: VecRestoreArrayRead(snes->xu,&xu);
125: VecRestoreArrayRead(snes->vec_sol,&x);
126: MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
127: MPI_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
128: MPI_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
129: fnorm = PetscSqrtReal(fnorm);
131: PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
132: if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
133: else tmp = 0.0;
134: 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);
136: PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
137: return(0);
138: }
140: /*
141: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
142: || 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
143: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
144: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
145: */
148: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
149: {
150: PetscReal a1;
152: PetscBool hastranspose;
155: *ismin = PETSC_FALSE;
156: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
157: if (hastranspose) {
158: /* Compute || J^T F|| */
159: MatMultTranspose(A,F,W);
160: VecNorm(W,NORM_2,&a1);
161: PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
162: if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
163: } else {
164: Vec work;
165: PetscScalar result;
166: PetscReal wnorm;
168: VecSetRandom(W,NULL);
169: VecNorm(W,NORM_2,&wnorm);
170: VecDuplicate(W,&work);
171: MatMult(A,W,work);
172: VecDot(F,work,&result);
173: VecDestroy(&work);
174: a1 = PetscAbsScalar(result)/(fnorm*wnorm);
175: PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
176: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
177: }
178: return(0);
179: }
181: /*
182: Checks if J^T(F - J*X) = 0
183: */
186: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
187: {
188: PetscReal a1,a2;
190: PetscBool hastranspose;
193: MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
194: if (hastranspose) {
195: MatMult(A,X,W1);
196: VecAXPY(W1,-1.0,F);
198: /* Compute || J^T W|| */
199: MatMultTranspose(A,W1,W2);
200: VecNorm(W1,NORM_2,&a1);
201: VecNorm(W2,NORM_2,&a2);
202: if (a1 != 0.0) {
203: PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
204: }
205: }
206: return(0);
207: }
209: /*
210: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
212: Notes:
213: The convergence criterion currently implemented is
214: merit < abstol
215: merit < rtol*merit_initial
216: */
219: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
220: {
227: *reason = SNES_CONVERGED_ITERATING;
229: if (!it) {
230: /* set parameter for default relative tolerance convergence test */
231: snes->ttol = fnorm*snes->rtol;
232: }
233: if (fnorm != fnorm) {
234: PetscInfo(snes,"Failed to converged, function norm is NaN\n");
235: *reason = SNES_DIVERGED_FNORM_NAN;
236: } else if (fnorm < snes->abstol) {
237: PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
238: *reason = SNES_CONVERGED_FNORM_ABS;
239: } else if (snes->nfuncs >= snes->max_funcs) {
240: PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
241: *reason = SNES_DIVERGED_FUNCTION_COUNT;
242: }
244: if (it && !*reason) {
245: if (fnorm < snes->ttol) {
246: PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
247: *reason = SNES_CONVERGED_FNORM_RELATIVE;
248: }
249: }
250: return(0);
251: }
254: /* -------------------------------------------------------------------------- */
255: /*
256: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
258: Input Parameters:
259: . SNES - nonlinear solver context
261: Output Parameters:
262: . X - Bound projected X
264: */
268: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
269: {
270: PetscErrorCode ierr;
271: const PetscScalar *xl,*xu;
272: PetscScalar *x;
273: PetscInt i,n;
276: VecGetLocalSize(X,&n);
277: VecGetArray(X,&x);
278: VecGetArrayRead(snes->xl,&xl);
279: VecGetArrayRead(snes->xu,&xu);
281: for (i = 0; i<n; i++) {
282: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
283: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
284: }
285: VecRestoreArray(X,&x);
286: VecRestoreArrayRead(snes->xl,&xl);
287: VecRestoreArrayRead(snes->xu,&xu);
288: return(0);
289: }
294: /*
295: SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
297: Input parameter
298: . snes - the SNES context
299: . X - the snes solution vector
300: . F - the nonlinear function vector
302: Output parameter
303: . ISact - active set index set
304: */
305: PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact)
306: {
307: PetscErrorCode ierr;
308: Vec Xl=snes->xl,Xu=snes->xu;
309: const PetscScalar *x,*f,*xl,*xu;
310: PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
313: VecGetLocalSize(X,&nlocal);
314: VecGetOwnershipRange(X,&ilow,&ihigh);
315: VecGetArrayRead(X,&x);
316: VecGetArrayRead(Xl,&xl);
317: VecGetArrayRead(Xu,&xu);
318: VecGetArrayRead(F,&f);
319: /* Compute active set size */
320: for (i=0; i < nlocal;i++) {
321: 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++;
322: }
324: PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);
326: /* Set active set indices */
327: for (i=0; i < nlocal; i++) {
328: 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;
329: }
331: /* Create active set IS */
332: ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);
334: VecRestoreArrayRead(X,&x);
335: VecRestoreArrayRead(Xl,&xl);
336: VecRestoreArrayRead(Xu,&xu);
337: VecRestoreArrayRead(F,&f);
338: return(0);
339: }
343: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
344: {
346: PetscInt rstart,rend;
349: SNESVIGetActiveSetIS(snes,X,F,ISact);
350: VecGetOwnershipRange(X,&rstart,&rend);
351: ISComplement(*ISact,rstart,rend,ISinact);
352: return(0);
353: }
357: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
358: {
359: PetscErrorCode ierr;
360: const PetscScalar *x,*xl,*xu,*f;
361: PetscInt i,n;
362: PetscReal rnorm;
365: VecGetLocalSize(X,&n);
366: VecGetArrayRead(snes->xl,&xl);
367: VecGetArrayRead(snes->xu,&xu);
368: VecGetArrayRead(X,&x);
369: VecGetArrayRead(F,&f);
370: rnorm = 0.0;
371: for (i=0; i<n; i++) {
372: 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]);
373: }
374: VecRestoreArrayRead(F,&f);
375: VecRestoreArrayRead(snes->xl,&xl);
376: VecRestoreArrayRead(snes->xu,&xu);
377: VecRestoreArrayRead(X,&x);
378: MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
379: *fnorm = PetscSqrtReal(*fnorm);
380: return(0);
381: }
385: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
386: {
390: DMComputeVariableBounds(snes->dm, xl, xu);
391: return(0);
392: }
395: /* -------------------------------------------------------------------------- */
396: /*
397: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
398: of the SNESVI nonlinear solver.
400: Input Parameter:
401: . snes - the SNES context
402: . x - the solution vector
404: Application Interface Routine: SNESSetUp()
406: Notes:
407: For basic use of the SNES solvers, the user need not explicitly call
408: SNESSetUp(), since these actions will automatically occur during
409: the call to SNESSolve().
410: */
413: PetscErrorCode SNESSetUp_VI(SNES snes)
414: {
416: PetscInt i_start[3],i_end[3];
419: SNESSetWorkVecs(snes,1);
420: SNESSetUpMatrices(snes);
422: if (!snes->ops->computevariablebounds && snes->dm) {
423: PetscBool flag;
424: DMHasVariableBounds(snes->dm, &flag);
426: snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
427: }
428: if (!snes->usersetbounds) {
429: if (snes->ops->computevariablebounds) {
430: if (!snes->xl) {VecDuplicate(snes->vec_sol,&snes->xl);}
431: if (!snes->xu) {VecDuplicate(snes->vec_sol,&snes->xu);}
432: (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
433: } else if (!snes->xl && !snes->xu) {
434: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
435: VecDuplicate(snes->vec_sol, &snes->xl);
436: VecSet(snes->xl,SNES_VI_NINF);
437: VecDuplicate(snes->vec_sol, &snes->xu);
438: VecSet(snes->xu,SNES_VI_INF);
439: } else {
440: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
441: VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
442: VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
443: VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
444: 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]))
445: 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.");
446: }
447: }
448: return(0);
449: }
450: /* -------------------------------------------------------------------------- */
453: PetscErrorCode SNESReset_VI(SNES snes)
454: {
458: VecDestroy(&snes->xl);
459: VecDestroy(&snes->xu);
460: snes->usersetbounds = PETSC_FALSE;
461: return(0);
462: }
464: /*
465: SNESDestroy_VI - Destroys the private SNES_VI context that was created
466: with SNESCreate_VI().
468: Input Parameter:
469: . snes - the SNES context
471: Application Interface Routine: SNESDestroy()
472: */
475: PetscErrorCode SNESDestroy_VI(SNES snes)
476: {
480: PetscFree(snes->data);
482: /* clear composed functions */
483: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
484: PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetMonitor_C",NULL);
485: return(0);
486: }
490: /*@
491: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
493: Input Parameters:
494: . snes - the SNES context.
495: . xl - lower bound.
496: . xu - upper bound.
498: Notes:
499: If this routine is not called then the lower and upper bounds are set to
500: SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp().
502: Level: advanced
504: @*/
505: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
506: {
507: PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
513: PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
514: if (!f) {SNESSetType(snes,SNESVINEWTONRSLS);}
515: PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
516: snes->usersetbounds = PETSC_TRUE;
517: return(0);
518: }
522: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
523: {
524: PetscErrorCode ierr;
525: const PetscScalar *xxl,*xxu;
526: PetscInt i,n, cnt = 0;
529: SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
530: if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
531: {
532: PetscInt xlN,xuN,N;
533: VecGetSize(xl,&xlN);
534: VecGetSize(xu,&xuN);
535: VecGetSize(snes->vec_func,&N);
536: if (xlN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xlN,N);
537: if (xuN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xuN,N);
538: }
539: SNESSetType(snes,SNESVINEWTONRSLS);
540: PetscObjectReference((PetscObject)xl);
541: PetscObjectReference((PetscObject)xu);
542: VecDestroy(&snes->xl);
543: VecDestroy(&snes->xu);
544: snes->xl = xl;
545: snes->xu = xu;
546: VecGetLocalSize(xl,&n);
547: VecGetArrayRead(xl,&xxl);
548: VecGetArrayRead(xu,&xxu);
549: for (i=0; i<n; i++) cnt += ((xxl[i] != SNES_VI_NINF) || (xxu[i] != SNES_VI_INF));
551: MPI_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
552: VecRestoreArrayRead(xl,&xxl);
553: VecRestoreArrayRead(xu,&xxu);
554: return(0);
555: }
559: PetscErrorCode SNESSetFromOptions_VI(SNES snes)
560: {
562: PetscBool flg;
563: SNESLineSearch linesearch;
566: PetscOptionsHead("SNES VI options");
567: PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);
568: if (flg) {
569: SNESMonitorSet(snes,SNESMonitorVI,0,0);
570: }
571: if (!snes->linesearch) {
572: SNESGetLineSearch(snes, &linesearch);
573: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
574: SNESLineSearchBTSetAlpha(linesearch, 0.0);
575: }
576: PetscOptionsTail();
577: return(0);
578: }