Actual source code: vi.c

  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 (*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec));

 21:   PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",&f);
 22:   if (f) PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));
 23:   else SNESVISetComputeVariableBounds_VI(snes,compute);
 24:   return 0;
 25: }

 27: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
 28: {
 29:   snes->ops->computevariablebounds = compute;
 30:   return 0;
 31: }

 33: /* --------------------------------------------------------------------------------------------------------*/

 35: PetscErrorCode  SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
 36: {
 37:   Vec            X, F, Finactive;
 38:   IS             isactive;
 39:   PetscViewer    viewer = (PetscViewer) dummy;

 41:   SNESGetFunction(snes,&F,NULL,NULL);
 42:   SNESGetSolution(snes,&X);
 43:   SNESVIGetActiveSetIS(snes,X,F,&isactive);
 44:   VecDuplicate(F,&Finactive);
 45:   VecCopy(F,Finactive);
 46:   VecISSet(Finactive,isactive,0.0);
 47:   ISDestroy(&isactive);
 48:   VecView(Finactive,viewer);
 49:   VecDestroy(&Finactive);
 50:   return 0;
 51: }

 53: PetscErrorCode  SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
 54: {
 55:   PetscViewer       viewer = (PetscViewer) dummy;
 56:   const PetscScalar *x,*xl,*xu,*f;
 57:   PetscInt          i,n,act[2] = {0,0},fact[2],N;
 58:   /* Number of components that actually hit the bounds (c.f. active variables) */
 59:   PetscInt          act_bound[2] = {0,0},fact_bound[2];
 60:   PetscReal         rnorm,fnorm,zerotolerance = snes->vizerotolerance;
 61:   double            tmp;

 64:   VecGetLocalSize(snes->vec_sol,&n);
 65:   VecGetSize(snes->vec_sol,&N);
 66:   VecGetArrayRead(snes->xl,&xl);
 67:   VecGetArrayRead(snes->xu,&xu);
 68:   VecGetArrayRead(snes->vec_sol,&x);
 69:   VecGetArrayRead(snes->vec_func,&f);

 71:   rnorm = 0.0;
 72:   for (i=0; i<n; i++) {
 73:     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]);
 74:     else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
 75:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
 76:     else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
 77:   }

 79:   for (i=0; i<n; i++) {
 80:     if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
 81:     else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
 82:   }
 83:   VecRestoreArrayRead(snes->vec_func,&f);
 84:   VecRestoreArrayRead(snes->xl,&xl);
 85:   VecRestoreArrayRead(snes->xu,&xu);
 86:   VecRestoreArrayRead(snes->vec_sol,&x);
 87:   MPIU_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
 88:   MPIU_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
 89:   MPIU_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
 90:   fnorm = PetscSqrtReal(fnorm);

 92:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
 93:   if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
 94:   else tmp = 0.0;
 95:   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);

 97:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
 98:   return 0;
 99: }

101: /*
102:      Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
103:     || 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
104:     0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
105:     for this trick. One assumes that the probability that W is in the null space of J is very, very small.
106: */
107: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
108: {
109:   PetscReal      a1;
110:   PetscBool      hastranspose;

112:   *ismin = PETSC_FALSE;
113:   MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
114:   if (hastranspose) {
115:     /* Compute || J^T F|| */
116:     MatMultTranspose(A,F,W);
117:     VecNorm(W,NORM_2,&a1);
118:     PetscInfo(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));
119:     if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
120:   } else {
121:     Vec         work;
122:     PetscScalar result;
123:     PetscReal   wnorm;

125:     VecSetRandom(W,NULL);
126:     VecNorm(W,NORM_2,&wnorm);
127:     VecDuplicate(W,&work);
128:     MatMult(A,W,work);
129:     VecDot(F,work,&result);
130:     VecDestroy(&work);
131:     a1   = PetscAbsScalar(result)/(fnorm*wnorm);
132:     PetscInfo(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);
133:     if (a1 < 1.e-4) *ismin = PETSC_TRUE;
134:   }
135:   return 0;
136: }

138: /*
139:      Checks if J^T(F - J*X) = 0
140: */
141: PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
142: {
143:   PetscReal      a1,a2;
144:   PetscBool      hastranspose;

146:   MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);
147:   if (hastranspose) {
148:     MatMult(A,X,W1);
149:     VecAXPY(W1,-1.0,F);

151:     /* Compute || J^T W|| */
152:     MatMultTranspose(A,W1,W2);
153:     VecNorm(W1,NORM_2,&a1);
154:     VecNorm(W2,NORM_2,&a2);
155:     if (a1 != 0.0) {
156:       PetscInfo(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));
157:     }
158:   }
159:   return 0;
160: }

