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;
13: VecDestroy(&tron->X_New);
14: VecDestroy(&tron->G_New);
15: VecDestroy(&tron->Work);
16: VecDestroy(&tron->DXFree);
17: VecDestroy(&tron->R);
18: VecDestroy(&tron->diag);
19: VecScatterDestroy(&tron->scatter);
20: ISDestroy(&tron->Free_Local);
21: MatDestroy(&tron->H_sub);
22: MatDestroy(&tron->Hpre_sub);
23: PetscFree(tao->data);
24: return(0);
25: }
27: /*------------------------------------------------------------*/
28: static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao)
29: {
30: TAO_TRON *tron = (TAO_TRON *)tao->data;
32: PetscBool flg;
35: PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");
36: PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
37: PetscOptionsTail();
38: TaoLineSearchSetFromOptions(tao->linesearch);
39: KSPSetFromOptions(tao->ksp);
40: return(0);
41: }
43: /*------------------------------------------------------------*/
44: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
45: {
46: TAO_TRON *tron = (TAO_TRON *)tao->data;
47: PetscBool isascii;
48: PetscErrorCode ierr;
51: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
52: if (isascii) {
53: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
54: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
55: }
56: return(0);
57: }
59: /* ---------------------------------------------------------- */
60: static PetscErrorCode TaoSetup_TRON(Tao tao)
61: {
63: TAO_TRON *tron = (TAO_TRON *)tao->data;
67: /* Allocate some arrays */
68: VecDuplicate(tao->solution, &tron->diag);
69: VecDuplicate(tao->solution, &tron->X_New);
70: VecDuplicate(tao->solution, &tron->G_New);
71: VecDuplicate(tao->solution, &tron->Work);
72: VecDuplicate(tao->solution, &tao->gradient);
73: VecDuplicate(tao->solution, &tao->stepdirection);
74: if (!tao->XL) {
75: VecDuplicate(tao->solution, &tao->XL);
76: VecSet(tao->XL, PETSC_NINFINITY);
77: }
78: if (!tao->XU) {
79: VecDuplicate(tao->solution, &tao->XU);
80: VecSet(tao->XU, PETSC_INFINITY);
81: }
82: return(0);
83: }
85: static PetscErrorCode TaoSolve_TRON(Tao tao)
86: {
87: TAO_TRON *tron = (TAO_TRON *)tao->data;
88: PetscErrorCode ierr;
89: PetscInt its;
90: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
91: PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
94: tron->pgstepsize = 1.0;
95: tao->trust = tao->trust0;
96: /* Project the current point onto the feasible set */
97: TaoComputeVariableBounds(tao);
98: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
100: /* Project the initial point onto the feasible region */
101: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
103: /* Compute the objective function and gradient */
104: TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
105: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
106: if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
108: /* Project the gradient and calculate the norm */
109: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
110: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
112: /* Initialize trust region radius */
113: tao->trust=tao->trust0;
114: if (tao->trust <= 0) {
115: tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
116: }
118: /* Initialize step sizes for the line searches */
119: tron->pgstepsize=1.0;
120: tron->stepsize=tao->trust;
122: tao->reason = TAO_CONTINUE_ITERATING;
123: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
124: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);
125: (*tao->ops->convergencetest)(tao,tao->cnvP);
126: while (tao->reason==TAO_CONTINUE_ITERATING) {
127: /* Call general purpose update function */
128: if (tao->ops->update) {
129: (*tao->ops->update)(tao, tao->niter, tao->user_update);
130: }
132: /* Perform projected gradient iterations */
133: TronGradientProjections(tao,tron);
135: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
136: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
138: tao->ksp_its=0;
139: f=tron->f; delta=tao->trust;
140: tron->n_free_last = tron->n_free;
141: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
143: /* Generate index set (IS) of which bound constraints are active */
144: ISDestroy(&tron->Free_Local);
145: VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
146: ISGetSize(tron->Free_Local, &tron->n_free);
148: /* If no free variables */
149: if (tron->n_free == 0) {
150: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
151: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
152: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);
153: (*tao->ops->convergencetest)(tao,tao->cnvP);
154: if (!tao->reason) {
155: tao->reason = TAO_CONVERGED_STEPTOL;
156: }
157: break;
158: }
159: /* use free_local to mask/submat gradient, hessian, stepdirection */
160: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
161: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
162: VecSet(tron->DXFree,0.0);
163: VecScale(tron->R, -1.0);
164: TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
165: if (tao->hessian == tao->hessian_pre) {
166: MatDestroy(&tron->Hpre_sub);
167: PetscObjectReference((PetscObject)(tron->H_sub));
168: tron->Hpre_sub = tron->H_sub;
169: } else {
170: TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
171: }
172: KSPReset(tao->ksp);
173: KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
174: while (1) {
176: /* Approximately solve the reduced linear system */
177: KSPCGSetRadius(tao->ksp,delta);
179: KSPSolve(tao->ksp, tron->R, tron->DXFree);
180: KSPGetIterationNumber(tao->ksp,&its);
181: tao->ksp_its+=its;
182: tao->ksp_tot_its+=its;
183: VecSet(tao->stepdirection,0.0);
185: /* Add dxfree matrix to compute step direction vector */
186: VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);
188: VecDot(tao->gradient, tao->stepdirection, &gdx);
189: VecCopy(tao->solution, tron->X_New);
190: VecCopy(tao->gradient, tron->G_New);
192: stepsize=1.0;f_new=f;
194: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
195: TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
196: TaoAddLineSearchCounts(tao);
198: MatMult(tao->hessian, tao->stepdirection, tron->Work);
199: VecAYPX(tron->Work, 0.5, tao->gradient);
200: VecDot(tao->stepdirection, tron->Work, &prered);
201: actred = f_new - f;
202: if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
203: rhok = 1.0;
204: } else if (actred<0) {
205: rhok=PetscAbs(-actred/prered);
206: } else {
207: rhok=0.0;
208: }
210: /* Compare actual improvement to the quadratic model */
211: if (rhok > tron->eta1) { /* Accept the point */
212: /* d = x_new - x */
213: VecCopy(tron->X_New, tao->stepdirection);
214: VecAXPY(tao->stepdirection, -1.0, tao->solution);
216: VecNorm(tao->stepdirection, NORM_2, &xdiff);
217: xdiff *= stepsize;
219: /* Adjust trust region size */
220: if (rhok < tron->eta2) {
221: delta = PetscMin(xdiff,delta)*tron->sigma1;
222: } else if (rhok > tron->eta4) {
223: delta= PetscMin(xdiff,delta)*tron->sigma3;
224: } else if (rhok > tron->eta3) {
225: delta=PetscMin(xdiff,delta)*tron->sigma2;
226: }
227: VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
228: ISDestroy(&tron->Free_Local);
229: VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
230: f=f_new;
231: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
232: VecCopy(tron->X_New, tao->solution);
233: VecCopy(tron->G_New, tao->gradient);
234: break;
235: }
236: else if (delta <= 1e-30) {
237: break;
238: }
239: else {
240: delta /= 4.0;
241: }
242: } /* end linear solve loop */
244: tron->f=f; tron->actred=actred; tao->trust=delta;
245: tao->niter++;
246: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
247: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);
248: (*tao->ops->convergencetest)(tao,tao->cnvP);
249: } /* END MAIN LOOP */
250: return(0);
251: }
253: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
254: {
255: PetscErrorCode ierr;
256: PetscInt i;
257: TaoLineSearchConvergedReason ls_reason;
258: PetscReal actred=-1.0,actred_max=0.0;
259: PetscReal f_new;
260: /*
261: The gradient and function value passed into and out of this
262: routine should be current and correct.
264: The free, active, and binding variables should be already identified
265: */
268: for (i=0;i<tron->maxgpits;++i) {
270: if (-actred <= (tron->pg_ftol)*actred_max) break;
272: ++tron->gp_iterates;
273: ++tron->total_gp_its;
274: f_new=tron->f;
276: VecCopy(tao->gradient,tao->stepdirection);
277: VecScale(tao->stepdirection,-1.0);
278: TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
279: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
280: &tron->pgstepsize, &ls_reason);
281: TaoAddLineSearchCounts(tao);
283: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
284: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
286: /* Update the iterate */
287: actred = f_new - tron->f;
288: actred_max = PetscMax(actred_max,-(f_new - tron->f));
289: tron->f = f_new;
290: }
291: return(0);
292: }
294: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU)
295: {
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: }