Actual source code: tron.c
petsc-3.6.4 2016-04-12
1: #include <../src/tao/bound/impls/tron/tron.h>
2: #include <petsc/private/kspimpl.h>
3: #include <petsc/private/matimpl.h>
4: #include <../src/tao/matrix/submatfree.h>
7: /* TRON Routines */
8: static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*);
9: /*------------------------------------------------------------*/
12: static PetscErrorCode TaoDestroy_TRON(Tao tao)
13: {
14: TAO_TRON *tron = (TAO_TRON *)tao->data;
18: VecDestroy(&tron->X_New);
19: VecDestroy(&tron->G_New);
20: VecDestroy(&tron->Work);
21: VecDestroy(&tron->DXFree);
22: VecDestroy(&tron->R);
23: VecDestroy(&tron->diag);
24: VecScatterDestroy(&tron->scatter);
25: ISDestroy(&tron->Free_Local);
26: MatDestroy(&tron->H_sub);
27: MatDestroy(&tron->Hpre_sub);
28: PetscFree(tao->data);
29: return(0);
30: }
32: /*------------------------------------------------------------*/
35: static PetscErrorCode TaoSetFromOptions_TRON(PetscOptions *PetscOptionsObject,Tao tao)
36: {
37: TAO_TRON *tron = (TAO_TRON *)tao->data;
39: PetscBool flg;
42: PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");
43: PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
44: PetscOptionsTail();
45: TaoLineSearchSetFromOptions(tao->linesearch);
46: KSPSetFromOptions(tao->ksp);
47: return(0);
48: }
50: /*------------------------------------------------------------*/
53: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
54: {
55: TAO_TRON *tron = (TAO_TRON *)tao->data;
56: PetscBool isascii;
57: PetscErrorCode ierr;
60: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
61: if (isascii) {
62: PetscViewerASCIIPushTab(viewer);
63: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
64: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
65: PetscViewerASCIIPopTab(viewer);
66: }
67: return(0);
68: }
71: /* ---------------------------------------------------------- */
74: static PetscErrorCode TaoSetup_TRON(Tao tao)
75: {
77: TAO_TRON *tron = (TAO_TRON *)tao->data;
81: /* Allocate some arrays */
82: VecDuplicate(tao->solution, &tron->diag);
83: VecDuplicate(tao->solution, &tron->X_New);
84: VecDuplicate(tao->solution, &tron->G_New);
85: VecDuplicate(tao->solution, &tron->Work);
86: VecDuplicate(tao->solution, &tao->gradient);
87: VecDuplicate(tao->solution, &tao->stepdirection);
88: if (!tao->XL) {
89: VecDuplicate(tao->solution, &tao->XL);
90: VecSet(tao->XL, PETSC_NINFINITY);
91: }
92: if (!tao->XU) {
93: VecDuplicate(tao->solution, &tao->XU);
94: VecSet(tao->XU, PETSC_INFINITY);
95: }
96: return(0);
97: }
103: static PetscErrorCode TaoSolve_TRON(Tao tao)
104: {
105: TAO_TRON *tron = (TAO_TRON *)tao->data;
106: PetscErrorCode ierr;
107: PetscInt its;
108: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
109: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
110: PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;
113: tron->pgstepsize=1.0;
114: tao->trust = tao->trust0;
115: /* Project the current point onto the feasible set */
116: TaoComputeVariableBounds(tao);
117: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
118: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
120: TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
121: ISDestroy(&tron->Free_Local);
123: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
125: /* Project the gradient and calculate the norm */
126: VecBoundGradientProjection(tao->gradient,tao->solution, tao->XL, tao->XU, tao->gradient);
127: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
129: if (PetscIsInfOrNanReal(tron->f) || PetscIsInfOrNanReal(tron->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
130: if (tao->trust <= 0) {
131: tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
132: }
134: tron->stepsize=tao->trust;
135: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize, &reason);
136: while (reason==TAO_CONTINUE_ITERATING){
137: tao->ksp_its=0;
138: TronGradientProjections(tao,tron);
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: ISGetSize(tron->Free_Local, &tron->n_free);
145: /* If no free variables */
146: if (tron->n_free == 0) {
147: actred=0;
148: PetscInfo(tao,"No free variables in tron iteration.\n");
149: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
150: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);
151: if (!reason) {
152: reason = TAO_CONVERGED_STEPTOL;
153: TaoSetConvergedReason(tao,reason);
154: }
156: 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: KSPSTCGSetRadius(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);
187: if (0) {
188: PetscReal rhs,stepnorm;
189: VecNorm(tron->R,NORM_2,&rhs);
190: VecNorm(tron->DXFree,NORM_2,&stepnorm);
191: PetscPrintf(PETSC_COMM_WORLD,"|rhs|=%g\t|s|=%g\n",(double)rhs,(double)stepnorm);
192: }
195: VecDot(tao->gradient, tao->stepdirection, &gdx);
196: PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",(double)gdx);
198: VecCopy(tao->solution, tron->X_New);
199: VecCopy(tao->gradient, tron->G_New);
201: stepsize=1.0;f_new=f;
203: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
204: TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
205: TaoAddLineSearchCounts(tao);
207: MatMult(tao->hessian, tao->stepdirection, tron->Work);
208: VecAYPX(tron->Work, 0.5, tao->gradient);
209: VecDot(tao->stepdirection, tron->Work, &prered);
210: actred = f_new - f;
211: if (actred<0) {
212: rhok=PetscAbs(-actred/prered);
213: } else {
214: rhok=0.0;
215: }
217: /* Compare actual improvement to the quadratic model */
218: if (rhok > tron->eta1) { /* Accept the point */
219: /* d = x_new - x */
220: VecCopy(tron->X_New, tao->stepdirection);
221: VecAXPY(tao->stepdirection, -1.0, tao->solution);
223: VecNorm(tao->stepdirection, NORM_2, &xdiff);
224: xdiff *= stepsize;
226: /* Adjust trust region size */
227: if (rhok < tron->eta2 ){
228: delta = PetscMin(xdiff,delta)*tron->sigma1;
229: } else if (rhok > tron->eta4 ){
230: delta= PetscMin(xdiff,delta)*tron->sigma3;
231: } else if (rhok > tron->eta3 ){
232: delta=PetscMin(xdiff,delta)*tron->sigma2;
233: }
234: VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
235: if (tron->Free_Local) {
236: ISDestroy(&tron->Free_Local);
237: }
238: VecWhichBetween(tao->XL, tron->X_New, tao->XU, &tron->Free_Local);
239: f=f_new;
240: VecNorm(tao->gradient,NORM_2,&tron->gnorm);
241: VecCopy(tron->X_New, tao->solution);
242: VecCopy(tron->G_New, tao->gradient);
243: break;
244: }
245: else if (delta <= 1e-30) {
246: break;
247: }
248: else {
249: delta /= 4.0;
250: }
251: } /* end linear solve loop */
254: tron->f=f; tron->actred=actred; tao->trust=delta;
255: tao->niter++;
256: TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta, &reason);
257: } /* END MAIN LOOP */
259: return(0);
260: }
265: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
266: {
267: PetscErrorCode ierr;
268: PetscInt i;
269: TaoLineSearchConvergedReason ls_reason;
270: PetscReal actred=-1.0,actred_max=0.0;
271: PetscReal f_new;
272: /*
273: The gradient and function value passed into and out of this
274: routine should be current and correct.
276: The free, active, and binding variables should be already identified
277: */
279: if (tron->Free_Local) {
280: ISDestroy(&tron->Free_Local);
281: }
282: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
284: for (i=0;i<tron->maxgpits;i++){
286: if ( -actred <= (tron->pg_ftol)*actred_max) break;
288: tron->gp_iterates++; tron->total_gp_its++;
289: f_new=tron->f;
291: VecCopy(tao->gradient, tao->stepdirection);
292: VecScale(tao->stepdirection, -1.0);
293: TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
294: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
295: &tron->pgstepsize, &ls_reason);
296: TaoAddLineSearchCounts(tao);
299: /* Update the iterate */
300: actred = f_new - tron->f;
301: actred_max = PetscMax(actred_max,-(f_new - tron->f));
302: tron->f = f_new;
303: if (tron->Free_Local) {
304: ISDestroy(&tron->Free_Local);
305: }
306: VecWhichBetween(tao->XL,tao->solution,tao->XU,&tron->Free_Local);
307: }
309: return(0);
310: }
314: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) {
316: TAO_TRON *tron = (TAO_TRON *)tao->data;
323: if (!tron->Work || !tao->gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
325: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);
326: VecCopy(tron->Work,DXL);
327: VecAXPY(DXL,-1.0,tao->gradient);
328: VecSet(DXU,0.0);
329: VecPointwiseMax(DXL,DXL,DXU);
331: VecCopy(tao->gradient,DXU);
332: VecAXPY(DXU,-1.0,tron->Work);
333: VecSet(tron->Work,0.0);
334: VecPointwiseMin(DXU,tron->Work,DXU);
335: return(0);
336: }
338: /*------------------------------------------------------------*/
339: /*MC
340: TAOTRON - The TRON algorithm is an active-set Newton trust region method
341: for bound-constrained minimization.
343: Options Database Keys:
344: + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
345: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
347: Level: beginner
348: M*/
351: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
352: {
353: TAO_TRON *tron;
355: const char *morethuente_type = TAOLINESEARCHMT;
358: tao->ops->setup = TaoSetup_TRON;
359: tao->ops->solve = TaoSolve_TRON;
360: tao->ops->view = TaoView_TRON;
361: tao->ops->setfromoptions = TaoSetFromOptions_TRON;
362: tao->ops->destroy = TaoDestroy_TRON;
363: tao->ops->computedual = TaoComputeDual_TRON;
365: PetscNewLog(tao,&tron);
366: tao->data = (void*)tron;
368: /* Override default settings (unless already changed) */
369: if (!tao->max_it_changed) tao->max_it = 50;
371: #if defined(PETSC_USE_REAL_SINGLE)
372: if (!tao->fatol_changed) tao->fatol = 1e-6;
373: if (!tao->frtol_changed) tao->frtol = 1e-6;
374: if (!tao->steptol_changed) tao->steptol = 1.0e-6;
375: #else
376: if (!tao->fatol_changed) tao->fatol = 1e-12;
377: if (!tao->frtol_changed) tao->frtol = 1e-12;
378: if (!tao->steptol_changed) tao->steptol = 1.0e-12;
379: #endif
381: if (!tao->trust0_changed) tao->trust0 = 1.0;
383: /* Initialize pointers and variables */
384: tron->n = 0;
385: tron->maxgpits = 3;
386: tron->pg_ftol = 0.001;
388: tron->eta1 = 1.0e-4;
389: tron->eta2 = 0.25;
390: tron->eta3 = 0.50;
391: tron->eta4 = 0.90;
393: tron->sigma1 = 0.5;
394: tron->sigma2 = 2.0;
395: tron->sigma3 = 4.0;
397: tron->gp_iterates = 0; /* Cumulative number */
398: tron->total_gp_its = 0;
399: tron->n_free = 0;
401: tron->DXFree=NULL;
402: tron->R=NULL;
403: tron->X_New=NULL;
404: tron->G_New=NULL;
405: tron->Work=NULL;
406: tron->Free_Local=NULL;
407: tron->H_sub=NULL;
408: tron->Hpre_sub=NULL;
409: tao->subset_type = TAO_SUBSET_SUBVEC;
411: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
412: TaoLineSearchSetType(tao->linesearch,morethuente_type);
413: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
414: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
416: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
417: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
418: KSPSetType(tao->ksp,KSPSTCG);
419: return(0);
420: }