162: /*
163:   SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.

165:   Notes:
166:   The convergence criterion currently implemented is
167:   merit < abstol
168:   merit < rtol*merit_initial
169: */
170: PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
171: {

175:   *reason = SNES_CONVERGED_ITERATING;

177:   if (!it) {
178:     /* set parameter for default relative tolerance convergence test */
179:     snes->ttol = fnorm*snes->rtol;
180:   }
181:   if (fnorm != fnorm) {
182:     PetscInfo(snes,"Failed to converged, function norm is NaN\n");
183:     *reason = SNES_DIVERGED_FNORM_NAN;
184:   } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
185:     PetscInfo(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);
186:     *reason = SNES_CONVERGED_FNORM_ABS;
187:   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
188:     PetscInfo(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
189:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
190:   }

192:   if (it && !*reason) {
193:     if (fnorm < snes->ttol) {
194:       PetscInfo(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
195:       *reason = SNES_CONVERGED_FNORM_RELATIVE;
196:     }
197:   }
198:   return 0;
199: }

201: /* -------------------------------------------------------------------------- */
202: /*
203:    SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.

205:    Input Parameters:
206: .  SNES - nonlinear solver context

208:    Output Parameters:
209: .  X - Bound projected X

211: */

213: PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
214: {
215:   const PetscScalar *xl,*xu;
216:   PetscScalar       *x;
217:   PetscInt          i,n;

219:   VecGetLocalSize(X,&n);
220:   VecGetArray(X,&x);
221:   VecGetArrayRead(snes->xl,&xl);
222:   VecGetArrayRead(snes->xu,&xu);

224:   for (i = 0; i<n; i++) {
225:     if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
226:     else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
227:   }
228:   VecRestoreArray(X,&x);
229:   VecRestoreArrayRead(snes->xl,&xl);
230:   VecRestoreArrayRead(snes->xu,&xu);
231:   return 0;
232: }

234: /*
235:    SNESVIGetActiveSetIndices - Gets the global indices for the active set variables

237:    Input parameter:
238: .  snes - the SNES context
239: .  X    - the snes solution vector
240: .  F    - the nonlinear function vector

242:    Output parameter:
243: .  ISact - active set index set
244:  */
245: PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact)
246: {
247:   Vec               Xl=snes->xl,Xu=snes->xu;
248:   const PetscScalar *x,*f,*xl,*xu;
249:   PetscInt          *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
250:   PetscReal         zerotolerance = snes->vizerotolerance;

252:   VecGetLocalSize(X,&nlocal);
253:   VecGetOwnershipRange(X,&ilow,&ihigh);
254:   VecGetArrayRead(X,&x);
255:   VecGetArrayRead(Xl,&xl);
256:   VecGetArrayRead(Xu,&xu);
257:   VecGetArrayRead(F,&f);
258:   /* Compute active set size */
259:   for (i=0; i < nlocal;i++) {
260:     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++;
261:   }

263:   PetscMalloc1(nloc_isact,&idx_act);

265:   /* Set active set indices */
266:   for (i=0; i < nlocal; i++) {
267:     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;
268:   }

270:   /* Create active set IS */
271:   ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);

273:   VecRestoreArrayRead(X,&x);
274:   VecRestoreArrayRead(Xl,&xl);
275:   VecRestoreArrayRead(Xu,&xu);
276:   VecRestoreArrayRead(F,&f);
277:   return 0;
278: }

280: PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
281: {
282:   PetscInt       rstart,rend;

284:   SNESVIGetActiveSetIS(snes,X,F,ISact);
285:   VecGetOwnershipRange(X,&rstart,&rend);
286:   ISComplement(*ISact,rstart,rend,ISinact);
287:   return 0;
288: }

290: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
291: {
292:   const PetscScalar *x,*xl,*xu,*f;
293:   PetscInt          i,n;
294:   PetscReal         rnorm,zerotolerance = snes->vizerotolerance;

296:   VecGetLocalSize(X,&n);
297:   VecGetArrayRead(snes->xl,&xl);
298:   VecGetArrayRead(snes->xu,&xu);
299:   VecGetArrayRead(X,&x);
300:   VecGetArrayRead(F,&f);
301:   rnorm = 0.0;
302:   for (i=0; i<n; i++) {
303:     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]);
304:   }
305:   VecRestoreArrayRead(F,&f);
306:   VecRestoreArrayRead(snes->xl,&xl);
307:   VecRestoreArrayRead(snes->xu,&xu);
308:   VecRestoreArrayRead(X,&x);
309:   MPIU_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
310:   *fnorm = PetscSqrtReal(*fnorm);
311:   return 0;
312: }

314: PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
315: {
316:   DMComputeVariableBounds(snes->dm, xl, xu);
317:   return 0;
318: }

