Actual source code: gpcg.c

  1: #include "private/kspimpl.h"
  2: #include "gpcg.h"        /*I "gpcg.h" I*/


  5: static PetscErrorCode GPCGGradProjections(TaoSolver tao);
  6: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);

  8: /*------------------------------------------------------------*/
 11: static PetscErrorCode TaoDestroy_GPCG(TaoSolver tao)
 12: {
 13:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
 14:   PetscErrorCode      ierr;
 15:   /* Free allocated memory in GPCG structure */
 17:   
 18:   VecDestroy(&gpcg->B);
 19:   VecDestroy(&gpcg->Work);
 20:   VecDestroy(&gpcg->X_New);
 21:   VecDestroy(&gpcg->G_New);
 22:   VecDestroy(&gpcg->DXFree);
 23:   VecDestroy(&gpcg->R);
 24:   VecDestroy(&gpcg->PG);
 25:   MatDestroy(&gpcg->Hsub); 
 26:   MatDestroy(&gpcg->Hsub_pre); 
 27:   ISDestroy(&gpcg->Free_Local);
 28:   PetscFree(tao->data); 
 29:   tao->data = PETSC_NULL;

 31:   return(0);
 32: }

 34: /*------------------------------------------------------------*/
 37: static PetscErrorCode TaoSetFromOptions_GPCG(TaoSolver tao)
 38: {
 39:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
 40:   PetscErrorCode      ierr;
 41:   PetscBool flg;
 43:   PetscOptionsHead("Gradient Projection, Conjugate Gradient method for bound constrained optimization");

 45:   ierr=PetscOptionsInt("-gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",0,gpcg->maxgpits,&gpcg->maxgpits,&flg);
 46:   

 48:   PetscOptionsTail();

 50:   TaoLineSearchSetFromOptions(tao->linesearch); 

 52:   return(0);
 53: }

 55: /*------------------------------------------------------------*/
 58: static PetscErrorCode TaoView_GPCG(TaoSolver tao, PetscViewer viewer)
 59: {
 60:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
 61:   PetscBool           isascii;
 62:   PetscErrorCode      ierr;

 65:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 66:   if (isascii) {
 67:     PetscViewerASCIIPushTab(viewer); 
 68:     PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",gpcg->total_gp_its);
 69:     PetscViewerASCIIPrintf(viewer,"PG tolerance: %G \n",gpcg->pg_ftol);
 70:     PetscViewerASCIIPopTab(viewer); 
 71:   } else {
 72:     SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO GPCG",((PetscObject)viewer)->type_name);
 73:   }

 75:   TaoLineSearchView(tao->linesearch,viewer);

 77:   return(0);
 78: }

 80: /* GPCGObjectiveAndGradient()
 81:    Compute f=0.5 * x'Hx + b'x + c
 82:            g=Hx + b
 83: */
 86: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void*tptr){
 87:   TaoSolver tao = (TaoSolver)tptr;
 88:   TAO_GPCG *gpcg = (TAO_GPCG*)tao->data;
 90:   PetscReal f1,f2;

 93:   MatMult(tao->hessian,X,G); 
 94:   VecDot(G,X,&f1); 
 95:   VecDot(gpcg->B,X,&f2); 
 96:   VecAXPY(G,1.0,gpcg->B); 
 97:   *f=f1/2.0 + f2 + gpcg->c;
 98:   return(0);
 99: }

101: /* ---------------------------------------------------------- */
104: static PetscErrorCode TaoSetup_GPCG(TaoSolver tao) {

106:   PetscErrorCode      ierr;
107:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;


111:   /* Allocate some arrays */
112:   if (!tao->gradient) {
113:       VecDuplicate(tao->solution, &tao->gradient);
114:       
115:   }
116:   if (!tao->stepdirection) {
117:       VecDuplicate(tao->solution, &tao->stepdirection);
118:       
119:   }
120:   if (!tao->XL) {
121:       VecDuplicate(tao->solution,&tao->XL); 
122:       VecSet(tao->XL,TAO_NINFINITY); 
123:   }
124:   if (!tao->XU) {
125:       VecDuplicate(tao->solution,&tao->XU); 
126:       VecSet(tao->XU,TAO_INFINITY); 
127:   }

129:   ierr=VecDuplicate(tao->solution,&gpcg->B); 
130:   ierr=VecDuplicate(tao->solution,&gpcg->Work); 
131:   ierr=VecDuplicate(tao->solution,&gpcg->X_New); 
132:   ierr=VecDuplicate(tao->solution,&gpcg->G_New); 
133:   ierr=VecDuplicate(tao->solution,&gpcg->DXFree); 
134:   ierr=VecDuplicate(tao->solution,&gpcg->R); 
135:   ierr=VecDuplicate(tao->solution,&gpcg->PG); 
136:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp); 
137:   KSPSetType(tao->ksp,KSPNASH); 
138:   KSPSetFromOptions(tao->ksp); 
139:   /*
140:     if (gpcg->ksp_type == GPCG_KSP_NASH) {
141:         KSPSetType(tao->ksp,KSPNASH); 
142:       }        else if (gpcg->ksp_type == GPCG_KSP_STCG) {
143:         KSPSetType(tao->ksp,KSPSTCG); 
144:       } else {
145:         KSPSetType(tao->ksp,KSPGLTR); 
146:       }          
147:       if (tao->ksp->ops->setfromoptions) {
148:         (*tao->ksp->ops->setfromoptions)(tao->ksp);
149:       }
150:       
151:     }
152:   */


