Actual source code: tron.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <../src/tao/bound/impls/tron/tron.h>
  2:  #include <../src/tao/matrix/submatfree.h>


  5: /* TRON Routines */
  6: static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*);
  7: /*------------------------------------------------------------*/
  8: static PetscErrorCode TaoDestroy_TRON(Tao tao)
  9: {
 10:   TAO_TRON       *tron = (TAO_TRON *)tao->data;

 14:   VecDestroy(&tron->X_New);
 15:   VecDestroy(&tron->G_New);
 16:   VecDestroy(&tron->Work);
 17:   VecDestroy(&tron->DXFree);
 18:   VecDestroy(&tron->R);
 19:   VecDestroy(&tron->diag);
 20:   VecScatterDestroy(&tron->scatter);
 21:   ISDestroy(&tron->Free_Local);
 22:   MatDestroy(&tron->H_sub);
 23:   MatDestroy(&tron->Hpre_sub);
 24:   PetscFree(tao->data);
 25:   return(0);
 26: }

 28: /*------------------------------------------------------------*/
 29: static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao)
 30: {
 31:   TAO_TRON       *tron = (TAO_TRON *)tao->data;
 33:   PetscBool      flg;

 36:   PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");
 37:   PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
 38:   PetscOptionsTail();
 39:   TaoLineSearchSetFromOptions(tao->linesearch);
 40:   KSPSetFromOptions(tao->ksp);
 41:   return(0);
 42: }

 44: /*------------------------------------------------------------*/
 45: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
 46: {
 47:   TAO_TRON         *tron = (TAO_TRON *)tao->data;
 48:   PetscBool        isascii;
 49:   PetscErrorCode   ierr;

 52:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 53:   if (isascii) {
 54:     PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
 55:     PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
 56:   }
 57:   return(0);
 58: }

 60: /* ---------------------------------------------------------- */
 61: static PetscErrorCode TaoSetup_TRON(Tao tao)
 62: {
 64:   TAO_TRON       *tron = (TAO_TRON *)tao->data;


 68:   /* Allocate some arrays */
 69:   VecDuplicate(tao->solution, &tron->diag);
 70:   VecDuplicate(tao->solution, &tron->X_New);
 71:   VecDuplicate(tao->solution, &tron->G_New);
 72:   VecDuplicate(tao->solution, &tron->Work);
 73:   VecDuplicate(tao->solution, &tao->gradient);
 74:   VecDuplicate(tao->solution, &tao->stepdirection);
 75:   if (!tao->XL) {
 76:     VecDuplicate(tao->solution, &tao->XL);
 77:     VecSet(tao->XL, PETSC_NINFINITY);
 78:   }
 79:   if (!tao->XU) {
 80:     VecDuplicate(tao->solution, &tao->XU);
 81:     VecSet(tao->XU, PETSC_INFINITY);
 82:   }
 83:   return(0);
 84: }

 86: static PetscErrorCode TaoSolve_TRON(Tao tao)
 87: {
 88:   TAO_TRON                     *tron = (TAO_TRON *)tao->data;
 89:   PetscErrorCode               ierr;
 90:   PetscInt                     its;
 91:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
 92:   PetscReal                    prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;

 95:   tron->pgstepsize = 1.0;
 96:   tao->trust = tao->trust0;
 97:   /*   Project the current point onto the feasible set */
 98:   TaoComputeVariableBounds(tao);
 99:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);

101:   /* Project the initial point onto the feasible region */
102:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);

104:   /* Compute the objective function and gradient */
105:   TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
106:   VecNorm(tao->gradient,NORM_2,&tron->gnorm);
107:   if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

109:   /* Project the gradient and calculate the norm */
110:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
111:   VecNorm(tao->gradient,NORM_2,&tron->gnorm);

113:   /* Initialize trust region radius */
114:   tao->trust=tao->trust0;
115:   if (tao->trust <= 0) {
116:     tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
117:   }

119:   /* Initialize step sizes for the line searches */
120:   tron->pgstepsize=1.0;
121:   tron->stepsize=tao->trust;

