Actual source code: gpcg.c

  1: #include <petscksp.h>
  2: #include <../src/tao/quadratic/impls/gpcg/gpcg.h>

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

  7: /*------------------------------------------------------------*/
  8: static PetscErrorCode TaoDestroy_GPCG(Tao tao)
  9: {
 10:   TAO_GPCG       *gpcg = (TAO_GPCG *)tao->data;

 12:   /* Free allocated memory in GPCG structure */
 13:   VecDestroy(&gpcg->B);
 14:   VecDestroy(&gpcg->Work);
 15:   VecDestroy(&gpcg->X_New);
 16:   VecDestroy(&gpcg->G_New);
 17:   VecDestroy(&gpcg->DXFree);
 18:   VecDestroy(&gpcg->R);
 19:   VecDestroy(&gpcg->PG);
 20:   MatDestroy(&gpcg->Hsub);
 21:   MatDestroy(&gpcg->Hsub_pre);
 22:   ISDestroy(&gpcg->Free_Local);
 23:   PetscFree(tao->data);
 24:   return 0;
 25: }

 27: /*------------------------------------------------------------*/
 28: static PetscErrorCode TaoSetFromOptions_GPCG(PetscOptionItems *PetscOptionsObject,Tao tao)
 29: {
 30:   TAO_GPCG       *gpcg = (TAO_GPCG *)tao->data;
 31:   PetscBool      flg;

 33:   PetscOptionsHead(PetscOptionsObject,"Gradient Projection, Conjugate Gradient method for bound constrained optimization");
 34:   PetscOptionsInt("-tao_gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",NULL,gpcg->maxgpits,&gpcg->maxgpits,&flg);
 35:   PetscOptionsTail();
 36:   KSPSetFromOptions(tao->ksp);
 37:   TaoLineSearchSetFromOptions(tao->linesearch);
 38:   return 0;
 39: }

 41: /*------------------------------------------------------------*/
 42: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
 43: {
 44:   TAO_GPCG       *gpcg = (TAO_GPCG *)tao->data;
 45:   PetscBool      isascii;

 47:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 48:   if (isascii) {
 49:     PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",gpcg->total_gp_its);
 50:     PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)gpcg->pg_ftol);
 51:   }
 52:   TaoLineSearchView(tao->linesearch,viewer);
 53:   return 0;
 54: }

 56: /* GPCGObjectiveAndGradient()
 57:    Compute f=0.5 * x'Hx + b'x + c
 58:            g=Hx + b
 59: */
 60: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void*tptr)
 61: {
 62:   Tao            tao = (Tao)tptr;
 63:   TAO_GPCG       *gpcg = (TAO_GPCG*)tao->data;
 64:   PetscReal      f1,f2;

 66:   MatMult(tao->hessian,X,G);
 67:   VecDot(G,X,&f1);
 68:   VecDot(gpcg->B,X,&f2);
 69:   VecAXPY(G,1.0,gpcg->B);
 70:   *f=f1/2.0 + f2 + gpcg->c;
 71:   return 0;
 72: }

 74: /* ---------------------------------------------------------- */
 75: static PetscErrorCode TaoSetup_GPCG(Tao tao)
 76: {
 77:   TAO_GPCG       *gpcg = (TAO_GPCG *)tao->data;

 79:   /* Allocate some arrays */
 80:   if (!tao->gradient) {
 81:       VecDuplicate(tao->solution, &tao->gradient);
 82:   }
 83:   if (!tao->stepdirection) {
 84:       VecDuplicate(tao->solution, &tao->stepdirection);
 85:   }
 86:   if (!tao->XL) {
 87:       VecDuplicate(tao->solution,&tao->XL);
 88:       VecSet(tao->XL,PETSC_NINFINITY);
 89:   }
 90:   if (!tao->XU) {
 91:       VecDuplicate(tao->solution,&tao->XU);
 92:       VecSet(tao->XU,PETSC_INFINITY);
 93:   }

 95:   VecDuplicate(tao->solution,&gpcg->B);
 96:   VecDuplicate(tao->solution,&gpcg->Work);
 97:   VecDuplicate(tao->solution,&gpcg->X_New);
 98:   VecDuplicate(tao->solution,&gpcg->G_New);
 99:   VecDuplicate(tao->solution,&gpcg->DXFree);
100:   VecDuplicate(tao->solution,&gpcg->R);
101:   VecDuplicate(tao->solution,&gpcg->PG);
102:   /*
103:     if (gpcg->ksp_type == GPCG_KSP_NASH) {
104:         KSPSetType(tao->ksp,KSPNASH);
105:       } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
106:         KSPSetType(tao->ksp,KSPSTCG);
107:       } else {
108:         KSPSetType(tao->ksp,KSPGLTR);
109:       }
110:       if (tao->ksp->ops->setfromoptions) {
111:         (*tao->ksp->ops->setfromoptions)(tao->ksp);
112:       }

114:     }
115:   */
116:   return 0;
117: }

