Actual source code: tron.c
1: /*$Id$*/
3: #include "tron.h" /*I "tao_solver.h" I*/
5: /* TRON Routines */
6: static int TaoGradProjections(TAO_SOLVER,TAO_TRON *);
7: static int TronCheckOptimalFace(TaoVec*, TaoVec*, TaoVec*, TaoVec*, TaoVec*,
8: TaoIndexSet*, TaoIndexSet*, TaoTruth *optimal);
10: /*------------------------------------------------------------*/
13: static int TaoSetDown_TRON(TAO_SOLVER tao, void*solver)
14: {
15: TAO_TRON *tron = (TAO_TRON *)solver;
16: int info;
17: /* Free allocated memory in TRON structure */
18: TaoFunctionBegin;
19:
20: info = TaoVecDestroy(tron->X_New);CHKERRQ(info);
21: info = TaoVecDestroy(tron->G_New);CHKERRQ(info);
22: info = TaoVecDestroy(tron->DX);CHKERRQ(info);tron->DX=0;
23: info = TaoVecDestroy(tron->Work);CHKERRQ(info);
24: info = TaoVecDestroy(tron->DXFree);CHKERRQ(info);
25: info = TaoVecDestroy(tron->R);CHKERRQ(info);
26: info = TaoVecDestroy(tron->G);CHKERRQ(info);
27: info = TaoVecDestroy(tron->PG);CHKERRQ(info);
28: info = TaoVecDestroy(tron->XL);CHKERRQ(info);
29: info = TaoVecDestroy(tron->XU);CHKERRQ(info);
30:
31: info = TaoIndexSetDestroy(tron->Free_Local);CHKERRQ(info);
32: info = TaoIndexSetDestroy(tron->TT);CHKERRQ(info);
33: info = TaoMatDestroy(tron->Hsub);CHKERRQ(info);
35: info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
36: info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);
37: info = TaoSetVariableBounds(tao,0,0);CHKERRQ(info);
39: TaoFunctionReturn(0);
40: }
42: /*------------------------------------------------------------*/
45: static int TaoSetOptions_TRON(TAO_SOLVER tao, void*solver)
46: {
47: TAO_TRON *tron = (TAO_TRON *)solver;
48: int info;
49: TaoInt ival;
50: TaoTruth flg;
52: TaoFunctionBegin;
54: info = TaoOptionsHead("Newton Trust Region Method for bound constrained optimization");CHKERRQ(info);
56: info = TaoOptionInt("-tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
57: CHKERRQ(info);
59: info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);
61: info = TaoOptionName("-submatrixfree","Mask full matrix instead of extract submatrices","TaoPetscISType",&flg); CHKERRQ(info);
63: info = TaoOptionsTail();CHKERRQ(info);
64: info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
66: TaoFunctionReturn(0);
67: }
69: /*------------------------------------------------------------*/
72: static int TaoView_TRON(TAO_SOLVER tao,void*solver)
73: {
74: TAO_TRON *tron = (TAO_TRON *)solver;
75: int info;
77: TaoFunctionBegin;
78: /*
79: info = TaoPrintf1(tao," Variables, Total: %d,",tron->n);
80: info = TaoPrintf3(tao,"Free: %d, Binding: %d \n",
81: tron->n_free, tron->n - tron->n_free,
82: tron->n_bind);CHKERRQ(info);
83: info = TaoPrintf1(tao," Equal lower bound: %d,",
84: tron->n_lower);CHKERRQ(info);
85: info = TaoPrintf1(tao," Equal upper bound: %d \n",
86: tron->n_upper);CHKERRQ(info);
87: */
88: info = TaoPrintInt(tao," Total PG its: %d,",tron->total_gp_its);CHKERRQ(info);
89: info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",tron->pg_ftol);CHKERRQ(info);
90: info = TaoLineSearchView(tao);CHKERRQ(info);
91: info = TaoPrintStatement(tao," Linear Solver minimizes quadratic over Trust Region: \n");CHKERRQ(info);
93: TaoFunctionReturn(0);
94: }
97: /* ---------------------------------------------------------- */
100: static int TaoSetUp_TRON(TAO_SOLVER tao, void*solver){
102: int info;
103: TAO_TRON *tron = (TAO_TRON *)solver;
104: TaoVec* X;
105: TaoMat *HH;
106: TaoIndexSet *TIS;
108: TaoFunctionBegin;
109: info = TaoGetSolution(tao,&tron->X);CHKERRQ(info); X=tron->X;
110: info = TaoGetHessian(tao,&tron->H);CHKERRQ(info); HH=tron->H;
112: /* Allocate some arrays */
113: info = X->Clone(&tron->DX); CHKERRQ(info);
114: info = X->Clone(&tron->X_New); CHKERRQ(info);
115: info = X->Clone(&tron->G_New); CHKERRQ(info);
116: info = X->Clone(&tron->Work); CHKERRQ(info);
117: info = X->Clone(&tron->DXFree); CHKERRQ(info);
118: info = X->Clone(&tron->R); CHKERRQ(info);
119: info = X->Clone(&tron->G); CHKERRQ(info);
120: info = X->Clone(&tron->PG); CHKERRQ(info);
121: info = X->Clone(&tron->XL); CHKERRQ(info);
122: info = X->Clone(&tron->XU); CHKERRQ(info);
124: info = TaoSetLagrangianGradientVector(tao,tron->PG);CHKERRQ(info);
125: info = TaoSetStepDirectionVector(tao,tron->DX);CHKERRQ(info);
126: info = TaoSetVariableBounds(tao,tron->XL,tron->XU);CHKERRQ(info);
128: info = X->GetDimension(&tron->n); CHKERRQ(info);
129:
130: info = X->CreateIndexSet(&tron->Free_Local); CHKERRQ(info);
131: info = tron->Free_Local->Duplicate(&tron->TT); CHKERRQ(info);
133: TIS=tron->Free_Local;
134: info = tron->H->CreateReducedMatrix(TIS,TIS,&tron->Hsub); CHKERRQ(info);
135: info = TaoCreateLinearSolver(tao,HH,220,0); CHKERRQ(info);
137: info = TaoCheckFGH(tao);CHKERRQ(info);
139: TaoFunctionReturn(0);
140: }
146: static int TaoSolve_TRON(TAO_SOLVER tao, void*solver){
148: TAO_TRON *tron = (TAO_TRON *)solver;
149: int info;
150: TaoInt lsflag,iter=0;
151: TaoTerminateReason reason;
152: TaoTruth optimal_face=TAO_FALSE,success;
153: double prered,actred,delta,f,f_new,f_full,rhok,gnorm,gdx,xdiff,stepsize;
154: TaoVec *XU, *XL;
155: TaoVec *X, *G;
156: TaoVec *PG=tron->PG;
157: TaoVec *R=tron->R, *DXFree=tron->DXFree;
158: TaoVec *X_New=tron->X_New, *G_New=tron->G_New;
159: TaoVec *DX=tron->DX, *Work=tron->Work;
160: TaoMat *H, *Hsub=tron->Hsub;
161: TaoIndexSet *Free_Local = tron->Free_Local, *TIS=tron->TT;
163: TaoFunctionBegin;
165: // Get initial trust region radius
166: info = TaoGetInitialTrustRegionRadius(tao, &tron->delta); CHKERRQ(info);
167: if (tron->delta <= 0) {
168: SETERRQ(1, "Initial trust region radius must be positive");
169: }
171: // Get vectors we will need
172: info = TaoGetSolution(tao, &X); CHKERRQ(info);
173: info = TaoGetGradient(tao, &G); CHKERRQ(info);
174: info = TaoGetHessian(tao, &H); CHKERRQ(info);
175: info = TaoGetVariableBounds(tao, &XL, &XU); CHKERRQ(info);
177: // Check that upper bound greater than lower bound
178: info = TaoEvaluateVariableBounds(tao, XL, XU); CHKERRQ(info);
180: tron->pgstepsize=1.0;
182: /* Project the current point onto the feasible set */
183: info = X->Median(XL,X,XU); CHKERRQ(info);
184:
185: info = TaoComputeMeritFunctionGradient(tao,X,&tron->f,G);CHKERRQ(info);
186: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
187:
188: /* Project the gradient and calculate the norm */
189: // info = G_New->CopyFrom(G);CHKERRQ(info);
190: info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
191: info = PG->Norm2(&tron->gnorm); CHKERRQ(info);
193: if (tron->delta <= 0.0){
194: tron->delta=TaoMax(tron->gnorm*tron->gnorm,1.0);
195: // tron->delta = TAO_INFINITY;
196: }
198: tron->stepsize=tron->delta;
200: info = TaoMonitor(tao,iter++,tron->f,tron->gnorm,0.0,tron->delta,&reason);
201: CHKERRQ(info);
203: while (reason==TAO_CONTINUE_ITERATING){
204:
205: info = TaoGradProjections(tao,tron); CHKERRQ(info);
207: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
208: info = Free_Local->GetSize(&tron->n_free); CHKERRQ(info);
209: f=tron->f; delta=tron->delta; gnorm=tron->gnorm;
211: if (tron->n_free > 0){
212:
213: info = TaoComputeHessian(tao,X,H);CHKERRQ(info);
215: /* Create a reduced linear system */
216: info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
217: info = R->Negate(); CHKERRQ(info);
218: info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
219: info = DXFree->SetToZero(); CHKERRQ(info);
221: info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);
223: info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
225: while (1) {
227: /* Approximately solve the reduced linear system */
229: info = TaoLinearSolveTrustRegion(tao, Hsub, R, DXFree, delta, &success); CHKERRQ(info);
231: info=DX->SetToZero(); CHKERRQ(info);
232: info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);
234: info = G->Dot(DX,&gdx); CHKERRQ(info);
235: info = PetscInfo1(tao,"Expected decrease in function value: %14.12e\n",gdx); CHKERRQ(info);
237: stepsize=1.0; f_new=f;
238: info = X_New->CopyFrom(X); CHKERRQ(info);
239: info = G_New->CopyFrom(G); CHKERRQ(info);
240:
241: info = TaoLineSearchApply(tao,X_New,G_New,DX,Work,
242: &f_new,&f_full,&stepsize,&lsflag);
243: CHKERRQ(info);
244: info = H->Multiply(DX,Work); CHKERRQ(info);
245: info = Work->Aypx(0.5,G); CHKERRQ(info);
246: info = Work->Dot(DX,&prered); CHKERRQ(info);
247: actred = f_new - f;
248:
249: if (actred<0) rhok=TaoAbsScalar(-actred/prered);
250: else rhok=0.0;
252: /* Compare actual improvement to the quadratic model */
253: if (rhok > tron->eta1) { /* Accept the point */
255: info = DX->Waxpby(1.0,X_New,-1.0, X); CHKERRQ(info);
256: info = DX->Norm2(&xdiff); CHKERRQ(info);
257: xdiff*=stepsize;
259: /* Adjust trust region size */
260: if (rhok < tron->eta2 ){
261: delta = TaoMin(xdiff,delta)*tron->sigma1;
262: } else if (rhok > tron->eta4 ){
263: delta= TaoMin(xdiff,delta)*tron->sigma3;
264: } else if (rhok > tron->eta3 ){
265: delta=TaoMin(xdiff,delta)*tron->sigma2;
266: }
268: info = PG->BoundGradientProjection(G_New,XL,X_New,XU);
269: CHKERRQ(info);
270: info = PG->Norm2(&gnorm); CHKERRQ(info);
271: info = TronCheckOptimalFace(X_New,XL,XU,G_New,PG, Free_Local, TIS,
272: &optimal_face); CHKERRQ(info);
273: if (stepsize < 1 || optimal_face==TAO_FALSE || reason!=TAO_CONTINUE_ITERATING ){
274: f=f_new;
275: info = X->CopyFrom(X_New); CHKERRQ(info);
276: info = G->CopyFrom(G_New); CHKERRQ(info);
277: break;
278: }
279: if (delta<=1e-30){
280: break;
281: }
282: }
283: else if (delta <= 1e-30) {
284: break;
285: }
286: else {
287: delta /= 4.0;
288: }
289: } /* end linear solve loop */
290:
291: } else {
292:
293: actred=0;
294: info = Work->BoundGradientProjection(G,XL,X,XU);
295: CHKERRQ(info);
296: info = Work->Norm2(&gnorm); CHKERRQ(info);
297: /* if there were no free variables, no cg method */
299: }
301: tron->f=f;tron->gnorm=gnorm; tron->actred=actred; tron->delta=delta;
302: info = TaoMonitor(tao,iter,f,gnorm,0.0,delta,&reason); CHKERRQ(info);
303: if (reason!=TAO_CONTINUE_ITERATING) break;
304: iter++;
305:
306: } /* END MAIN LOOP */
308: TaoFunctionReturn(0);
309: }
314: static int TaoGradProjections(TAO_SOLVER tao,TAO_TRON *tron)
315: {
316: int info;
317: TaoInt lsflag=0,i;
318: TaoTruth sameface=TAO_FALSE;
319: double actred=-1.0,actred_max=0.0;
320: double f_new, f_full;
321: TaoVec *DX=tron->DX,*XL=tron->XL,*XU=tron->XU,*Work=tron->Work;
322: TaoVec *X=tron->X,*G=tron->G;
323: TaoIndexSet *TT1=tron->Free_Local, *TT2=tron->TT, *TT3;
324: /*
325: The gradient and function value passed into and out of this
326: routine should be current and correct.
327:
328: The free, active, and binding variables should be already identified
329: */
330:
331: TaoFunctionBegin;
332:
333: info = TaoGetSolution(tao,&X);CHKERRQ(info);
334: info = TaoGetGradient(tao,&G);CHKERRQ(info);
335: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
337: info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);
339: for (i=0;i<tron->maxgpits;i++){
341: if ( -actred <= (tron->pg_ftol)*actred_max) break;
342:
343: tron->gp_iterates++; tron->total_gp_its++;
344: f_new=tron->f;
346: info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
348: info = TaoLineSearchApply(tao,X,G,DX,Work,
349: &f_new,&f_full,&tron->pgstepsize,&lsflag);
350: CHKERRQ(info);
352: /* Update the iterate */
353: actred = f_new - tron->f;
354: actred_max = TaoMax(actred_max,-(f_new - tron->f));
355: tron->f = f_new;
357: info = TT2->WhichBetween(XL,X,XU); CHKERRQ(info);
358: info = TT2->IsSame(TT1,&sameface); CHKERRQ(info);
359: if (sameface==TAO_TRUE) {
360: break;
361: } else {
362: // info = TT1->WhichBetween(XL,X,XU); CHKERRQ(info);
363: TT3=TT2;
364: TT2=TT1;
365: TT1=TT3;
366: }
368: }
369:
370: TaoFunctionReturn(0);
371: }
376: static int TronCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
377: TaoIndexSet*Free_Local, TaoIndexSet*TT,
378: TaoTruth *optimal)
379: {
380: int info;
381: TaoInt n_free;
382: double rr;
383: TaoTruth same;
385: TaoFunctionBegin;
386: *optimal = TAO_FALSE;
388: /* Also check to see if the active set is the same */
390: info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
391: info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
392: info = Free_Local->GetSize(&n_free); CHKERRQ(info);
393: if (same == TAO_FALSE){
394: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
395: *optimal = TAO_FALSE;
396: TaoFunctionReturn(0);
397: } else {
398: *optimal = TAO_TRUE;
399: }
401: info = W->CopyFrom(PG); CHKERRQ(info);
402: info = W->Negate(); CHKERRQ(info);
404: info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
405: info = W->Axpy(1.0,PG); CHKERRQ(info);
407: info = W->Norm2(&rr); CHKERRQ(info);
408: if (rr>0) *optimal = TAO_FALSE;
410: *optimal = TAO_FALSE;
411: /*
412: info = tron->TT->whichNonNegative(W); CHKERRQ(info);
413: info = tron->TT->GetSize(&n); CHKERRQ(info);
414: if (n==0) *optimal = TAO_TRUE;
415: */
416: TaoFunctionReturn(0);
417: }
423: int TaoDefaultMonitor_TRON(TAO_SOLVER tao,void *dummy)
424: {
425: int info;
426: TaoInt its,nfree,nbind;
427: double fct,gnorm;
428: TAO_TRON *tron;
430: TaoFunctionBegin;
431: info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
432: info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
433: if (tron){
434: nfree=tron->n_free;
435: nbind=tron->n_bind;
436: info=TaoPrintInt(tao,"iter = %d,",its); CHKERRQ(info);
437: info=TaoPrintDouble(tao," Function value: %g,",fct); CHKERRQ(info);
438: info=TaoPrintDouble(tao," Residual: %g \n",gnorm);CHKERRQ(info);
439:
440: info=TaoPrintInt(tao," free vars = %d,",nfree); CHKERRQ(info);
441: info=TaoPrintInt(tao," binding vars = %d\n",nbind); CHKERRQ(info);
442: }
443: TaoFunctionReturn(0);
444: }
446: int TaoSetMaxGPIts(TAO_SOLVER tao, int its){
447: int info;
448: TAO_TRON *tron;
450: TaoFunctionBegin;
452: info = TaoGetSolverContext(tao,"tao_tron",(void**)&tron); CHKERRQ(info);
453: if (tron){
454: tron->maxgpits = its;
455: }
456: TaoFunctionReturn(0);
457: }
462: static int TaoGetDualVariables_TRON(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void *solver){
464: TAO_TRON *tron = (TAO_TRON *) solver;
465: TaoVec *G=tron->G,*GP=tron->Work;
466: TaoVec *X,*XL,*XU;
467: int info;
469: TaoFunctionBegin;
470: info = TaoGetSolution(tao,&X); CHKERRQ(info);
471: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
472: info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
474: info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
475: info = DXU->SetToZero(); CHKERRQ(info);
476: info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);
478: info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
479: info = GP->SetToZero(); CHKERRQ(info);
480: info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);
482: TaoFunctionReturn(0);
483: }
485: /*------------------------------------------------------------*/
489: int TaoCreate_TRON(TAO_SOLVER tao)
490: {
491: TAO_TRON *tron;
492: int info;
494: TaoFunctionBegin;
496: info = TaoNew(TAO_TRON,&tron); CHKERRQ(info);
497: info = PetscLogObjectMemory(tao,sizeof(TAO_TRON)); CHKERRQ(info);
499: info=TaoSetTaoSolveRoutine(tao,TaoSolve_TRON,(void*)tron); CHKERRQ(info);
500: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_TRON,TaoSetDown_TRON); CHKERRQ(info);
501: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_TRON); CHKERRQ(info);
502: info=TaoSetTaoViewRoutine(tao,TaoView_TRON); CHKERRQ(info);
503: info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_TRON); CHKERRQ(info);
505: info = TaoSetMaximumIterates(tao, 50); CHKERRQ(info);
506: info = TaoSetTolerances(tao, 1e-10, 1e-10, 0, 0); CHKERRQ(info);
508: info = TaoSetTrustRegionRadius(tao, 1.0); CHKERRQ(info);
509: info = TaoSetTrustRegionTolerance(tao, 1.0e-12); CHKERRQ(info);
511: /* Initialize pointers and variables */
512: tron->n = 0;
513: tron->delta = -1.0;
514: tron->maxgpits = 3;
515: tron->pg_ftol = 0.001;
517: tron->eta1 = 1.0e-4;
518: tron->eta2 = 0.25;
519: tron->eta3 = 0.50;
520: tron->eta4 = 0.90;
522: tron->sigma1 = 0.5;
523: tron->sigma2 = 2.0;
524: tron->sigma3 = 4.0;
526: tron->gp_iterates = 0; /* Cumulative number */
527: tron->cgits = 0; /* Current iteration */
528: tron->total_gp_its = 0;
529: tron->cg_iterates = 0;
530: tron->total_cgits = 0;
531:
532: tron->n_bind = 0;
533: tron->n_free = 0;
534: tron->n_upper = 0;
535: tron->n_lower = 0;
537: tron->DX=0;
538: tron->DXFree=0;
539: tron->R=0;
540: tron->X_New=0;
541: tron->G_New=0;
542: tron->Work=0;
543: tron->Free_Local=0;
544: tron->TT=0;
545: tron->Hsub=0;
547: info = TaoCreateMoreThuenteBoundLineSearch(tao, 0 , 0); CHKERRQ(info);
548: TaoFunctionReturn(0);
549: }