123:   tao->reason = TAO_CONTINUE_ITERATING;
124:   TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
125:   TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);
126:   (*tao->ops->convergencetest)(tao,tao->cnvP);
127:   while (tao->reason==TAO_CONTINUE_ITERATING){
128:     /* Call general purpose update function */
129:     if (tao->ops->update) {
130:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
131:     }
132:     
133:     /* Perform projected gradient iterations */
134:     TronGradientProjections(tao,tron);

136:     VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
137:     VecNorm(tao->gradient,NORM_2,&tron->gnorm);

139:     tao->ksp_its=0;
140:     f=tron->f; delta=tao->trust;
141:     tron->n_free_last = tron->n_free;
142:     TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);

144:     /* Generate index set (IS) of which bound constraints are active */
145:     ISDestroy(&tron->Free_Local);
146:     VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
147:     ISGetSize(tron->Free_Local, &tron->n_free);

149:     /* If no free variables */
150:     if (tron->n_free == 0) {
151:       VecNorm(tao->gradient,NORM_2,&tron->gnorm);
152:       TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
153:       TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);
154:       (*tao->ops->convergencetest)(tao,tao->cnvP);
155:       if (!tao->reason) {
156:         tao->reason = TAO_CONVERGED_STEPTOL;
157:       }
158:       break;
159:     }
160:     /* use free_local to mask/submat gradient, hessian, stepdirection */
161:     TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
162:     TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
163:     VecSet(tron->DXFree,0.0);
164:     VecScale(tron->R, -1.0);
165:     TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
166:     if (tao->hessian == tao->hessian_pre) {
167:       MatDestroy(&tron->Hpre_sub);
168:       PetscObjectReference((PetscObject)(tron->H_sub));
169:       tron->Hpre_sub = tron->H_sub;
170:     } else {
171:       TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
172:     }
173:     KSPReset(tao->ksp);
174:     KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
175:     while (1) {

177:       /* Approximately solve the reduced linear system */
178:       KSPCGSetRadius(tao->ksp,delta);

180:       KSPSolve(tao->ksp, tron->R, tron->DXFree);
181:       KSPGetIterationNumber(tao->ksp,&its);
182:       tao->ksp_its+=its;
183:       tao->ksp_tot_its+=its;
184:       VecSet(tao->stepdirection,0.0);

186:       /* Add dxfree matrix to compute step direction vector */
187:       VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);

189:       VecDot(tao->gradient, tao->stepdirection, &gdx);
190:       VecCopy(tao->solution, tron->X_New);
191:       VecCopy(tao->gradient, tron->G_New);

193:       stepsize=1.0;f_new=f;

195:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
196:       TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
197:       TaoAddLineSearchCounts(tao);

199:       MatMult(tao->hessian, tao->stepdirection, tron->Work);
200:       VecAYPX(tron->Work, 0.5, tao->gradient);
201:       VecDot(tao->stepdirection, tron->Work, &prered);
202:       actred = f_new - f;
203:       if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
204:         rhok = 1.0;
205:       } else if (actred<0) {
206:         rhok=PetscAbs(-actred/prered);
207:       } else {
208:         rhok=0.0;
209:       }

211:       /* Compare actual improvement to the quadratic model */
212:       if (rhok > tron->eta1) { /* Accept the point */
213:         /* d = x_new - x */
214:         VecCopy(tron->X_New, tao->stepdirection);
215:         VecAXPY(tao->stepdirection, -1.0, tao->solution);

217:         VecNorm(tao->stepdirection, NORM_2, &xdiff);
218:         xdiff *= stepsize;

220:         /* Adjust trust region size */
221:         if (rhok < tron->eta2 ){
222:           delta = PetscMin(xdiff,delta)*tron->sigma1;
223:         } else if (rhok > tron->eta4 ){
224:           delta= PetscMin(xdiff,delta)*tron->sigma3;
225:         } else if (rhok > tron->eta3 ){
226:           delta=PetscMin(xdiff,delta)*tron->sigma2;
227:         }
228:         VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
229:         ISDestroy(&tron->Free_Local);
230:         VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
231:         f=f_new;
232:         VecNorm(tao->gradient,NORM_2,&tron->gnorm);
233:         VecCopy(tron->X_New, tao->solution);
234:         VecCopy(tron->G_New, tao->gradient);
235:         break;
236:       }
237:       else if (delta <= 1e-30) {
238:         break;
239:       }
240:       else {
241:         delta /= 4.0;
242:       }
243:     } /* end linear solve loop */