320: /* -------------------------------------------------------------------------- */
321: /*
322:    SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
323:    of the SNESVI nonlinear solver.

325:    Input Parameter:
326: .  snes - the SNES context

328:    Application Interface Routine: SNESSetUp()

330:    Notes:
331:    For basic use of the SNES solvers, the user need not explicitly call
332:    SNESSetUp(), since these actions will automatically occur during
333:    the call to SNESSolve().
334:  */
335: PetscErrorCode SNESSetUp_VI(SNES snes)
336: {
337:   PetscInt       i_start[3],i_end[3];

339:   SNESSetWorkVecs(snes,1);
340:   SNESSetUpMatrices(snes);

342:   if (!snes->ops->computevariablebounds && snes->dm) {
343:     PetscBool flag;
344:     DMHasVariableBounds(snes->dm, &flag);
345:     if (flag) {
346:       snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
347:     }
348:   }
349:   if (!snes->usersetbounds) {
350:     if (snes->ops->computevariablebounds) {
351:       if (!snes->xl) VecDuplicate(snes->vec_sol,&snes->xl);
352:       if (!snes->xu) VecDuplicate(snes->vec_sol,&snes->xu);
353:       (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);
354:     } else if (!snes->xl && !snes->xu) {
355:       /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
356:       VecDuplicate(snes->vec_sol, &snes->xl);
357:       VecSet(snes->xl,PETSC_NINFINITY);
358:       VecDuplicate(snes->vec_sol, &snes->xu);
359:       VecSet(snes->xu,PETSC_INFINITY);
360:     } else {
361:       /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
362:       VecGetOwnershipRange(snes->vec_sol,i_start,i_end);
363:       VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);
364:       VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);
365:       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]))
366:         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.");
367:     }
368:   }
369:   return 0;
370: }
371: /* -------------------------------------------------------------------------- */
372: PetscErrorCode SNESReset_VI(SNES snes)
373: {
374:   VecDestroy(&snes->xl);
375:   VecDestroy(&snes->xu);
376:   snes->usersetbounds = PETSC_FALSE;
377:   return 0;
378: }

380: /*
381:    SNESDestroy_VI - Destroys the private SNES_VI context that was created
382:    with SNESCreate_VI().

384:    Input Parameter:
385: .  snes - the SNES context

387:    Application Interface Routine: SNESDestroy()
388:  */
389: PetscErrorCode SNESDestroy_VI(SNES snes)
390: {
391:   PetscFree(snes->data);

393:   /* clear composed functions */
394:   PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);
395:   PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetDefaultMonitor_C",NULL);
396:   return 0;
397: }

399: /*@
400:    SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.

402:    Input Parameters:
403: +  snes - the SNES context.
404: .  xl   - lower bound.
405: -  xu   - upper bound.

407:    Notes:
408:    If this routine is not called then the lower and upper bounds are set to
409:    PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().

411:    Level: advanced

413: @*/
414: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
415: {
416:   PetscErrorCode (*f)(SNES,Vec,Vec);

421:   PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);
422:   if (f) PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));
423:   else SNESVISetVariableBounds_VI(snes, xl, xu);
424:   snes->usersetbounds = PETSC_TRUE;
425:   return 0;
426: }

428: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
429: {
430:   const PetscScalar *xxl,*xxu;
431:   PetscInt          i,n, cnt = 0;

433:   SNESGetFunction(snes,&snes->vec_func,NULL,NULL);
435:   {
436:     PetscInt xlN,xuN,N;
437:     VecGetSize(xl,&xlN);
438:     VecGetSize(xu,&xuN);
439:     VecGetSize(snes->vec_func,&N);
442:   }
443:   PetscObjectReference((PetscObject)xl);
444:   PetscObjectReference((PetscObject)xu);
445:   VecDestroy(&snes->xl);
446:   VecDestroy(&snes->xu);
447:   snes->xl = xl;
448:   snes->xu = xu;
449:   VecGetLocalSize(xl,&n);
450:   VecGetArrayRead(xl,&xxl);
451:   VecGetArrayRead(xu,&xxu);
452:   for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));

454:   MPIU_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));
455:   VecRestoreArrayRead(xl,&xxl);
456:   VecRestoreArrayRead(xu,&xxu);
457:   return 0;
458: }

460: PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems *PetscOptionsObject,SNES snes)
461: {
462:   PetscBool      flg = PETSC_FALSE;

464:   PetscOptionsHead(PetscOptionsObject,"SNES VI options");
465:   PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
466:   PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
467:   if (flg) {
468:     SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);
469:   }
470:   flg = PETSC_FALSE;
471:   PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
472:   if (flg) {
473:     SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
474:   }
475:   PetscOptionsTail();
476:   return 0;
477: }