Actual source code: gpcg.c
1: /*$Id: gpcg.c 1.93 05/05/09 11:09:17-05:00 sarich@zorak.(none) $*/
3: #include "gpcg.h" /*I "tao_solver.h" I*/
5: char gpcgname[]="tao_gpcg";
6: TaoMethod gpcgtypename = gpcgname;
8: static int TaoGradProjections(TAO_SOLVER, TAO_GPCG *);
9: static int GPCGCheckOptimalFace(TaoVec*, TaoVec*, TaoVec*, TaoVec*, TaoVec*,
10: TaoIndexSet*, TaoIndexSet*, TaoTruth *optimal);
12: /*------------------------------------------------------------*/
15: static int TaoSetDown_GPCG(TAO_SOLVER tao, void*solver)
16: {
17: TAO_GPCG *gpcg = (TAO_GPCG *)solver;
18: int info;
19: /* Free allocated memory in GPCG structure */
20: TaoFunctionBegin;
21:
22: info = TaoVecDestroy(gpcg->X_New);CHKERRQ(info);
23: info = TaoVecDestroy(gpcg->G_New);CHKERRQ(info);
24: info = TaoVecDestroy(gpcg->DX);CHKERRQ(info);gpcg->DX=0;
25: info = TaoVecDestroy(gpcg->Work);CHKERRQ(info);
26: info = TaoVecDestroy(gpcg->DXFree);CHKERRQ(info);
27: info = TaoVecDestroy(gpcg->R);CHKERRQ(info);
28: info = TaoVecDestroy(gpcg->B);CHKERRQ(info);
29: info = TaoVecDestroy(gpcg->G);CHKERRQ(info);
30: info = TaoVecDestroy(gpcg->PG);CHKERRQ(info);
31: info = TaoVecDestroy(gpcg->XL);CHKERRQ(info);
32: info = TaoVecDestroy(gpcg->XU);CHKERRQ(info);
33:
34: info = TaoIndexSetDestroy(gpcg->Free_Local);CHKERRQ(info);
35: info = TaoIndexSetDestroy(gpcg->TT);CHKERRQ(info);
36: info = TaoMatDestroy(gpcg->Hsub);CHKERRQ(info);
37:
38: info = TaoDestroyLinearSolver(tao);CHKERRQ(info);
40: TaoFunctionReturn(0);
41: }
43: /*------------------------------------------------------------*/
46: static int TaoSetFromOptions_GPCG(TAO_SOLVER tao, void*solver)
47: {
48: TAO_GPCG *gpcg = (TAO_GPCG *)solver;
49: int info,ival;
50: TaoTruth flg;
52: TaoFunctionBegin;
53: info = TaoOptionsHead("Gradient Projection, Conjugate Gradient method for bound constrained optimization");CHKERRQ(info);
55: info=TaoOptionInt("-gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",0,gpcg->maxgpits,&gpcg->maxgpits,&flg);
56: CHKERRQ(info);
58: info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);
60: 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);
65: /* info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */
67: TaoFunctionReturn(0);
68: }
70: /*------------------------------------------------------------*/
73: static int TaoView_GPCG(TAO_SOLVER tao,void*solver)
74: {
75: TAO_GPCG *gpcg = (TAO_GPCG *)solver;
76: int info;
78: TaoFunctionBegin;
80: info = TaoPrintInt(tao," Total PG its: %d,",gpcg->total_gp_its);CHKERRQ(info);
81: info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",gpcg->pg_ftol);CHKERRQ(info);
82: info = TaoLineSearchView(tao);CHKERRQ(info);
84: TaoFunctionReturn(0);
85: }
87: /* ---------------------------------------------------------- */
90: int TaoGPCGComputeFunctionGradient(TAO_SOLVER tao, TaoVec *XX, double *ff, TaoVec *GG){
91: TAO_GPCG *gpcg;
92: TaoMat *HH;
93: int info;
94: double f1,f2;
96: TaoFunctionBegin;
97: info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
98: if (gpcg==0){TaoFunctionReturn(0);}
99: info = TaoGetHessian(tao,&HH);CHKERRQ(info);
100: info = HH->Multiply(XX,GG);CHKERRQ(info);
101: info = XX->Dot(GG,&f1);CHKERRQ(info);
102: info = XX->Dot(gpcg->B,&f2);CHKERRQ(info);
103: info = GG->Axpy(1.0,gpcg->B);CHKERRQ(info);
104: *ff=f1/2.0 + f2 + gpcg->c;
105: TaoFunctionReturn(0);
106: }
108: /* ---------------------------------------------------------- */
111: static int TaoSetUp_GPCG(TAO_SOLVER tao,void*solver){
113: int n,info;
114: TAO_GPCG *gpcg = (TAO_GPCG *) solver;
115: TaoVec *X;
116: TaoMat *HH;
117: TaoIndexSet *TIS;
118: TaoLinearSolver *ksp;
120: TaoFunctionBegin;
122: info = TaoGetSolution(tao,&X);CHKERRQ(info); gpcg->X=X;
123: info = TaoGetHessian(tao,&HH);CHKERRQ(info); gpcg->H=HH;
125: /* Allocate some arrays */
126: info=X->Clone(&gpcg->DX); CHKERRQ(info);
127: info=X->Clone(&gpcg->B); CHKERRQ(info);
128: info=X->Clone(&gpcg->Work); CHKERRQ(info);
129: info=X->Clone(&gpcg->X_New); CHKERRQ(info);
130: info=X->Clone(&gpcg->G_New); CHKERRQ(info);
131: info=X->Clone(&gpcg->DXFree); CHKERRQ(info);
132: info=X->Clone(&gpcg->R); CHKERRQ(info);
133: info=X->Clone(&gpcg->G); CHKERRQ(info);
134: info=X->Clone(&gpcg->PG); CHKERRQ(info);
135: info=X->Clone(&gpcg->XL); CHKERRQ(info);
136: info=X->Clone(&gpcg->XU); CHKERRQ(info);
138: info = TaoSetLagrangianGradientVector(tao,gpcg->PG);CHKERRQ(info);
139: info = TaoSetStepDirectionVector(tao,gpcg->DX);CHKERRQ(info);
140: info = TaoSetVariableBounds(tao,gpcg->XL,gpcg->XU);CHKERRQ(info);
142: info = X->GetDimension(&n); CHKERRQ(info);
143: gpcg->n=n;
144: info = TaoCreateLinearSolver(tao,HH,300,0); CHKERRQ(info);
145: info = TaoGetLinearSolver(tao,&ksp); CHKERRQ(info);
146: info = ksp->SetTolerances(1e-2,1e-30,1e30,n);CHKERRQ(info);
148: info = X->CreateIndexSet(&TIS); CHKERRQ(info);
149: gpcg->Free_Local = TIS;
150: info = gpcg->Free_Local->Duplicate(&gpcg->TT); CHKERRQ(info);
152: info = HH->CreateReducedMatrix(TIS,TIS,&gpcg->Hsub); CHKERRQ(info);
154: TaoFunctionReturn(0);
155: }
159: static int TaoSolve_GPCG(TAO_SOLVER tao, void *solver)
160: {
161: TAO_GPCG *gpcg = (TAO_GPCG *)solver ;
162: int info,lsflag,iter=0;
163: TaoTruth optimal_face=TAO_FALSE,success;
164: double actred,f,f_new,gnorm,gdx,stepsize;
165: double c;
166: TaoVec *XU, *XL;
167: TaoVec *X, *G=gpcg->G , *B=gpcg->B, *PG=gpcg->PG;
168: TaoVec *R=gpcg->R, *DXFree=gpcg->DXFree;
169: TaoVec *X_New=gpcg->X_New, *G_New=gpcg->G_New;
170: TaoVec *DX=gpcg->DX, *Work=gpcg->Work;
171: TaoMat *H, *Hsub=gpcg->Hsub;
172: TaoIndexSet *Free_Local = gpcg->Free_Local, *TIS=gpcg->TT;
173: TaoTerminateReason reason;
175: TaoFunctionBegin;
177: /* Check if upper bound greater than lower bound. */
178: info = TaoGetSolution(tao,&X);CHKERRQ(info);
179: info = TaoGetHessian(tao,&H);CHKERRQ(info);
181: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
182: info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
183: info = X->Median(XL,X,XU); CHKERRQ(info);
185: info = TaoComputeHessian(tao,X,H); CHKERRQ(info);
186: info = TaoComputeFunctionGradient(tao,X,&f,B);
187: CHKERRQ(info);
189: /* Compute quadratic representation */
190: info = H->Multiply(X,Work); CHKERRQ(info);
191: info = X->Dot(Work,&c); CHKERRQ(info);
192: info = B->Axpy(-1.0,Work); CHKERRQ(info);
193: info = X->Dot(B,&stepsize); CHKERRQ(info);
194: gpcg->c=f-c/2.0-stepsize;
196: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
197:
198: info = TaoGPCGComputeFunctionGradient(tao, X, &gpcg->f , G);
199:
200: /* Project the gradient and calculate the norm */
201: info = G_New->CopyFrom(G);CHKERRQ(info);
202: info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
203: info = PG->Norm2(&gpcg->gnorm); CHKERRQ(info);
204: gpcg->step=1.0;
206: /* Check Stopping Condition */
207: info=TaoMonitor(tao,iter++,gpcg->f,gpcg->gnorm,0,gpcg->step,&reason); CHKERRQ(info);
209: while (reason == TAO_CONTINUE_ITERATING){
211: info = TaoGradProjections(tao, gpcg); CHKERRQ(info);
213: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
214: info = Free_Local->GetSize(&gpcg->n_free); CHKERRQ(info);
215: f=gpcg->f; gnorm=gpcg->gnorm;
217: if (gpcg->n_free > 0){
218:
219: /* Create a reduced linear system */
220: info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
221: info = R->Negate(); CHKERRQ(info);
222: info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
223: info = DXFree->SetToZero(); CHKERRQ(info);
225: info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);
227: info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
229: while (1){
230:
231: /* Approximately solve the reduced linear system */
232: info = TaoLinearSolve(tao,Hsub,R,DXFree,&success);CHKERRQ(info);
233:
234: info=DX->SetToZero(); CHKERRQ(info);
235: info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);
237: info = G->Dot(DX,&gdx); CHKERRQ(info);
239: stepsize=1.0; f_new=f;
240: info = X_New->CopyFrom(X); CHKERRQ(info);
241: info = G_New->CopyFrom(G); CHKERRQ(info);
242:
243: info = TaoLineSearchApply(tao,X_New,G_New,DX,Work,
244: &f_new,&stepsize,&lsflag);
245: CHKERRQ(info);
247: actred = f_new - f;
248:
249: /* Evaluate the function and gradient at the new point */
250: info = PG->BoundGradientProjection(G_New,XL,X_New,XU);
251: CHKERRQ(info);
252: info = PG->Norm2(&gnorm); CHKERRQ(info);
254: info = GPCGCheckOptimalFace(X_New,XL,XU,PG,Work, Free_Local, TIS,
255: &optimal_face); CHKERRQ(info);
257: info = TaoMonitor(tao,iter,f_new,gnorm,0.0,stepsize,&reason); CHKERRQ(info);
259: if (stepsize < 1 || optimal_face==TAO_FALSE || reason!=TAO_CONTINUE_ITERATING){
260: f=f_new;
261: info = X->CopyFrom(X_New); CHKERRQ(info);
262: info = G->CopyFrom(G_New); CHKERRQ(info);
263: break;
264: }
266: } /* end linear solve loop */
267:
268:
269: } else {
270:
271: actred = 0;
272: /* if there were no free variables, no cg method */
273: info = TaoMonitor(tao,iter,f,gnorm,0.0,1.0,&reason); CHKERRQ(info);
275: }
277: gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
278: if (reason!=TAO_CONTINUE_ITERATING) break;
279: iter++;
282: } /* END MAIN LOOP */
284: TaoFunctionReturn(0);
285: }
289: static int TaoGradProjections(TAO_SOLVER tao, TAO_GPCG *gpcg)
290: {
291: int i,lsflag=0,info;
292: TaoTruth optimal_face=TAO_FALSE;
293: double actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
294: double f_new,gdx;
295: TaoMat *H;
296: TaoVec *DX=gpcg->DX,*XL=gpcg->XL,*XU=gpcg->XU,*Work=gpcg->Work;
297: TaoVec *X=gpcg->X,*G=gpcg->G;
298: /*
299: The gradient and function value passed into and out of this
300: routine should be current and correct.
301:
302: The free, active, and binding variables should be already identified
303: */
304:
305: TaoFunctionBegin;
306:
307: info = TaoGetSolution(tao,&X);CHKERRQ(info);
308: info = TaoGetHessian(tao,&H);CHKERRQ(info);
309: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
311: for (i=0;i<gpcg->maxgpits;i++){
313: if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
314:
315: info = DX->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
316: info = DX->Negate(); CHKERRQ(info);
317: info = DX->Dot(G,&gdx); CHKERRQ(info);
319: info= H->Multiply(DX,Work); CHKERRQ(info);
320: info= DX->Dot(Work,&gAg); CHKERRQ(info);
321:
322: gpcg->gp_iterates++; gpcg->total_gp_its++;
323:
324: gtg=-gdx;
325: alpha = TaoAbsDouble(gtg/gAg);
326: gpcg->stepsize = alpha; f_new=gpcg->f;
328: info = TaoLineSearchApply(tao,X,G,DX,Work,
329: &f_new,&gpcg->stepsize,&lsflag);
330: CHKERRQ(info);
332: /* Update the iterate */
333: actred = f_new - gpcg->f;
334: actred_max = TaoMax(actred_max,-(f_new - gpcg->f));
335: gpcg->f = f_new;
336: info = GPCGCheckOptimalFace(X,XL,XU,G,Work,gpcg->Free_Local,gpcg->TT,
337: &optimal_face); CHKERRQ(info);
339: if ( optimal_face == TAO_TRUE ) break;
341: }
342:
343: gpcg->gnorm=gtg;
344: TaoFunctionReturn(0);
346: } /* End gradient projections */
351: static int GPCGCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
352: TaoIndexSet*Free_Local, TaoIndexSet*TT,
353: TaoTruth *optimal)
354: {
355: int info,n_free;
356: double rr;
357: TaoTruth same;
359: TaoFunctionBegin;
360: *optimal = TAO_FALSE;
362: /* Also check to see if the active set is the same */
364: info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
365: info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
366: info = Free_Local->GetSize(&n_free); CHKERRQ(info);
367: if (same == TAO_FALSE){
368: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
369: *optimal = TAO_FALSE;
370: TaoFunctionReturn(0);
371: } else {
372: *optimal = TAO_TRUE;
373: }
376: info = W->CopyFrom(PG); CHKERRQ(info);
377: info = W->Negate(); CHKERRQ(info);
379: info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
380: info = W->Axpy(1.0,PG); CHKERRQ(info);
382: info = W->Norm2(&rr); CHKERRQ(info);
383: if (rr>0) *optimal = TAO_FALSE;
385: *optimal = TAO_FALSE;
386: /*
387: info = gpcg->TT->whichNonNegative(W); CHKERRQ(info);
388: info = gpcg->TT->GetSize(&n); CHKERRQ(info);
389: if (n==0) *optimal = TAO_TRUE;
390: */
391: TaoFunctionReturn(0);
392: }
394: int TaoDefaultMonitor_GPCG(TAO_SOLVER tao,void *dummy)
395: {
396: TAO_GPCG *gpcg;
397: double fct,gnorm;
398: int its,nfree,info;
400: TaoFunctionBegin;
401: info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
402: info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
403: if (gpcg==0){TaoFunctionReturn(0);}
404: nfree=gpcg->n_free;
405: info = TaoPrintInt(tao,"iter: %d,",its);CHKERRQ(info);
406: info = TaoPrintDouble(tao," Fcn value: %g,",fct);CHKERRQ(info);
407: info = TaoPrintDouble(tao," PGrad. norm: %g, ",gnorm);CHKERRQ(info);
408: info = TaoPrintInt(tao,"free vars:%d \n",nfree);CHKERRQ(info);
409: TaoFunctionReturn(0);
410: }
414: int TaoGetDualVariables_GPCG(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
415: {
417: TAO_GPCG *gpcg = (TAO_GPCG *) solver;
418: TaoVec *G=gpcg->G,*GP=gpcg->Work;
419: TaoVec *X,*XL,*XU;
420: int info;
422: TaoFunctionBegin;
423: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
424: info = TaoGetSolution(tao,&X); CHKERRQ(info);
425: info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
427: info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
428: info = DXU->SetToZero(); CHKERRQ(info);
429: info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);
431: info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
432: info = GP->SetToZero(); CHKERRQ(info);
433: info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);
435: TaoFunctionReturn(0);
436: }
438: /*------------------------------------------------------------*/
442: int TaoCreate_GPCG(TAO_SOLVER tao)
443: {
444: TAO_GPCG *gpcg;
445: int info;
447: TaoFunctionBegin;
449: info = TaoNew(TAO_GPCG,&gpcg); CHKERRQ(info);
450: info = PetscLogObjectMemory(tao,sizeof(TAO_GPCG)); CHKERRQ(info);
452: info=TaoSetTaoSolveRoutine(tao,TaoSolve_GPCG,(void*)gpcg); CHKERRQ(info);
453: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_GPCG,TaoSetDown_GPCG); CHKERRQ(info);
454: info=TaoSetTaoOptionsRoutine(tao,TaoSetFromOptions_GPCG); CHKERRQ(info);
455: info=TaoSetTaoViewRoutine(tao,TaoView_GPCG); CHKERRQ(info);
456: info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_GPCG); CHKERRQ(info);
458: info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
459: info = TaoSetMaximumFunctionEvaluations(tao,100000); CHKERRQ(info);
460: info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);
462: /* Initialize pointers and variables */
463: gpcg->n=0;
464: gpcg->maxgpits = 8;
465: gpcg->pg_ftol = 0.1;
467: gpcg->gp_iterates=0; /* Cumulative number */
468: gpcg->total_gp_its = 0;
469:
470: /* Initialize pointers and variables */
471: gpcg->n_bind=0;
472: gpcg->n_free = 0;
473: gpcg->n_upper=0;
474: gpcg->n_lower=0;
476: // info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
477: info = TaoGPCGCreateLineSearch(tao); CHKERRQ(info);
479: TaoFunctionReturn(0);
480: }