Actual source code: tron.c

  1: #include "tron.h"
  2: #include "private/kspimpl.h"
  3: #include "private/matimpl.h"
  4: #include "src/matrix/submatfree.h"


  7: /* TRON Routines */
  8: static PetscErrorCode TronGradientProjections(TaoSolver,TAO_TRON*);

 10: /*------------------------------------------------------------*/
 13: static PetscErrorCode TaoDestroy_TRON(TaoSolver tao)
 14: {
 15:   TAO_TRON *tron = (TAO_TRON *)tao->data;


 20:   VecDestroy(&tron->X_New);
 21:   VecDestroy(&tron->G_New);
 22:   VecDestroy(&tron->Work);
 23:   VecDestroy(&tron->DXFree);
 24:   VecDestroy(&tron->R); 
 25:   VecDestroy(&tron->diag); 
 26:   VecScatterDestroy(&tron->scatter); 
 27:   ISDestroy(&tron->Free_Local); 
 28:   MatDestroy(&tron->H_sub); 
 29:   MatDestroy(&tron->Hpre_sub); 
 30:   PetscFree(tao->data);
 31:   tao->data = PETSC_NULL;

 33:   return(0);
 34: }

 36: /*------------------------------------------------------------*/
 39: static PetscErrorCode TaoSetFromOptions_TRON(TaoSolver tao)
 40: {
 41:   TAO_TRON  *tron = (TAO_TRON *)tao->data;
 42:   PetscErrorCode        ierr;
 43:   PetscBool flg;


 47:   PetscOptionsHead("Newton Trust Region Method for bound constrained optimization");
 48:   
 49:   PetscOptionsInt("-tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
 50:   

 52:   PetscOptionsTail();
 53:   TaoLineSearchSetFromOptions(tao->linesearch);
 54:   KSPSetFromOptions(tao->ksp); 

 56:   return(0);
 57: }

 59: /*------------------------------------------------------------*/
 62: static PetscErrorCode TaoView_TRON(TaoSolver tao, PetscViewer viewer)
 63: {
 64:   TAO_TRON  *tron = (TAO_TRON *)tao->data;
 65:   PetscBool isascii;
 66:   PetscErrorCode   ierr;

 69:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 70:   if (isascii) {
 71:     PetscViewerASCIIPushTab(viewer); 
 72:     PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
 73:     PetscViewerASCIIPrintf(viewer,"PG tolerance: %G \n",tron->pg_ftol);
 74:     PetscViewerASCIIPopTab(viewer); 
 75:   } else {
 76:     SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO TRON",((PetscObject)viewer)->type_name);
 77:   }    
 78:   return(0);
 79: }


 82: /* ---------------------------------------------------------- */
 85: static PetscErrorCode TaoSetup_TRON(TaoSolver tao)
 86: {
 88:   TAO_TRON *tron = (TAO_TRON *)tao->data;


 92:   /* Allocate some arrays */
 93:   VecDuplicate(tao->solution, &tron->diag); 
 94:   VecDuplicate(tao->solution, &tron->X_New); 
 95:   VecDuplicate(tao->solution, &tron->G_New); 
 96:   VecDuplicate(tao->solution, &tron->Work); 
 97:   VecDuplicate(tao->solution, &tao->gradient); 
 98:   VecDuplicate(tao->solution, &tao->stepdirection); 
 99:   if (!tao->XL) {
100:       VecDuplicate(tao->solution, &tao->XL); 
101:       VecSet(tao->XL, TAO_NINFINITY); 
102:   }
103:   if (!tao->XU) {
104:       VecDuplicate(tao->solution, &tao->XU); 
105:       VecSet(tao->XU, TAO_INFINITY); 
106:   }

108:   return(0);
109: }



115: static PetscErrorCode TaoSolve_TRON(TaoSolver tao){

117:   TAO_TRON *tron = (TAO_TRON *)tao->data;
119:   PetscInt iter=0,its;

121:   TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
122:   TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
123:   PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;

126:   tron->pgstepsize=1.0;


129:   /*   Project the current point onto the feasible set */
130:   TaoComputeVariableBounds(tao); 
131:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); 
132:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU); 

134:   
135:   TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
136:   if (tron->Free_Local) {
137:     ISDestroy(&tron->Free_Local); 
138:     tron->Free_Local = PETSC_NULL;
139:   }
140:   VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); 
141:   
142:   /* Project the gradient and calculate the norm */
143:   VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient); 
144:   VecNorm(tao->gradient,NORM_2,&tron->gnorm); 

146:   if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) {
147:     SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
148:   }

150:   if (tao->trust <= 0) {
151:     tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
152:   }

