Actual source code: tron.c
petsc-3.8.4 2018-03-24
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: PetscViewerASCIIPushTab(viewer);
55: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
56: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
57: PetscViewerASCIIPopTab(viewer);
58: }
59: return(0);
60: }
63: /* ---------------------------------------------------------- */
64: static PetscErrorCode TaoSetup_TRON(Tao tao)
65: {
67: TAO_TRON *tron = (TAO_TRON *)tao->data;
71: /* Allocate some arrays */
72: VecDuplicate(tao->solution, &tron->diag);
73: VecDuplicate(tao->solution, &tron->X_New);
74: VecDuplicate(tao->solution, &tron->G_New);
75: VecDuplicate(tao->solution, &tron->Work);
76: VecDuplicate(tao->solution, &tao->gradient);
77: VecDuplicate(tao->solution, &tao->stepdirection);
78: if (!tao->XL) {
79: VecDuplicate(tao->solution, &tao->XL);
80: VecSet(tao->XL, PETSC_NINFINITY);
81: }
82: if (!tao->XU) {
83: VecDuplicate(tao->solution, &tao->XU);
84: VecSet(tao->XU, PETSC_INFINITY);
85: }
86: return(0);
87: }
91: static PetscErrorCode TaoSolve_TRON(Tao tao)
92: {
93: TAO_TRON *tron = (TAO_TRON *)tao->data;
94: PetscErrorCode ierr;
95: PetscInt its;
96: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
97: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
98: PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
101: tron->pgstepsize=1.0;
102: tao->trust = tao->trust0;
103: /* Project the current point onto the feasible set */
104: TaoComputeVariableBounds(tao);
105: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
106: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
108: TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
109: ISDestroy(&tron->Free_Local);
111: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
113: /* Project the gradient and calculate the norm */
114: VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);
115: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
117: if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
118: if (tao->trust <= 0) {
119: tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
120: }
122: tron->stepsize=tao->trust;
123: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);
124: while (reason==TAO_CONTINUE_ITERATING){
125: tao->ksp_its=0;
126: TronGradientProjections(tao,tron);
127: f=tron->f; delta=tao->trust;
128: tron->n_free_last = tron->n_free;
129: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
131: ISGetSize(tron->Free_Local, &tron->n_free);
133: /* If no free variables */
134: if (tron->n_free == 0) {
135: PetscInfo(tao,"No free variables in tron iteration.\n");
136: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
137: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);
138: if (!reason) {
139: reason = TAO_CONVERGED_STEPTOL;
140: TaoSetConvergedReason(tao,reason);
141: }
142: break;
143: }
144: /* use free_local to mask/submat gradient, hessian, stepdirection */
145: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
146: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
147: VecSet(tron->DXFree,0.0);
148: VecScale(tron->R, -1.0);
149: TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
150: if (tao->hessian == tao->hessian_pre) {
151: MatDestroy(&tron->Hpre_sub);
152: PetscObjectReference((PetscObject)(tron->H_sub));
153: tron->Hpre_sub = tron->H_sub;
154: } else {
155: TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
156: }
157: KSPReset(tao->ksp);
158: KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
159: while (1) {
161: /* Approximately solve the reduced linear system */
162: KSPCGSetRadius(tao->ksp,delta);
164: KSPSolve(tao->ksp, tron->R, tron->DXFree);
165: KSPGetIterationNumber(tao->ksp,&its);
166: tao->ksp_its+=its;
167: tao->ksp_tot_its+=its;
168: VecSet(tao->stepdirection,0.0);
170: /* Add dxfree matrix to compute step direction vector */
171: VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);
173: VecDot(tao->gradient, tao->stepdirection, &gdx);
174: PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",(double)gdx);
176: VecCopy(tao->solution, tron->X_New);
177: VecCopy(tao->gradient, tron->G_New);
179: stepsize=1.0;f_new=f;
181: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
182: TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
183: TaoAddLineSearchCounts(tao);
185: MatMult(tao->hessian, tao->stepdirection, tron->Work);
186: VecAYPX(tron->Work, 0.5, tao->gradient);
187: VecDot(tao->stepdirection, tron->Work, &prered);
188: actred = f_new - f;
189: if (actred<0) {
190: rhok=PetscAbs(-actred/prered);
191: } else {
192: rhok=0.0;
193: }
195: /* Compare actual improvement to the quadratic model */
196: if (rhok > tron->eta1) { /* Accept the point */
197: /* d = x_new - x */
198: VecCopy(tron->X_New, tao->stepdirection);
199: VecAXPY(tao->stepdirection, -1.0, tao->solution);
201: VecNorm(tao->stepdirection, NORM_2, &xdiff);
202: xdiff *= stepsize;
204: /* Adjust trust region size */
205: if (rhok < tron->eta2 ){
206: delta = PetscMin(xdiff,delta)*tron->sigma1;
207: } else if (rhok > tron->eta4 ){
208: delta= PetscMin(xdiff,delta)*tron->sigma3;
209: } else if (rhok > tron->eta3 ){
210: delta=PetscMin(xdiff,delta)*tron->sigma2;
211: }
212: VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
213: ISDestroy(&tron->Free_Local);
214: VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);
215: f=f_new;
216: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
217: VecCopy(tron->X_New, tao->solution);
218: VecCopy(tron->G_New, tao->gradient);
219: break;
220: }
221: else if (delta <= 1e-30) {
222: break;
223: }
224: else {
225: delta /= 4.0;
226: }
227: } /* end linear solve loop */
229: tron->f=f; tron->actred=actred; tao->trust=delta;
230: tao->niter++;
231: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);
232: } /* END MAIN LOOP */
234: return(0);
235: }
238: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
239: {
240: PetscErrorCode ierr;
241: PetscInt i;
242: TaoLineSearchConvergedReason ls_reason;
243: PetscReal actred=-1.0,actred_max=0.0;
244: PetscReal f_new;
245: /*
246: The gradient and function value passed into and out of this
247: routine should be current and correct.
249: The free, active, and binding variables should be already identified
250: */
252: ISDestroy(&tron->Free_Local);
253: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
255: for (i=0;i<tron->maxgpits;i++){
257: if ( -actred <= (tron->pg_ftol)*actred_max) break;
259: tron->gp_iterates++; tron->total_gp_its++;
260: f_new=tron->f;
262: VecCopy(tao->gradient, tao->stepdirection);
263: VecScale(tao->stepdirection, -1.0);
264: TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
265: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
266: &tron->pgstepsize, &ls_reason);
267: TaoAddLineSearchCounts(tao);
270: /* Update the iterate */
271: actred = f_new - tron->f;
272: actred_max = PetscMax(actred_max,-(f_new - tron->f));
273: tron->f = f_new;
274: ISDestroy(&tron->Free_Local);
275: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
276: }
278: return(0);
279: }
281: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) {
283: TAO_TRON *tron = (TAO_TRON *)tao->data;
290: if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
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;
320: const char *morethuente_type = TAOLINESEARCHMT;
323: tao->ops->setup = TaoSetup_TRON;
324: tao->ops->solve = TaoSolve_TRON;
325: tao->ops->view = TaoView_TRON;
326: tao->ops->setfromoptions = TaoSetFromOptions_TRON;
327: tao->ops->destroy = TaoDestroy_TRON;
328: tao->ops->computedual = TaoComputeDual_TRON;
330: PetscNewLog(tao,&tron);
331: tao->data = (void*)tron;
333: /* Override default settings (unless already changed) */
334: if (!tao->max_it_changed) tao->max_it = 50;
335: if (!tao->trust0_changed) tao->trust0 = 1.0;
336: if (!tao->steptol_changed) tao->steptol = 0.0;
338: /* Initialize pointers and variables */
339: tron->n = 0;
340: tron->maxgpits = 3;
341: tron->pg_ftol = 0.001;
343: tron->eta1 = 1.0e-4;
344: tron->eta2 = 0.25;
345: tron->eta3 = 0.50;
346: tron->eta4 = 0.90;
348: tron->sigma1 = 0.5;
349: tron->sigma2 = 2.0;
350: tron->sigma3 = 4.0;
352: tron->gp_iterates = 0; /* Cumulative number */
353: tron->total_gp_its = 0;
354: tron->n_free = 0;
356: tron->DXFree=NULL;
357: tron->R=NULL;
358: tron->X_New=NULL;
359: tron->G_New=NULL;
360: tron->Work=NULL;
361: tron->Free_Local=NULL;
362: tron->H_sub=NULL;
363: tron->Hpre_sub=NULL;
364: tao->subset_type = TAO_SUBSET_SUBVEC;
366: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
367: TaoLineSearchSetType(tao->linesearch,morethuente_type);
368: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
369: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
371: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
372: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
373: KSPSetType(tao->ksp,KSPCGSTCG);
374: return(0);
375: }