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: }