154:   tron->stepsize=tao->trust;
155:   TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason); 
156:   while (reason==TAO_CONTINUE_ITERATING){

158:     TronGradientProjections(tao,tron); 
159:     f=tron->f; delta=tao->trust; 
160:     
161:     tron->n_free_last = tron->n_free;
162:     TaoComputeHessian(tao,tao->solution,&tao->hessian, &tao->hessian_pre, &tron->matflag);

164:     ISGetSize(tron->Free_Local, &tron->n_free);  

166:     /* If no free variables */
167:     if (tron->n_free == 0) {
168:       actred=0;
169:       PetscInfo(tao,"No free variables in tron iteration.");
170:       break;

172:     }
173:     /* use free_local to mask/submat gradient, hessian, stepdirection */
174:     VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
175:     
176:     VecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
177:     
178:     VecSet(tron->DXFree,0.0); 
179:     VecScale(tron->R, -1.0); 
180:     MatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub); 
181:     if (tao->hessian == tao->hessian_pre) {
182:       MatDestroy(&tron->Hpre_sub); 
183:       PetscObjectReference((PetscObject)(tron->H_sub)); 
184:       tron->Hpre_sub = tron->H_sub;
185:     } else {
186:       MatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
187:       
188:     }
189:     KSPReset(tao->ksp); 
190:     KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub, tron->matflag); 
191:     while (1) {

193:       /* Approximately solve the reduced linear system */
194:       KSPSolve(tao->ksp, tron->R, tron->DXFree); 
195:       KSPGetIterationNumber(tao->ksp,&its); 
196:       tao->ksp_its+=its;
197:       VecSet(tao->stepdirection,0.0); 
198:       
199:       /* Add dxfree matrix to compute step direction vector */
200:       VecReducedXPY(tao->stepdirection,tron->DXFree,tron->Free_Local); 
201:       if (0) { 
202:         PetscReal rhs,stepnorm;
203:         VecNorm(tron->R,NORM_2,&rhs); 
204:         VecNorm(tron->DXFree,NORM_2,&stepnorm); 
205:         PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",rhs,stepnorm); 
206:       }
207:       
208:       
209:       VecDot(tao->gradient, tao->stepdirection, &gdx); 
210:       PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",gdx); 
211:       
212:       
213:       VecCopy(tao->solution, tron->X_New); 
214:       VecCopy(tao->gradient, tron->G_New); 
215:       
216:       stepsize=1.0;f_new=f;

218:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 
219:       TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,
220:                                 &stepsize,&ls_reason);  
221:       TaoAddLineSearchCounts(tao); 
222:       
223:       MatMult(tao->hessian, tao->stepdirection, tron->Work); 
224:       VecAYPX(tron->Work, 0.5, tao->gradient); 
225:       VecDot(tao->stepdirection, tron->Work, &prered); 
226:       actred = f_new - f;
227:       if (actred<0) {
228:         rhok=PetscAbs(-actred/prered);
229:       } else {
230:         rhok=0.0;
231:       }
232:       
233:       /* Compare actual improvement to the quadratic model */
234:       if (rhok > tron->eta1) { /* Accept the point */
235:         /* d = x_new - x */
236:         VecCopy(tron->X_New, tao->stepdirection); 
237:         VecAXPY(tao->stepdirection, -1.0, tao->solution); 
238:         
239:         VecNorm(tao->stepdirection, NORM_2, &xdiff); 
240:         xdiff *= stepsize;

242:         /* Adjust trust region size */
243:         if (rhok < tron->eta2 ){
244:           delta = PetscMin(xdiff,delta)*tron->sigma1;
245:         } else if (rhok > tron->eta4 ){
246:           delta= PetscMin(xdiff,delta)*tron->sigma3;
247:         } else if (rhok > tron->eta3 ){
248:           delta=PetscMin(xdiff,delta)*tron->sigma2;
249:         }
250:         VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient); 
251:         if (tron->Free_Local) {
252:           ISDestroy(&tron->Free_Local); 
253:           tron->Free_Local=PETSC_NULL;
254:         }
255:         VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local); 
256:         f=f_new;
257:         VecNorm(tao->gradient,NORM_2,&tron->gnorm); 
258:         VecCopy(tron->X_New, tao->solution); 
259:         VecCopy(tron->G_New, tao->gradient); 
260:         break;
261:       } 
262:       else if (delta <= 1e-30) {
263:         break;
264:       }
265:       else {
266:         delta /= 4.0;
267:       }
268:     } /* end linear solve loop */


271:     tron->f=f; tron->actred=actred; tao->trust=delta;
272:     iter++;
273:     TaoMonitor(tao, iter, tron->f, tron->gnorm, 0.0, delta, &reason); 
274:   }  /* END MAIN LOOP  */

