Actual source code: tron.c
1: #include <../src/tao/bound/impls/tron/tron.h>
2: #include <../src/tao/matrix/submatfree.h>
4: /* TRON Routines */
5: static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*);
6: /*------------------------------------------------------------*/
7: static PetscErrorCode TaoDestroy_TRON(Tao tao)
8: {
9: TAO_TRON *tron = (TAO_TRON *)tao->data;
11: VecDestroy(&tron->X_New);
12: VecDestroy(&tron->G_New);
13: VecDestroy(&tron->Work);
14: VecDestroy(&tron->DXFree);
15: VecDestroy(&tron->R);
16: VecDestroy(&tron->diag);
17: VecScatterDestroy(&tron->scatter);
18: ISDestroy(&tron->Free_Local);
19: MatDestroy(&tron->H_sub);
20: MatDestroy(&tron->Hpre_sub);
21: PetscFree(tao->data);
22: return 0;
23: }
25: /*------------------------------------------------------------*/
26: static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao)
27: {
28: TAO_TRON *tron = (TAO_TRON *)tao->data;
29: PetscBool flg;
31: PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");
32: PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
33: PetscOptionsTail();
34: KSPSetFromOptions(tao->ksp);
35: return 0;
36: }
38: /*------------------------------------------------------------*/
39: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
40: {
41: TAO_TRON *tron = (TAO_TRON *)tao->data;
42: PetscBool isascii;
44: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
45: if (isascii) {
46: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
47: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
48: }
49: return 0;
50: }
52: /* ---------------------------------------------------------- */
53: static PetscErrorCode TaoSetup_TRON(Tao tao)
54: {
55: TAO_TRON *tron = (TAO_TRON *)tao->data;
58: /* Allocate some arrays */
59: VecDuplicate(tao->solution, &tron->diag);
60: VecDuplicate(tao->solution, &tron->X_New);
61: VecDuplicate(tao->solution, &tron->G_New);
62: VecDuplicate(tao->solution, &tron->Work);
63: VecDuplicate(tao->solution, &tao->gradient);
64: VecDuplicate(tao->solution, &tao->stepdirection);
65: if (!tao->XL) {
66: VecDuplicate(tao->solution, &tao->XL);
67: VecSet(tao->XL, PETSC_NINFINITY);
68: }
69: if (!tao->XU) {
70: VecDuplicate(tao->solution, &tao->XU);
71: VecSet(tao->XU, PETSC_INFINITY);
72: }
73: return 0;
74: }
76: static PetscErrorCode TaoSolve_TRON(Tao tao)
77: {
78: TAO_TRON *tron = (TAO_TRON *)tao->data;
79: PetscInt its;
80: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
81: PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
83: tron->pgstepsize = 1.0;
84: tao->trust = tao->trust0;
85: /* Project the current point onto the feasible set */
86: TaoComputeVariableBounds(tao);
87: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
89: /* Project the initial point onto the feasible region */
90: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
92: /* Compute the objective function and gradient */
93: TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
94: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
97: /* Project the gradient and calculate the norm */
98: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
99: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
101: /* Initialize trust region radius */
102: tao->trust=tao->trust0;
103: if (tao->trust <= 0) {
104: tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
105: }
107: /* Initialize step sizes for the line searches */
108: tron->pgstepsize=1.0;
109: tron->stepsize=tao->trust;
111: tao->reason = TAO_CONTINUE_ITERATING;
112: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
113: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);
114: (*tao->ops->convergencetest)(tao,tao->cnvP);
115: while (tao->reason==TAO_CONTINUE_ITERATING) {
116: /* Call general purpose update function */
117: if (tao->ops->update) {
118: (*tao->ops->update)(tao, tao->niter, tao->user_update);
119: }
121: /* Perform projected gradient iterations */
122: TronGradientProjections(tao,tron);
124: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
125: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
127: tao->ksp_its=0;
128: f=tron->f; delta=tao->trust;
129: tron->n_free_last = tron->n_free;
130: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
132: /* Generate index set (IS) of which bound constraints are active */
133: ISDestroy(&tron->Free_Local);
134: VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
135: ISGetSize(tron->Free_Local, &tron->n_free);
137: /* If no free variables */
138: if (tron->n_free == 0) {
139: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
140: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
141: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);
142: (*tao->ops->convergencetest)(tao,tao->cnvP);
143: if (!tao->reason) {
144: tao->reason = TAO_CONVERGED_STEPTOL;
145: }
146: break;
147: }
148: /* use free_local to mask/submat gradient, hessian, stepdirection */
149: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
150: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
151: VecSet(tron->DXFree,0.0);
152: VecScale(tron->R, -1.0);
153: TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
154: if (tao->hessian == tao->hessian_pre) {
155: MatDestroy(&tron->Hpre_sub);
156: PetscObjectReference((PetscObject)(tron->H_sub));
157: tron->Hpre_sub = tron->H_sub;
158: } else {
159: TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
160: }
161: KSPReset(tao->ksp);
162: KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
163: while (1) {
165: /* Approximately solve the reduced linear system */
166: KSPCGSetRadius(tao->ksp,delta);
168: KSPSolve(tao->ksp, tron->R, tron->DXFree);
169: KSPGetIterationNumber(tao->ksp,&its);
170: tao->ksp_its+=its;
171: tao->ksp_tot_its+=its;
172: VecSet(tao->stepdirection,0.0);
174: /* Add dxfree matrix to compute step direction vector */
175: VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);
177: VecDot(tao->gradient, tao->stepdirection, &gdx);
178: VecCopy(tao->solution, tron->X_New);
179: VecCopy(tao->gradient, tron->G_New);
181: stepsize=1.0;f_new=f;
183: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
184: TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
185: TaoAddLineSearchCounts(tao);
187: MatMult(tao->hessian, tao->stepdirection, tron->Work);
188: VecAYPX(tron->Work, 0.5, tao->gradient);
189: VecDot(tao->stepdirection, tron->Work, &prered);
190: actred = f_new - f;
191: if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
192: rhok = 1.0;
193: } else if (actred<0) {
194: rhok=PetscAbs(-actred/prered);
195: } else {
196: rhok=0.0;
197: }
199: /* Compare actual improvement to the quadratic model */
200: if (rhok > tron->eta1) { /* Accept the point */
201: /* d = x_new - x */
202: VecCopy(tron->X_New, tao->stepdirection);
203: VecAXPY(tao->stepdirection, -1.0, tao->solution);
205: VecNorm(tao->stepdirection, NORM_2, &xdiff);
206: xdiff *= stepsize;
208: /* Adjust trust region size */
209: if (rhok < tron->eta2) {
210: delta = PetscMin(xdiff,delta)*tron->sigma1;
211: } else if (rhok > tron->eta4) {
212: delta= PetscMin(xdiff,delta)*tron->sigma3;
213: } else if (rhok > tron->eta3) {
214: delta=PetscMin(xdiff,delta)*tron->sigma2;
215: }
216: VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
217: ISDestroy(&tron->Free_Local);
218: VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
219: f=f_new;
220: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
221: VecCopy(tron->X_New, tao->solution);
222: VecCopy(tron->G_New, tao->gradient);
223: break;
224: }
225: else if (delta <= 1e-30) {
226: break;
227: }
228: else {
229: delta /= 4.0;
230: }
231: } /* end linear solve loop */
233: tron->f=f; tron->actred=actred; tao->trust=delta;
234: tao->niter++;
235: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
236: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);
237: (*tao->ops->convergencetest)(tao,tao->cnvP);
238: } /* END MAIN LOOP */
239: return 0;
240: }
242: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
243: {
244: PetscErrorCode ierr;
245: PetscInt i;
246: TaoLineSearchConvergedReason ls_reason;
247: PetscReal actred=-1.0,actred_max=0.0;
248: PetscReal f_new;
249: /*
250: The gradient and function value passed into and out of this
251: routine should be current and correct.
253: The free, active, and binding variables should be already identified
254: */
256: for (i=0;i<tron->maxgpits;++i) {
258: if (-actred <= (tron->pg_ftol)*actred_max) break;
260: ++tron->gp_iterates;
261: ++tron->total_gp_its;
262: f_new=tron->f;
264: VecCopy(tao->gradient,tao->stepdirection);
265: VecScale(tao->stepdirection,-1.0);
266: TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
267: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
268: &tron->pgstepsize, &ls_reason);
269: TaoAddLineSearchCounts(tao);
271: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
272: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
274: /* Update the iterate */
275: actred = f_new - tron->f;
276: actred_max = PetscMax(actred_max,-(f_new - tron->f));
277: tron->f = f_new;
278: }
279: return 0;
280: }
282: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU)
283: {
285: TAO_TRON *tron = (TAO_TRON *)tao->data;
292: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);
293: VecCopy(tron->Work,DXL);
294: VecAXPY(DXL,-1.0,tao->gradient);
295: VecSet(DXU,0.0);
296: VecPointwiseMax(DXL,DXL,DXU);
298: VecCopy(tao->gradient,DXU);
299: VecAXPY(DXU,-1.0,tron->Work);
300: VecSet(tron->Work,0.0);
301: VecPointwiseMin(DXU,tron->Work,DXU);
302: return 0;
303: }
305: /*------------------------------------------------------------*/
306: /*MC
307: TAOTRON - The TRON algorithm is an active-set Newton trust region method
308: for bound-constrained minimization.
310: Options Database Keys:
311: + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
312: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
314: Level: beginner
315: M*/
316: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
317: {
318: TAO_TRON *tron;
319: const char *morethuente_type = TAOLINESEARCHMT;
321: tao->ops->setup = TaoSetup_TRON;
322: tao->ops->solve = TaoSolve_TRON;
323: tao->ops->view = TaoView_TRON;
324: tao->ops->setfromoptions = TaoSetFromOptions_TRON;
325: tao->ops->destroy = TaoDestroy_TRON;
326: tao->ops->computedual = TaoComputeDual_TRON;
328: PetscNewLog(tao,&tron);
329: tao->data = (void*)tron;
331: /* Override default settings (unless already changed) */
332: if (!tao->max_it_changed) tao->max_it = 50;
333: if (!tao->trust0_changed) tao->trust0 = 1.0;
334: if (!tao->steptol_changed) tao->steptol = 0.0;
336: /* Initialize pointers and variables */
337: tron->n = 0;
338: tron->maxgpits = 3;
339: tron->pg_ftol = 0.001;
341: tron->eta1 = 1.0e-4;
342: tron->eta2 = 0.25;
343: tron->eta3 = 0.50;
344: tron->eta4 = 0.90;
346: tron->sigma1 = 0.5;
347: tron->sigma2 = 2.0;
348: tron->sigma3 = 4.0;
350: tron->gp_iterates = 0; /* Cumulative number */
351: tron->total_gp_its = 0;
352: tron->n_free = 0;
354: tron->DXFree=NULL;
355: tron->R=NULL;
356: tron->X_New=NULL;
357: tron->G_New=NULL;
358: tron->Work=NULL;
359: tron->Free_Local=NULL;
360: tron->H_sub=NULL;
361: tron->Hpre_sub=NULL;
362: tao->subset_type = TAO_SUBSET_SUBVEC;
364: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
365: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
366: TaoLineSearchSetType(tao->linesearch,morethuente_type);
367: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
368: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
370: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
371: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
372: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
373: KSPSetType(tao->ksp,KSPSTCG);
374: return 0;
375: }