Actual source code: vi.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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;

510:   PetscOptionsHead(PetscOptionsObject,"SNES VI options");
511:   PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);
512:   PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);
513:   if (flg) {
514:     SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);
515:   }
516:   flg = PETSC_FALSE;
517:   PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);
518:   if (flg) {
519:     SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);
520:   }
521:   PetscOptionsTail();
522:   return(0);
523: }