Actual source code: taocg.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/unconstrained/impls/cg/taocg.h>

  4: #define CG_FletcherReeves       0
  5: #define CG_PolakRibiere         1
  6: #define CG_PolakRibierePlus     2
  7: #define CG_HestenesStiefel      3
  8: #define CG_DaiYuan              4
  9: #define CG_Types                5

 11: static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};

 13: static PetscErrorCode TaoSolve_CG(Tao tao)
 14: {
 15:   TAO_CG                       *cgP = (TAO_CG*)tao->data;
 16:   PetscErrorCode               ierr;
 17:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 18:   PetscReal                    step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta;
 19:   PetscReal                    gd_old,gnorm2_old,f_old;

 22:   if (tao->XL || tao->XU || tao->ops->computebounds) {
 23:     PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
 24:   }

 26:   /*  Check convergence criteria */
 27:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 28:   VecNorm(tao->gradient,NORM_2,&gnorm);
 29:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

 31:   tao->reason = TAO_CONTINUE_ITERATING;
 32:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 33:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
 34:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 35:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

 37:   /*  Set initial direction to -gradient */
 38:   VecCopy(tao->gradient, tao->stepdirection);
 39:   VecScale(tao->stepdirection, -1.0);
 40:   gnorm2 = gnorm*gnorm;

 42:   /*  Set initial scaling for the function */
 43:   if (f != 0.0) {
 44:     delta = 2.0*PetscAbsScalar(f) / gnorm2;
 45:     delta = PetscMax(delta,cgP->delta_min);
 46:     delta = PetscMin(delta,cgP->delta_max);
 47:   } else {
 48:     delta = 2.0 / gnorm2;
 49:     delta = PetscMax(delta,cgP->delta_min);
 50:     delta = PetscMin(delta,cgP->delta_max);
 51:   }
 52:   /*  Set counter for gradient and reset steps */
 53:   cgP->ngradsteps = 0;
 54:   cgP->nresetsteps = 0;

 56:   while (1) {
 57:     /* Call general purpose update function */
 58:     if (tao->ops->update) {
 59:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 60:     }

 62:     /*  Save the current gradient information */
 63:     f_old = f;
 64:     gnorm2_old = gnorm2;
 65:     VecCopy(tao->solution, cgP->X_old);
 66:     VecCopy(tao->gradient, cgP->G_old);
 67:     VecDot(tao->gradient, tao->stepdirection, &gd);
 68:     if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
 69:       ++cgP->ngradsteps;
 70:       if (f != 0.0) {
 71:         delta = 2.0*PetscAbsScalar(f) / gnorm2;
 72:         delta = PetscMax(delta,cgP->delta_min);
 73:         delta = PetscMin(delta,cgP->delta_max);
 74:       } else {
 75:         delta = 2.0 / gnorm2;
 76:         delta = PetscMax(delta,cgP->delta_min);
 77:         delta = PetscMin(delta,cgP->delta_max);
 78:       }

 80:       VecCopy(tao->gradient, tao->stepdirection);
 81:       VecScale(tao->stepdirection, -1.0);
 82:     }

 84:     /*  Search direction for improving point */
 85:     TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
 86:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
 87:     TaoAddLineSearchCounts(tao);
 88:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 89:       /*  Linesearch failed */
 90:       /*  Reset factors and use scaled gradient step */
 91:       ++cgP->nresetsteps;
 92:       f = f_old;
 93:       gnorm2 = gnorm2_old;
 94:       VecCopy(cgP->X_old, tao->solution);
 95:       VecCopy(cgP->G_old, tao->gradient);

 97:       if (f != 0.0) {
 98:         delta = 2.0*PetscAbsScalar(f) / gnorm2;
 99:         delta = PetscMax(delta,cgP->delta_min);
100:         delta = PetscMin(delta,cgP->delta_max);
101:       } else {
102:         delta = 2.0 / gnorm2;
103:         delta = PetscMax(delta,cgP->delta_min);
104:         delta = PetscMin(delta,cgP->delta_max);
105:       }

107:       VecCopy(tao->gradient, tao->stepdirection);
108:       VecScale(tao->stepdirection, -1.0);

110:       TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
111:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
112:       TaoAddLineSearchCounts(tao);

114:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
115:         /*  Linesearch failed again */
116:         /*  switch to unscaled gradient */
117:         f = f_old;
118:         VecCopy(cgP->X_old, tao->solution);
119:         VecCopy(cgP->G_old, tao->gradient);
120:         delta = 1.0;
121:         VecCopy(tao->solution, tao->stepdirection);
122:         VecScale(tao->stepdirection, -1.0);

124:         TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
125:         TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
126:         TaoAddLineSearchCounts(tao);
127:         if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {

129:           /*  Line search failed for last time -- give up */
130:           f = f_old;
131:           VecCopy(cgP->X_old, tao->solution);
132:           VecCopy(cgP->G_old, tao->gradient);
133:           step = 0.0;
134:           tao->reason = TAO_DIVERGED_LS_FAILURE;
135:         }
136:       }
137:     }