245:     tron->f=f; tron->actred=actred; tao->trust=delta;
246:     tao->niter++;
247:     TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
248:     TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);
249:     (*tao->ops->convergencetest)(tao,tao->cnvP);
250:   }  /* END MAIN LOOP  */
251:   return(0);
252: }

254: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
255: {
256:   PetscErrorCode               ierr;
257:   PetscInt                     i;
258:   TaoLineSearchConvergedReason ls_reason;
259:   PetscReal                    actred=-1.0,actred_max=0.0;
260:   PetscReal                    f_new;
261:   /*
262:      The gradient and function value passed into and out of this
263:      routine should be current and correct.

265:      The free, active, and binding variables should be already identified
266:   */

269:   for (i=0;i<tron->maxgpits;++i){

271:     if (-actred <= (tron->pg_ftol)*actred_max) break;

273:     ++tron->gp_iterates;
274:     ++tron->total_gp_its;
275:     f_new=tron->f;

277:     VecCopy(tao->gradient,tao->stepdirection);
278:     VecScale(tao->stepdirection,-1.0);
279:     TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
280:     TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
281:                               &tron->pgstepsize, &ls_reason);
282:     TaoAddLineSearchCounts(tao);

284:     VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
285:     VecNorm(tao->gradient,NORM_2,&tron->gnorm);

287:     /* Update the iterate */
288:     actred = f_new - tron->f;
289:     actred_max = PetscMax(actred_max,-(f_new - tron->f));
290:     tron->f = f_new;
291:   }
292:   return(0);
293: }

295: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) {

297:   TAO_TRON       *tron = (TAO_TRON *)tao->data;

304:   if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");

306:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);
307:   VecCopy(tron->Work,DXL);
308:   VecAXPY(DXL,-1.0,tao->gradient);
309:   VecSet(DXU,0.0);
310:   VecPointwiseMax(DXL,DXL,DXU);

312:   VecCopy(tao->gradient,DXU);
313:   VecAXPY(DXU,-1.0,tron->Work);
314:   VecSet(tron->Work,0.0);
315:   VecPointwiseMin(DXU,tron->Work,DXU);
316:   return(0);
317: }

319: /*------------------------------------------------------------*/
320: /*MC
321:   TAOTRON - The TRON algorithm is an active-set Newton trust region method
322:   for bound-constrained minimization.

324:   Options Database Keys:
325: + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
326: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets

328:   Level: beginner
329: M*/
330: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
331: {
332:   TAO_TRON       *tron;
334:   const char     *morethuente_type = TAOLINESEARCHMT;

337:   tao->ops->setup          = TaoSetup_TRON;
338:   tao->ops->solve          = TaoSolve_TRON;
339:   tao->ops->view           = TaoView_TRON;
340:   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
341:   tao->ops->destroy        = TaoDestroy_TRON;
342:   tao->ops->computedual    = TaoComputeDual_TRON;

344:   PetscNewLog(tao,&tron);
345:   tao->data = (void*)tron;

347:   /* Override default settings (unless already changed) */
348:   if (!tao->max_it_changed) tao->max_it = 50;
349:   if (!tao->trust0_changed) tao->trust0 = 1.0;
350:   if (!tao->steptol_changed) tao->steptol = 0.0;

352:   /* Initialize pointers and variables */
353:   tron->n            = 0;
354:   tron->maxgpits     = 3;
355:   tron->pg_ftol      = 0.001;

357:   tron->eta1         = 1.0e-4;
358:   tron->eta2         = 0.25;
359:   tron->eta3         = 0.50;
360:   tron->eta4         = 0.90;

362:   tron->sigma1       = 0.5;
363:   tron->sigma2       = 2.0;
364:   tron->sigma3       = 4.0;

366:   tron->gp_iterates  = 0; /* Cumulative number */
367:   tron->total_gp_its = 0;
368:   tron->n_free       = 0;

370:   tron->DXFree=NULL;
371:   tron->R=NULL;
372:   tron->X_New=NULL;
373:   tron->G_New=NULL;
374:   tron->Work=NULL;
375:   tron->Free_Local=NULL;
376:   tron->H_sub=NULL;
377:   tron->Hpre_sub=NULL;
378:   tao->subset_type = TAO_SUBSET_SUBVEC;

380:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
381:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
382:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
383:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
384:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

386:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
387:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
388:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
389:   KSPSetType(tao->ksp,KSPSTCG);
390:   return(0);
391: }