155:   
156:   return(0);
157: }

161: static PetscErrorCode TaoSolve_GPCG(TaoSolver tao)
162: {
163:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
165:   PetscInt iter=0,its;
166:   PetscReal actred,f,f_new,gnorm,gdx,stepsize,xtb;
167:   PetscReal xtHx;
168:   MatStructure structure;
169:   TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
170:   TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

173:   gpcg->Hsub=PETSC_NULL;
174:   gpcg->Hsub_pre=PETSC_NULL;
175:   
176:   TaoComputeVariableBounds(tao); 
177:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); 
178:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU); 
179:   
180:   /* Using f = .5*x'Hx + x'b + c and g=Hx + b,  compute b,c */
181:   TaoComputeHessian(tao,tao->solution,&tao->hessian, &tao->hessian_pre,&structure); 
182:   TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);  
183:   VecCopy(tao->gradient, gpcg->B); 
184:   MatMult(tao->hessian,tao->solution,gpcg->Work); 
185:   VecDot(gpcg->Work, tao->solution, &xtHx); 
186:   VecAXPY(gpcg->B,-1.0,gpcg->Work); 
187:   VecDot(gpcg->B,tao->solution,&xtb); 
188:   gpcg->c=f-xtHx/2.0-xtb;
189:   if (gpcg->Free_Local) {
190:       ISDestroy(&gpcg->Free_Local); 
191:   }
192:   VecWhichBetween(tao->XL,tao->solution,tao->XU,&gpcg->Free_Local); 
193:   
194:   /* Project the gradient and calculate the norm */
195:   VecCopy(tao->gradient,gpcg->G_New); 
196:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,gpcg->PG); 
197:   VecNorm(gpcg->PG,NORM_2,&gpcg->gnorm);
198:   tao->step=1.0;
199:   gpcg->f = f;

201:     /* Check Stopping Condition      */
202:   ierr=TaoMonitor(tao,iter,f,gpcg->gnorm,0.0,tao->step,&reason); 

204:   while (reason == TAO_CONTINUE_ITERATING){

206:     GPCGGradProjections(tao); 
207:     ISGetSize(gpcg->Free_Local,&gpcg->n_free); 

209:     f=gpcg->f; gnorm=gpcg->gnorm; 

211:     KSPReset(tao->ksp); 

213:     if (gpcg->n_free > 0){
214:       
215:       /* Create a reduced linear system */
216:       VecDestroy(&gpcg->R); 
217:       VecDestroy(&gpcg->DXFree); 
218:       VecGetSubVec(tao->gradient,gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R); 
219:       VecScale(gpcg->R, -1.0); 
220:       VecGetSubVec(tao->stepdirection,gpcg->Free_Local,tao->subset_type, 0.0, &gpcg->DXFree); 
221:       VecSet(gpcg->DXFree,0.0); 

223:       
224:       MatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub); 

226:       if (tao->hessian_pre == tao->hessian) {
227:         MatDestroy(&gpcg->Hsub_pre); 
228:         PetscObjectReference((PetscObject)gpcg->Hsub); 
229:         gpcg->Hsub_pre = gpcg->Hsub;
230:       }         else {
231:       MatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre); 
232:       }

234:       KSPReset(tao->ksp); 
235:       KSPSetOperators(tao->ksp,gpcg->Hsub,gpcg->Hsub_pre,DIFFERENT_NONZERO_PATTERN);  

237:       KSPSolve(tao->ksp,gpcg->R,gpcg->DXFree); 
238:       KSPGetIterationNumber(tao->ksp,&its); 
239:       tao->ksp_its+=its;

241:       VecSet(tao->stepdirection,0.0); 
242:       VecReducedXPY(tao->stepdirection,gpcg->DXFree,gpcg->Free_Local);

