Actual source code: gpcg.c
1: /*$Id$*/
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: TaoFunctionReturn(0);
39: }
41: /*------------------------------------------------------------*/
44: static int TaoSetFromOptions_GPCG(TAO_SOLVER tao, void*solver)
45: {
46: TAO_GPCG *gpcg = (TAO_GPCG *)solver;
47: int info;
48: TaoInt ival;
49: TaoTruth flg;
51: TaoFunctionBegin;
52: info = TaoOptionsHead("Gradient Projection, Conjugate Gradient method for bound constrained optimization");CHKERRQ(info);
54: info=TaoOptionInt("-gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",0,gpcg->maxgpits,&gpcg->maxgpits,&flg);
55: CHKERRQ(info);
57: info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);
59: info = TaoOptionName("-submatrixfree","Mask full matrix instead of extract submatrices","TaoPetscISType",&flg); CHKERRQ(info);
62: info = TaoOptionsTail();CHKERRQ(info);
63: info=TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
65: TaoFunctionReturn(0);
66: }
68: /*------------------------------------------------------------*/
71: static int TaoView_GPCG(TAO_SOLVER tao,void*solver)
72: {
73: TAO_GPCG *gpcg = (TAO_GPCG *)solver;
74: int info;
76: TaoFunctionBegin;
78: info = TaoPrintInt(tao," Total PG its: %d,",gpcg->total_gp_its);CHKERRQ(info);
79: info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",gpcg->pg_ftol);CHKERRQ(info);
80: info = TaoLineSearchView(tao);CHKERRQ(info);
82: TaoFunctionReturn(0);
83: }
85: /* ---------------------------------------------------------- */
88: int TaoGPCGComputeFunctionGradient(TAO_SOLVER tao, TaoVec *XX, double *ff, TaoVec *GG){
89: TAO_GPCG *gpcg;
90: TaoMat *HH;
91: int info;
92: double f1,f2;
94: TaoFunctionBegin;
95: info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
96: if (gpcg==0){TaoFunctionReturn(0);}
97: info = TaoGetHessian(tao,&HH);CHKERRQ(info);
98: info = HH->Multiply(XX,GG);CHKERRQ(info);
99: info = XX->Dot(GG,&f1);CHKERRQ(info);
100: info = XX->Dot(gpcg->B,&f2);CHKERRQ(info);
101: info = GG->Axpy(1.0,gpcg->B);CHKERRQ(info);
102: *ff=f1/2.0 + f2 + gpcg->c;
103: TaoFunctionReturn(0);
104: }
106: /* ---------------------------------------------------------- */
109: static int TaoSetUp_GPCG(TAO_SOLVER tao,void*solver){
111: int info;
112: TaoInt n;
113: TAO_GPCG *gpcg = (TAO_GPCG *) solver;
114: TaoVec *X;
115: TaoMat *HH;
116: TaoIndexSet *TIS;
118: TaoFunctionBegin;
120: info = TaoGetSolution(tao,&X);CHKERRQ(info); gpcg->X=X;
121: info = TaoGetHessian(tao,&HH);CHKERRQ(info); gpcg->H=HH;
123: /* Allocate some arrays */
124: info=X->Clone(&gpcg->DX); CHKERRQ(info);
125: info=X->Clone(&gpcg->B); CHKERRQ(info);
126: info=X->Clone(&gpcg->Work); CHKERRQ(info);
127: info=X->Clone(&gpcg->X_New); CHKERRQ(info);
128: info=X->Clone(&gpcg->G_New); CHKERRQ(info);
129: info=X->Clone(&gpcg->DXFree); CHKERRQ(info);
130: info=X->Clone(&gpcg->R); CHKERRQ(info);
131: info=X->Clone(&gpcg->G); CHKERRQ(info);
132: info=X->Clone(&gpcg->PG); CHKERRQ(info);
133: info=X->Clone(&gpcg->XL); CHKERRQ(info);
134: info=X->Clone(&gpcg->XU); CHKERRQ(info);
136: info = TaoSetLagrangianGradientVector(tao,gpcg->PG);CHKERRQ(info);
137: info = TaoSetStepDirectionVector(tao,gpcg->DX);CHKERRQ(info);
138: info = TaoSetVariableBounds(tao,gpcg->XL,gpcg->XU);CHKERRQ(info);
140: info = X->GetDimension(&n); CHKERRQ(info);
141: gpcg->n=n;
142: info = TaoCreateLinearSolver(tao,HH,300,0); CHKERRQ(info);
144: info = X->CreateIndexSet(&TIS); CHKERRQ(info);
145: gpcg->Free_Local = TIS;
146: info = gpcg->Free_Local->Duplicate(&gpcg->TT); CHKERRQ(info);
148: info = HH->CreateReducedMatrix(TIS,TIS,&gpcg->Hsub); CHKERRQ(info);
150: TaoFunctionReturn(0);
151: }
155: static int TaoSolve_GPCG(TAO_SOLVER tao, void *solver)
156: {
157: TAO_GPCG *gpcg = (TAO_GPCG *)solver ;
158: int info;
159: TaoInt lsflag,iter=0;
160: TaoTruth optimal_face=TAO_FALSE,success;
161: double actred,f,f_new,f_full,gnorm,gdx,stepsize;
162: double c;
163: TaoVec *XU, *XL;
164: TaoVec *X, *G=gpcg->G , *B=gpcg->B, *PG=gpcg->PG;
165: TaoVec *R=gpcg->R, *DXFree=gpcg->DXFree;
166: TaoVec *G_New=gpcg->G_New;
167: TaoVec *DX=gpcg->DX, *Work=gpcg->Work;
168: TaoMat *H, *Hsub=gpcg->Hsub;
169: TaoIndexSet *Free_Local = gpcg->Free_Local, *TIS=gpcg->TT;
170: TaoTerminateReason reason;
172: TaoFunctionBegin;
174: /* Check if upper bound greater than lower bound. */
175: info = TaoGetSolution(tao,&X);CHKERRQ(info);
176: info = TaoGetHessian(tao,&H);CHKERRQ(info);
178: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
179: info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
180: info = X->Median(XL,X,XU); CHKERRQ(info);
182: info = TaoComputeHessian(tao,X,H); CHKERRQ(info);
183: info = TaoComputeFunctionGradient(tao,X,&f,B);
184: CHKERRQ(info);
186: /* Compute quadratic representation */
187: info = H->Multiply(X,Work); CHKERRQ(info);
188: info = X->Dot(Work,&c); CHKERRQ(info);
189: info = B->Axpy(-1.0,Work); CHKERRQ(info);
190: info = X->Dot(B,&stepsize); CHKERRQ(info);
191: gpcg->c=f-c/2.0-stepsize;
193: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
194:
195: info = TaoGPCGComputeFunctionGradient(tao, X, &gpcg->f , G);
196:
197: /* Project the gradient and calculate the norm */
198: info = G_New->CopyFrom(G);CHKERRQ(info);
199: info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
200: info = PG->Norm2(&gpcg->gnorm); CHKERRQ(info);
201: gpcg->step=1.0;
203: /* Check Stopping Condition */
204: info=TaoMonitor(tao,iter++,gpcg->f,gpcg->gnorm,0,gpcg->step,&reason); CHKERRQ(info);
206: while (reason == TAO_CONTINUE_ITERATING){
208: info = TaoGradProjections(tao, gpcg); CHKERRQ(info);
210: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
211: info = Free_Local->GetSize(&gpcg->n_free); CHKERRQ(info);
212: f=gpcg->f; gnorm=gpcg->gnorm;
214: if (gpcg->n_free > 0){
215:
216: /* Create a reduced linear system */
217: info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
218: info = R->Negate(); CHKERRQ(info);
219: info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
220: info = DXFree->SetToZero(); CHKERRQ(info);
222: info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);
224: info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
226: /* Approximately solve the reduced linear system */
227: info = TaoLinearSolve(tao,Hsub,R,DXFree,&success);CHKERRQ(info);
228:
229: info=DX->SetToZero(); CHKERRQ(info);
230: info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);
231:
232: info = G->Dot(DX,&gdx); CHKERRQ(info);
233:
234: stepsize=1.0; f_new=f;
235: info = TaoLineSearchApply(tao,X,G,DX,Work,
236: &f_new,&f_full,&stepsize,&lsflag);
237: CHKERRQ(info);
238:
239: actred = f_new - f;
240:
241: /* Evaluate the function and gradient at the new point */
242: info = PG->BoundGradientProjection(G,XL,X,XU);
243: CHKERRQ(info);
244: info = PG->Norm2(&gnorm); CHKERRQ(info);
245: f=f_new;
246:
247: info = GPCGCheckOptimalFace(X,XL,XU,PG,Work, Free_Local, TIS,
248: &optimal_face); CHKERRQ(info);
249:
250: } else {
251:
252: actred = 0; stepsize=1.0;
253: /* if there were no free variables, no cg method */
255: }
257: info = TaoMonitor(tao,iter,f,gnorm,0.0,stepsize,&reason); CHKERRQ(info);
258: gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
259: if (reason!=TAO_CONTINUE_ITERATING) break;
260: iter++;
263: } /* END MAIN LOOP */
265: TaoFunctionReturn(0);
266: }
270: static int TaoGradProjections(TAO_SOLVER tao, TAO_GPCG *gpcg)
271: {
272: int info;
273: TaoInt lsflag=0,i;
274: TaoTruth optimal_face=TAO_FALSE;
275: double actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
276: double f_new,f_full,gdx;
277: TaoMat *H;
278: TaoVec *DX=gpcg->DX,*XL=gpcg->XL,*XU=gpcg->XU,*Work=gpcg->Work;
279: TaoVec *X=gpcg->X,*G=gpcg->G;
280: /*
281: The gradient and function value passed into and out of this
282: routine should be current and correct.
283:
284: The free, active, and binding variables should be already identified
285: */
286:
287: TaoFunctionBegin;
288:
289: info = TaoGetSolution(tao,&X);CHKERRQ(info);
290: info = TaoGetHessian(tao,&H);CHKERRQ(info);
291: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
293: for (i=0;i<gpcg->maxgpits;i++){
295: if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
296:
297: info = DX->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
298: info = DX->Negate(); CHKERRQ(info);
299: info = DX->Dot(G,&gdx); CHKERRQ(info);
301: info= H->Multiply(DX,Work); CHKERRQ(info);
302: info= DX->Dot(Work,&gAg); CHKERRQ(info);
303:
304: gpcg->gp_iterates++; gpcg->total_gp_its++;
305:
306: gtg=-gdx;
307: alpha = TaoAbsDouble(gtg/gAg);
308: gpcg->stepsize = alpha; f_new=gpcg->f;
310: info = TaoLineSearchApply(tao,X,G,DX,Work,
311: &f_new,&f_full,&gpcg->stepsize,&lsflag);
312: CHKERRQ(info);
314: /* Update the iterate */
315: actred = f_new - gpcg->f;
316: actred_max = TaoMax(actred_max,-(f_new - gpcg->f));
317: gpcg->f = f_new;
318: info = GPCGCheckOptimalFace(X,XL,XU,G,Work,gpcg->Free_Local,gpcg->TT,
319: &optimal_face); CHKERRQ(info);
321: if ( optimal_face == TAO_TRUE ) break;
323: }
324:
325: gpcg->gnorm=gtg;
326: TaoFunctionReturn(0);
328: } /* End gradient projections */
333: static int GPCGCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
334: TaoIndexSet*Free_Local, TaoIndexSet*TT,
335: TaoTruth *optimal)
336: {
337: int info;
338: TaoInt n_free;
339: double rr;
340: TaoTruth same;
342: TaoFunctionBegin;
343: *optimal = TAO_FALSE;
345: /* Also check to see if the active set is the same */
347: info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
348: info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
349: info = Free_Local->GetSize(&n_free); CHKERRQ(info);
350: if (same == TAO_FALSE){
351: info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
352: *optimal = TAO_FALSE;
353: TaoFunctionReturn(0);
354: } else {
355: *optimal = TAO_TRUE;
356: }
359: info = W->CopyFrom(PG); CHKERRQ(info);
360: info = W->Negate(); CHKERRQ(info);
362: info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
363: info = W->Axpy(1.0,PG); CHKERRQ(info);
365: info = W->Norm2(&rr); CHKERRQ(info);
366: if (rr>0) *optimal = TAO_FALSE;
368: *optimal = TAO_FALSE;
369: /*
370: info = gpcg->TT->whichNonNegative(W); CHKERRQ(info);
371: info = gpcg->TT->GetSize(&n); CHKERRQ(info);
372: if (n==0) *optimal = TAO_TRUE;
373: */
374: TaoFunctionReturn(0);
375: }
377: int TaoDefaultMonitor_GPCG(TAO_SOLVER tao,void *dummy)
378: {
379: TAO_GPCG *gpcg;
380: double fct,gnorm;
381: int info;
382: TaoInt its,nfree;
384: TaoFunctionBegin;
385: info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
386: info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
387: if (gpcg==0){TaoFunctionReturn(0);}
388: nfree=gpcg->n_free;
389: info = TaoPrintInt(tao,"iter: %d,",its);CHKERRQ(info);
390: info = TaoPrintDouble(tao," Fcn value: %g,",fct);CHKERRQ(info);
391: info = TaoPrintDouble(tao," PGrad. norm: %g, ",gnorm);CHKERRQ(info);
392: info = TaoPrintInt(tao,"free vars:%d \n",nfree);CHKERRQ(info);
393: TaoFunctionReturn(0);
394: }
398: int TaoGetDualVariables_GPCG(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
399: {
401: TAO_GPCG *gpcg = (TAO_GPCG *) solver;
402: TaoVec *G=gpcg->G,*GP=gpcg->Work;
403: TaoVec *X,*XL,*XU;
404: int info;
406: TaoFunctionBegin;
407: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
408: info = TaoGetSolution(tao,&X); CHKERRQ(info);
409: info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
411: info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
412: info = DXU->SetToZero(); CHKERRQ(info);
413: info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);
415: info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
416: info = GP->SetToZero(); CHKERRQ(info);
417: info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);
419: TaoFunctionReturn(0);
420: }
422: /*------------------------------------------------------------*/
426: int TaoCreate_GPCG(TAO_SOLVER tao)
427: {
428: TAO_GPCG *gpcg;
429: int info;
431: TaoFunctionBegin;
433: info = TaoNew(TAO_GPCG,&gpcg); CHKERRQ(info);
434: info = PetscLogObjectMemory(tao,sizeof(TAO_GPCG)); CHKERRQ(info);
436: info=TaoSetTaoSolveRoutine(tao,TaoSolve_GPCG,(void*)gpcg); CHKERRQ(info);
437: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_GPCG,TaoSetDown_GPCG); CHKERRQ(info);
438: info=TaoSetTaoOptionsRoutine(tao,TaoSetFromOptions_GPCG); CHKERRQ(info);
439: info=TaoSetTaoViewRoutine(tao,TaoView_GPCG); CHKERRQ(info);
440: info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_GPCG); CHKERRQ(info);
442: info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
443: info = TaoSetMaximumFunctionEvaluations(tao,100000); CHKERRQ(info);
444: info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);
446: /* Initialize pointers and variables */
447: gpcg->n=0;
448: gpcg->maxgpits = 8;
449: gpcg->pg_ftol = 0.1;
451: gpcg->gp_iterates=0; /* Cumulative number */
452: gpcg->total_gp_its = 0;
453:
454: /* Initialize pointers and variables */
455: gpcg->n_bind=0;
456: gpcg->n_free = 0;
457: gpcg->n_upper=0;
458: gpcg->n_lower=0;
460: // info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
461: info = TaoGPCGCreateLineSearch(tao); CHKERRQ(info);
463: TaoFunctionReturn(0);
464: }