Actual source code: tron.c
petsc-3.10.5 2019-03-28
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: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
55: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
56: }
57: return(0);
58: }
60: /* ---------------------------------------------------------- */
61: static PetscErrorCode TaoSetup_TRON(Tao tao)
62: {
64: TAO_TRON *tron = (TAO_TRON *)tao->data;
68: /* Allocate some arrays */
69: VecDuplicate(tao->solution, &tron->diag);
70: VecDuplicate(tao->solution, &tron->X_New);
71: VecDuplicate(tao->solution, &tron->G_New);
72: VecDuplicate(tao->solution, &tron->Work);
73: VecDuplicate(tao->solution, &tao->gradient);
74: VecDuplicate(tao->solution, &tao->stepdirection);
75: if (!tao->XL) {
76: VecDuplicate(tao->solution, &tao->XL);
77: VecSet(tao->XL, PETSC_NINFINITY);
78: }
79: if (!tao->XU) {
80: VecDuplicate(tao->solution, &tao->XU);
81: VecSet(tao->XU, PETSC_INFINITY);
82: }
83: return(0);
84: }
86: static PetscErrorCode TaoSolve_TRON(Tao tao)
87: {
88: TAO_TRON *tron = (TAO_TRON *)tao->data;
89: PetscErrorCode ierr;
90: PetscInt its;
91: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
92: PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
95: tron->pgstepsize = 1.0;
96: tao->trust = tao->trust0;
97: /* Project the current point onto the feasible set */
98: TaoComputeVariableBounds(tao);
99: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
101: /* Project the initial point onto the feasible region */
102: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
104: /* Compute the objective function and gradient */
105: TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
106: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
107: if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
109: /* Project the gradient and calculate the norm */
110: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
111: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
113: /* Initialize trust region radius */
114: tao->trust=tao->trust0;
115: if (tao->trust <= 0) {
116: tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
117: }
119: /* Initialize step sizes for the line searches */
120: tron->pgstepsize=1.0;
121: tron->stepsize=tao->trust;
123: tao->reason = TAO_CONTINUE_ITERATING;
124: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
125: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);
126: (*tao->ops->convergencetest)(tao,tao->cnvP);
127: while (tao->reason==TAO_CONTINUE_ITERATING){
129: /* Perform projected gradient iterations */
130: TronGradientProjections(tao,tron);
132: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
133: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
135: tao->ksp_its=0;
136: f=tron->f; delta=tao->trust;
137: tron->n_free_last = tron->n_free;
138: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
140: /* Generate index set (IS) of which bound constraints are active */
141: ISDestroy(&tron->Free_Local);
142: VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
143: ISGetSize(tron->Free_Local, &tron->n_free);
145: /* If no free variables */
146: if (tron->n_free == 0) {
147: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
148: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
149: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);
150: (*tao->ops->convergencetest)(tao,tao->cnvP);
151: if (!tao->reason) {
152: tao->reason = TAO_CONVERGED_STEPTOL;
153: }
154: break;
155: }
156: /* use free_local to mask/submat gradient, hessian, stepdirection */
157: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
158: TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
159: VecSet(tron->DXFree,0.0);
160: VecScale(tron->R, -1.0);
161: TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
162: if (tao->hessian == tao->hessian_pre) {
163: MatDestroy(&tron->Hpre_sub);
164: PetscObjectReference((PetscObject)(tron->H_sub));
165: tron->Hpre_sub = tron->H_sub;
166: } else {
167: TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
168: }
169: KSPReset(tao->ksp);
170: KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
171: while (1) {
173: /* Approximately solve the reduced linear system */
174: KSPCGSetRadius(tao->ksp,delta);
176: KSPSolve(tao->ksp, tron->R, tron->DXFree);
177: KSPGetIterationNumber(tao->ksp,&its);
178: tao->ksp_its+=its;
179: tao->ksp_tot_its+=its;
180: VecSet(tao->stepdirection,0.0);
182: /* Add dxfree matrix to compute step direction vector */
183: VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);
185: VecDot(tao->gradient, tao->stepdirection, &gdx);
186: VecCopy(tao->solution, tron->X_New);
187: VecCopy(tao->gradient, tron->G_New);
189: stepsize=1.0;f_new=f;
191: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
192: TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
193: TaoAddLineSearchCounts(tao);
195: MatMult(tao->hessian, tao->stepdirection, tron->Work);
196: VecAYPX(tron->Work, 0.5, tao->gradient);
197: VecDot(tao->stepdirection, tron->Work, &prered);
198: actred = f_new - f;
199: if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
200: rhok = 1.0;
201: } else if (actred<0) {
202: rhok=PetscAbs(-actred/prered);
203: } else {
204: rhok=0.0;
205: }
207: /* Compare actual improvement to the quadratic model */
208: if (rhok > tron->eta1) { /* Accept the point */
209: /* d = x_new - x */
210: VecCopy(tron->X_New, tao->stepdirection);
211: VecAXPY(tao->stepdirection, -1.0, tao->solution);
213: VecNorm(tao->stepdirection, NORM_2, &xdiff);
214: xdiff *= stepsize;
216: /* Adjust trust region size */
217: if (rhok < tron->eta2 ){
218: delta = PetscMin(xdiff,delta)*tron->sigma1;
219: } else if (rhok > tron->eta4 ){
220: delta= PetscMin(xdiff,delta)*tron->sigma3;
221: } else if (rhok > tron->eta3 ){
222: delta=PetscMin(xdiff,delta)*tron->sigma2;
223: }
224: VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
225: ISDestroy(&tron->Free_Local);
226: VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
227: f=f_new;
228: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
229: VecCopy(tron->X_New, tao->solution);
230: VecCopy(tron->G_New, tao->gradient);
231: break;
232: }
233: else if (delta <= 1e-30) {
234: break;
235: }
236: else {
237: delta /= 4.0;
238: }
239: } /* end linear solve loop */
241: tron->f=f; tron->actred=actred; tao->trust=delta;
242: tao->niter++;
243: TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
244: TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);
245: (*tao->ops->convergencetest)(tao,tao->cnvP);
246: } /* END MAIN LOOP */
247: return(0);
248: }
250: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
251: {
252: PetscErrorCode ierr;
253: PetscInt i;
254: TaoLineSearchConvergedReason ls_reason;
255: PetscReal actred=-1.0,actred_max=0.0;
256: PetscReal f_new;
257: /*
258: The gradient and function value passed into and out of this
259: routine should be current and correct.
261: The free, active, and binding variables should be already identified
262: */
265: for (i=0;i<tron->maxgpits;++i){
267: if (-actred <= (tron->pg_ftol)*actred_max) break;
269: ++tron->gp_iterates;
270: ++tron->total_gp_its;
271: f_new=tron->f;
273: VecCopy(tao->gradient,tao->stepdirection);
274: VecScale(tao->stepdirection,-1.0);
275: TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
276: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
277: &tron->pgstepsize, &ls_reason);
278: TaoAddLineSearchCounts(tao);
280: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
281: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
283: /* Update the iterate */
284: actred = f_new - tron->f;
285: actred_max = PetscMax(actred_max,-(f_new - tron->f));
286: tron->f = f_new;
287: }
288: return(0);
289: }
291: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) {
293: TAO_TRON *tron = (TAO_TRON *)tao->data;
300: if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
302: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);
303: VecCopy(tron->Work,DXL);
304: VecAXPY(DXL,-1.0,tao->gradient);
305: VecSet(DXU,0.0);
306: VecPointwiseMax(DXL,DXL,DXU);
308: VecCopy(tao->gradient,DXU);
309: VecAXPY(DXU,-1.0,tron->Work);
310: VecSet(tron->Work,0.0);
311: VecPointwiseMin(DXU,tron->Work,DXU);
312: return(0);
313: }
315: /*------------------------------------------------------------*/
316: /*MC
317: TAOTRON - The TRON algorithm is an active-set Newton trust region method
318: for bound-constrained minimization.
320: Options Database Keys:
321: + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
322: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
324: Level: beginner
325: M*/
326: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
327: {
328: TAO_TRON *tron;
330: const char *morethuente_type = TAOLINESEARCHMT;
333: tao->ops->setup = TaoSetup_TRON;
334: tao->ops->solve = TaoSolve_TRON;
335: tao->ops->view = TaoView_TRON;
336: tao->ops->setfromoptions = TaoSetFromOptions_TRON;
337: tao->ops->destroy = TaoDestroy_TRON;
338: tao->ops->computedual = TaoComputeDual_TRON;
340: PetscNewLog(tao,&tron);
341: tao->data = (void*)tron;
343: /* Override default settings (unless already changed) */
344: if (!tao->max_it_changed) tao->max_it = 50;
345: if (!tao->trust0_changed) tao->trust0 = 1.0;
346: if (!tao->steptol_changed) tao->steptol = 0.0;
348: /* Initialize pointers and variables */
349: tron->n = 0;
350: tron->maxgpits = 3;
351: tron->pg_ftol = 0.001;
353: tron->eta1 = 1.0e-4;
354: tron->eta2 = 0.25;
355: tron->eta3 = 0.50;
356: tron->eta4 = 0.90;
358: tron->sigma1 = 0.5;
359: tron->sigma2 = 2.0;
360: tron->sigma3 = 4.0;
362: tron->gp_iterates = 0; /* Cumulative number */
363: tron->total_gp_its = 0;
364: tron->n_free = 0;
366: tron->DXFree=NULL;
367: tron->R=NULL;
368: tron->X_New=NULL;
369: tron->G_New=NULL;
370: tron->Work=NULL;
371: tron->Free_Local=NULL;
372: tron->H_sub=NULL;
373: tron->Hpre_sub=NULL;
374: tao->subset_type = TAO_SUBSET_SUBVEC;
376: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
377: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
378: TaoLineSearchSetType(tao->linesearch,morethuente_type);
379: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
380: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
382: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
383: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
384: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
385: KSPSetType(tao->ksp,KSPCGSTCG);
386: return(0);
387: }