119: static PetscErrorCode TaoSolve_GPCG(Tao tao)
120: {
121:   TAO_GPCG                     *gpcg = (TAO_GPCG *)tao->data;
122:   PetscInt                     its;
123:   PetscReal                    actred,f,f_new,gnorm,gdx,stepsize,xtb;
124:   PetscReal                    xtHx;
125:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;


128:   TaoComputeVariableBounds(tao);
129:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
130:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);

132:   /* Using f = .5*x'Hx + x'b + c and g=Hx + b,  compute b,c */
133:   TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
134:   TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
135:   VecCopy(tao->gradient, gpcg->B);
136:   MatMult(tao->hessian,tao->solution,gpcg->Work);
137:   VecDot(gpcg->Work, tao->solution, &xtHx);
138:   VecAXPY(gpcg->B,-1.0,gpcg->Work);
139:   VecDot(gpcg->B,tao->solution,&xtb);
140:   gpcg->c=f-xtHx/2.0-xtb;
141:   if (gpcg->Free_Local) {
142:       ISDestroy(&gpcg->Free_Local);
143:   }
144:   VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local);

146:   /* Project the gradient and calculate the norm */
147:   VecCopy(tao->gradient,gpcg->G_New);
148:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,gpcg->PG);
149:   VecNorm(gpcg->PG,NORM_2,&gpcg->gnorm);
150:   tao->step=1.0;
151:   gpcg->f = f;

153:     /* Check Stopping Condition      */
154:   tao->reason = TAO_CONTINUE_ITERATING;
155:   TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its);
156:   TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step);
157:   (*tao->ops->convergencetest)(tao,tao->cnvP);

159:   while (tao->reason == TAO_CONTINUE_ITERATING) {
160:     /* Call general purpose update function */
161:     if (tao->ops->update) {
162:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
163:     }
164:     tao->ksp_its=0;

166:     GPCGGradProjections(tao);
167:     ISGetSize(gpcg->Free_Local,&gpcg->n_free);

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

171:     KSPReset(tao->ksp);

173:     if (gpcg->n_free > 0) {
174:       /* Create a reduced linear system */
175:       VecDestroy(&gpcg->R);
176:       VecDestroy(&gpcg->DXFree);
177:       TaoVecGetSubVec(tao->gradient,gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R);
178:       VecScale(gpcg->R, -1.0);
179:       TaoVecGetSubVec(tao->stepdirection,gpcg->Free_Local,tao->subset_type, 0.0, &gpcg->DXFree);
180:       VecSet(gpcg->DXFree,0.0);

182:       TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub);

184:       if (tao->hessian_pre == tao->hessian) {
185:         MatDestroy(&gpcg->Hsub_pre);
186:         PetscObjectReference((PetscObject)gpcg->Hsub);
187:         gpcg->Hsub_pre = gpcg->Hsub;
188:       }  else {
189:         TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre);
190:       }

192:       KSPReset(tao->ksp);
193:       KSPSetOperators(tao->ksp,gpcg->Hsub,gpcg->Hsub_pre);

195:       KSPSolve(tao->ksp,gpcg->R,gpcg->DXFree);
196:       KSPGetIterationNumber(tao->ksp,&its);
197:       tao->ksp_its+=its;
198:       tao->ksp_tot_its+=its;
199:       VecSet(tao->stepdirection,0.0);
200:       VecISAXPY(tao->stepdirection,gpcg->Free_Local,1.0,gpcg->DXFree);

202:       VecDot(tao->stepdirection,tao->gradient,&gdx);
203:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
204:       f_new=f;
205:       TaoLineSearchApply(tao->linesearch,tao->solution,&f_new,tao->gradient,tao->stepdirection,&stepsize,&ls_status);

207:       actred = f_new - f;

209:       /* Evaluate the function and gradient at the new point */
210:       VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU, gpcg->PG);
211:       VecNorm(gpcg->PG, NORM_2, &gnorm);
212:       f=f_new;
213:       ISDestroy(&gpcg->Free_Local);
214:       VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local);
215:     } else {
216:       actred = 0; gpcg->step=1.0;
217:       /* if there were no free variables, no cg method */
218:     }

220:     tao->niter++;
221:     gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
222:     TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its);
223:     TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step);
224:     (*tao->ops->convergencetest)(tao,tao->cnvP);
225:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
226:   }  /* END MAIN LOOP  */

228:   return 0;
229: }