244:       VecDot(tao->stepdirection,tao->gradient,&gdx); 
245:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 
246:       f_new=f;
247:       TaoLineSearchApply(tao->linesearch,tao->solution,&f_new,tao->gradient,tao->stepdirection,&stepsize,
248:                                 &ls_status);   
249:       
250:       actred = f_new - f;
251:       
252:       /* Evaluate the function and gradient at the new point */      
253:       VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU, gpcg->PG); 
254:       VecNorm(gpcg->PG, NORM_2, &gnorm); 
255:       f=f_new;
256:       ISDestroy(&gpcg->Free_Local); 
257:       VecWhichBetween(tao->XL,tao->solution,tao->XU,&gpcg->Free_Local); 
258:       
259:     } else {

261:       actred = 0; gpcg->step=1.0;
262:       /* if there were no free variables, no cg method */
263:     }

265:     iter++;
266:     TaoMonitor(tao,iter,f,gnorm,0.0,gpcg->step,&reason); 
267:     gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
268:     if (reason!=TAO_CONTINUE_ITERATING) break;


271:   }  /* END MAIN LOOP  */

273:   return(0);
274: }

278: static PetscErrorCode GPCGGradProjections(TaoSolver tao)
279: {
281:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
282:   PetscInt i;
283:   PetscReal actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
284:   PetscReal f_new,gdx,stepsize;
285:   Vec DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work;
286:   Vec X=tao->solution,G=tao->gradient;
287:   TaoLineSearchTerminationReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING;

289:   /*
290:      The free, active, and binding variables should be already identified
291:   */
292:   

295:   for (i=0;i<gpcg->maxgpits;i++){
296:     if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
297:     VecBoundGradientProjection(G,X,XL,XU,DX); 
298:     VecScale(DX,-1.0); 
299:     VecDot(DX,G,&gdx); 

301:     MatMult(tao->hessian,DX,Work); 
302:     VecDot(DX,Work,&gAg); 

304:     gpcg->gp_iterates++; 
305:     gpcg->total_gp_its++;    
306:   
307:     gtg=-gdx;
308:     alpha = PetscAbsReal(gtg/gAg);
309:     TaoLineSearchSetInitialStepLength(tao->linesearch,alpha); 
310:     f_new=gpcg->f;
311:     TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag);
312:     

314:     /* Update the iterate */
315:     actred = f_new - gpcg->f;
316:     actred_max = PetscMax(actred_max,-(f_new - gpcg->f));
317:     gpcg->f = f_new;
318:     ISDestroy(&gpcg->Free_Local); 
319:     VecWhichBetween(XL,X,XU,&gpcg->Free_Local); 
320:   }
321:   
322:   gpcg->gnorm=gtg;
323:   return(0);

325: } /* End gradient projections */




332: static PetscErrorCode TaoComputeDual_GPCG(TaoSolver tao, Vec DXL, Vec DXU)
333: {

335:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

339:   VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work); 

341:   VecCopy(gpcg->Work, DXL); 
342:   VecAXPY(DXL,-1.0,tao->gradient); 
343:   VecSet(DXU,0.0); 
344:   VecPointwiseMax(DXL,DXL,DXU); 

346:   VecCopy(tao->gradient,DXU); 
347:   VecAXPY(DXU,-1.0,gpcg->Work); 
348:   VecSet(gpcg->Work,0.0); 
349:   VecPointwiseMin(DXU,gpcg->Work,DXU); 

351:   return(0);
352: }

354: /*------------------------------------------------------------*/
358: PetscErrorCode TaoCreate_GPCG(TaoSolver tao)
359: {
360:   TAO_GPCG *gpcg;
361:   PetscErrorCode      ierr;

364:   tao->ops->setup = TaoSetup_GPCG;
365:   tao->ops->solve = TaoSolve_GPCG;
366:   tao->ops->view  = TaoView_GPCG;
367:   tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
368:   tao->ops->destroy = TaoDestroy_GPCG;
369:   tao->ops->computedual = TaoComputeDual_GPCG;

371:   PetscNewLog(tao, TAO_GPCG, &gpcg); 
372:   tao->data = (void*)gpcg;

374:   tao->max_it = 500;
375:   tao->max_funcs = 100000;
376:   tao->fatol = 1e-12;
377:   tao->frtol = 1e-12;

379:   /* Initialize pointers and variables */
380:   gpcg->n=0;
381:   gpcg->maxgpits = 8;
382:   gpcg->pg_ftol = 0.1;

384:   gpcg->gp_iterates=0; /* Cumulative number */
385:   gpcg->total_gp_its = 0;
386:  
387:   /* Initialize pointers and variables */
388:   gpcg->n_bind=0;
389:   gpcg->n_free = 0;
390:   gpcg->n_upper=0;
391:   gpcg->n_lower=0;
392:   gpcg->subset_type = TAO_SUBSET_MASK;
393:   //gpcg->ksp_type = GPCG_KSP_STCG;


396:       
397:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
398:   TaoLineSearchSetType(tao->linesearch, TAOLINESEARCH_GPCG); 
399:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao); 

401:   return(0);
402: }