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