139:     /*  Check for bad value */
140:     VecNorm(tao->gradient,NORM_2,&gnorm);
141:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER,"User-provided compute function generated Inf or NaN");

143:     /*  Check for termination */
144:     gnorm2 =gnorm * gnorm;
145:     tao->niter++;
146:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
147:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
148:     (*tao->ops->convergencetest)(tao,tao->cnvP);
149:     if (tao->reason != TAO_CONTINUE_ITERATING) {
150:       break;
151:     }

153:     /*  Check for restart condition */
154:     VecDot(tao->gradient, cgP->G_old, &ginner);
155:     if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
156:       /*  Gradients far from orthognal; use steepest descent direction */
157:       beta = 0.0;
158:     } else {
159:       /*  Gradients close to orthogonal; use conjugate gradient formula */
160:       switch (cgP->cg_type) {
161:       case CG_FletcherReeves:
162:         beta = gnorm2 / gnorm2_old;
163:         break;

165:       case CG_PolakRibiere:
166:         beta = (gnorm2 - ginner) / gnorm2_old;
167:         break;

169:       case CG_PolakRibierePlus:
170:         beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
171:         break;

173:       case CG_HestenesStiefel:
174:         VecDot(tao->gradient, tao->stepdirection, &gd);
175:         VecDot(cgP->G_old, tao->stepdirection, &gd_old);
176:         beta = (gnorm2 - ginner) / (gd - gd_old);
177:         break;

179:       case CG_DaiYuan:
180:         VecDot(tao->gradient, tao->stepdirection, &gd);
181:         VecDot(cgP->G_old, tao->stepdirection, &gd_old);
182:         beta = gnorm2 / (gd - gd_old);
183:         break;

185:       default:
186:         beta = 0.0;
187:         break;
188:       }
189:     }

191:     /*  Compute the direction d=-g + beta*d */
192:     VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);

194:     /*  update initial steplength choice */
195:     delta = 1.0;
196:     delta = PetscMax(delta, cgP->delta_min);
197:     delta = PetscMin(delta, cgP->delta_max);
198:   }
199:   return(0);
200: }

202: static PetscErrorCode TaoSetUp_CG(Tao tao)
203: {
204:   TAO_CG         *cgP = (TAO_CG*)tao->data;

208:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
209:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
210:   if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old);}
211:   if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old); }
212:   return(0);
213: }

215: static PetscErrorCode TaoDestroy_CG(Tao tao)
216: {
217:   TAO_CG         *cgP = (TAO_CG*) tao->data;

221:   if (tao->setupcalled) {
222:     VecDestroy(&cgP->X_old);
223:     VecDestroy(&cgP->G_old);
224:   }
225:   TaoLineSearchDestroy(&tao->linesearch);
226:   PetscFree(tao->data);
227:   return(0);
228: }

230: static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
231:  {
232:     TAO_CG         *cgP = (TAO_CG*)tao->data;

236:     TaoLineSearchSetFromOptions(tao->linesearch);
237:     PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
238:     PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);
239:     PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);
240:     PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);
241:     PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);
242:    PetscOptionsTail();
243:    return(0);
244: }

246: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
247: {
248:   PetscBool      isascii;
249:   TAO_CG         *cgP = (TAO_CG*)tao->data;

253:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
254:   if (isascii) {
255:     PetscViewerASCIIPushTab(viewer);
256:     PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
257:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
258:     ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
259:     PetscViewerASCIIPopTab(viewer);
260:   }
261:   return(0);
262: }

264: /*MC
265:      TAOCG -   Nonlinear conjugate gradient method is an extension of the
266: nonlinear conjugate gradient solver for nonlinear optimization.

268:    Options Database Keys:
269: +      -tao_cg_eta <r> - restart tolerance
270: .      -tao_cg_type <taocg_type> - cg formula
271: .      -tao_cg_delta_min <r> - minimum delta value
272: -      -tao_cg_delta_max <r> - maximum delta value

274:   Notes:
275:      CG formulas are:
276:          "fr" - Fletcher-Reeves
277:          "pr" - Polak-Ribiere
278:          "prp" - Polak-Ribiere-Plus
279:          "hs" - Hestenes-Steifel
280:          "dy" - Dai-Yuan
281:   Level: beginner
282: M*/


285: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
286: {
287:   TAO_CG         *cgP;
288:   const char     *morethuente_type = TAOLINESEARCHMT;

292:   tao->ops->setup = TaoSetUp_CG;
293:   tao->ops->solve = TaoSolve_CG;
294:   tao->ops->view = TaoView_CG;
295:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
296:   tao->ops->destroy = TaoDestroy_CG;

298:   /* Override default settings (unless already changed) */
299:   if (!tao->max_it_changed) tao->max_it = 2000;
300:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

302:   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
303:   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
304:   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
305:   /*  linesearch because it seems to work better. */
306:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
307:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
308:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
309:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
310:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

312:   PetscNewLog(tao,&cgP);
313:   tao->data = (void*)cgP;
314:   cgP->eta = 0.1;
315:   cgP->delta_min = 1e-7;
316:   cgP->delta_max = 100;
317:   cgP->cg_type = CG_PolakRibierePlus;
318:   return(0);
319: }