276:   return(0);
277: }


282: static PetscErrorCode TronGradientProjections(TaoSolver tao,TAO_TRON *tron)
283: {
285:   PetscInt i;
286:   TaoLineSearchTerminationReason ls_reason;
287:   PetscReal actred=-1.0,actred_max=0.0;
288:   PetscReal f_new;
289:   /*
290:      The gradient and function value passed into and out of this
291:      routine should be current and correct.
292:      
293:      The free, active, and binding variables should be already identified
294:   */
295:   
297:   if (tron->Free_Local) {
298:     ISDestroy(&tron->Free_Local); 
299:     tron->Free_Local = PETSC_NULL;
300:   }
301:   VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); 

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

305:     if ( -actred <= (tron->pg_ftol)*actred_max) break;
306:   
307:     tron->gp_iterates++; tron->total_gp_its++;      
308:     f_new=tron->f;

310:     VecCopy(tao->gradient, tao->stepdirection); 
311:     VecScale(tao->stepdirection, -1.0); 
312:     TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize); 
313:     TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
314:                               &tron->pgstepsize, &ls_reason); 
315:     TaoAddLineSearchCounts(tao); 


318:     /* Update the iterate */
319:     actred = f_new - tron->f;
320:     actred_max = PetscMax(actred_max,-(f_new - tron->f));
321:     tron->f = f_new;
322:     if (tron->Free_Local) {
323:       ISDestroy(&tron->Free_Local); 
324:       tron->Free_Local = PETSC_NULL;
325:     }
326:     VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local); 
327:   }
328:   
329:   return(0);
330: }


335: static PetscErrorCode TaoComputeDual_TRON(TaoSolver tao, Vec DXL, Vec DXU) {

337:   TAO_TRON *tron = (TAO_TRON *)tao->data;
338:   PetscErrorCode       ierr;



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

350:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work); 
351:   VecCopy(tron->Work,DXL); 
352:   VecAXPY(DXL,-1.0,tao->gradient); 
353:   VecSet(DXU,0.0); 
354:   VecPointwiseMax(DXL,DXL,DXU); 

356:   VecCopy(tao->gradient,DXU); 
357:   VecAXPY(DXU,-1.0,tron->Work); 
358:   VecSet(tron->Work,0.0); 
359:   VecPointwiseMin(DXU,tron->Work,DXU); 

361:   return(0);
362: }

364: /*------------------------------------------------------------*/
368: PetscErrorCode TaoCreate_TRON(TaoSolver tao)
369: {
370:   TAO_TRON *tron;
371:   PetscErrorCode   ierr;
372:   const char *morethuente_type = TAOLINESEARCH_MT;

375:   tao->ops->setup = TaoSetup_TRON;
376:   tao->ops->solve = TaoSolve_TRON;
377:   tao->ops->view = TaoView_TRON;
378:   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
379:   tao->ops->destroy = TaoDestroy_TRON;
380:   tao->ops->computedual = TaoComputeDual_TRON;

382:   PetscNewLog(tao,TAO_TRON,&tron); 

384:   tao->max_it = 50;
385:   tao->fatol = 1e-10;
386:   tao->frtol = 1e-10;
387:   tao->data = (void*)tron;
388:   tao->steptol = 1e-12;
389:   tao->trust        = 1.0;

391:   /* Initialize pointers and variables */
392:   tron->n            = 0;
393:   tron->maxgpits     = 3;
394:   tron->pg_ftol      = 0.001;

396:   tron->eta1         = 1.0e-4;
397:   tron->eta2         = 0.25;
398:   tron->eta3         = 0.50;
399:   tron->eta4         = 0.90;

401:   tron->sigma1       = 0.5;
402:   tron->sigma2       = 2.0;
403:   tron->sigma3       = 4.0;

405:   tron->gp_iterates  = 0; /* Cumulative number */
406:   tron->total_gp_its = 0;
407:  
408:   tron->n_free       = 0;

410:   tron->DXFree=PETSC_NULL;
411:   tron->R=PETSC_NULL;
412:   tron->X_New=PETSC_NULL;
413:   tron->G_New=PETSC_NULL;
414:   tron->Work=PETSC_NULL;
415:   tron->Free_Local=PETSC_NULL;
416:   tron->H_sub=PETSC_NULL;
417:   tron->Hpre_sub=PETSC_NULL;
418:   tao->subset_type = TAO_SUBSET_SUBVEC;

420:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
421:   TaoLineSearchSetType(tao->linesearch,morethuente_type); 
422:   TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao); 

424:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp); 


427:   return(0);
428: }