231: static PetscErrorCode GPCGGradProjections(Tao tao)
232: {
233:   TAO_GPCG                       *gpcg = (TAO_GPCG *)tao->data;
234:   PetscInt                       i;
235:   PetscReal                      actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
236:   PetscReal                      f_new,gdx,stepsize;
237:   Vec                            DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work;
238:   Vec                            X=tao->solution,G=tao->gradient;
239:   TaoLineSearchConvergedReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING;

241:   /*
242:      The free, active, and binding variables should be already identified
243:   */
244:   for (i=0;i<gpcg->maxgpits;i++) {
245:     if (-actred <= (gpcg->pg_ftol)*actred_max) break;
246:     VecBoundGradientProjection(G,X,XL,XU,DX);
247:     VecScale(DX,-1.0);
248:     VecDot(DX,G,&gdx);

250:     MatMult(tao->hessian,DX,Work);
251:     VecDot(DX,Work,&gAg);

253:     gpcg->gp_iterates++;
254:     gpcg->total_gp_its++;

256:     gtg=-gdx;
257:     if (PetscAbsReal(gAg) == 0.0) {
258:       alpha = 1.0;
259:     } else {
260:       alpha = PetscAbsReal(gtg/gAg);
261:     }
262:     TaoLineSearchSetInitialStepLength(tao->linesearch,alpha);
263:     f_new=gpcg->f;
264:     TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag);

266:     /* Update the iterate */
267:     actred = f_new - gpcg->f;
268:     actred_max = PetscMax(actred_max,-(f_new - gpcg->f));
269:     gpcg->f = f_new;
270:     ISDestroy(&gpcg->Free_Local);
271:     VecWhichInactive(XL,X,tao->gradient,XU,PETSC_TRUE,&gpcg->Free_Local);
272:   }

274:   gpcg->gnorm=gtg;
275:   return 0;
276: } /* End gradient projections */

278: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
279: {
280:   TAO_GPCG       *gpcg = (TAO_GPCG *)tao->data;

282:   VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work);
283:   VecCopy(gpcg->Work, DXL);
284:   VecAXPY(DXL,-1.0,tao->gradient);
285:   VecSet(DXU,0.0);
286:   VecPointwiseMax(DXL,DXL,DXU);

288:   VecCopy(tao->gradient,DXU);
289:   VecAXPY(DXU,-1.0,gpcg->Work);
290:   VecSet(gpcg->Work,0.0);
291:   VecPointwiseMin(DXU,gpcg->Work,DXU);
292:   return 0;
293: }

295: /*------------------------------------------------------------*/
296: /*MC
297:   TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
298:         conjugate-gradient based method for bound-constrained minimization

300:   Options Database Keys:
301: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
302: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets

304:   Level: beginner
305: M*/
306: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
307: {
308:   TAO_GPCG       *gpcg;

310:   tao->ops->setup = TaoSetup_GPCG;
311:   tao->ops->solve = TaoSolve_GPCG;
312:   tao->ops->view  = TaoView_GPCG;
313:   tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
314:   tao->ops->destroy = TaoDestroy_GPCG;
315:   tao->ops->computedual = TaoComputeDual_GPCG;

317:   PetscNewLog(tao,&gpcg);
318:   tao->data = (void*)gpcg;

320:   /* Override default settings (unless already changed) */
321:   if (!tao->max_it_changed) tao->max_it=500;
322:   if (!tao->max_funcs_changed) tao->max_funcs = 100000;
323: #if defined(PETSC_USE_REAL_SINGLE)
324:   if (!tao->gatol_changed) tao->gatol=1e-6;
325:   if (!tao->grtol_changed) tao->grtol=1e-6;
326: #else
327:   if (!tao->gatol_changed) tao->gatol=1e-12;
328:   if (!tao->grtol_changed) tao->grtol=1e-12;
329: #endif

331:   /* Initialize pointers and variables */
332:   gpcg->n=0;
333:   gpcg->maxgpits = 8;
334:   gpcg->pg_ftol = 0.1;

336:   gpcg->gp_iterates=0; /* Cumulative number */
337:   gpcg->total_gp_its = 0;

339:   /* Initialize pointers and variables */
340:   gpcg->n_bind=0;
341:   gpcg->n_free = 0;
342:   gpcg->n_upper=0;
343:   gpcg->n_lower=0;
344:   gpcg->subset_type = TAO_SUBSET_MASK;
345:   gpcg->Hsub=NULL;
346:   gpcg->Hsub_pre=NULL;

348:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
349:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
350:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
351:   KSPSetType(tao->ksp,KSPNASH);

353:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
354:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
355:   TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG);
356:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao);
357:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
358:   return 0;
359: }