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: }