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;
13: /* Free allocated memory in GPCG structure */
15: VecDestroy(&gpcg->B);
16: VecDestroy(&gpcg->Work);
17: VecDestroy(&gpcg->X_New);
18: VecDestroy(&gpcg->G_New);
19: VecDestroy(&gpcg->DXFree);
20: VecDestroy(&gpcg->R);
21: VecDestroy(&gpcg->PG);
22: MatDestroy(&gpcg->Hsub);
23: MatDestroy(&gpcg->Hsub_pre);
24: ISDestroy(&gpcg->Free_Local);
25: PetscFree(tao->data);
26: return(0);
27: }
29: /*------------------------------------------------------------*/
30: static PetscErrorCode TaoSetFromOptions_GPCG(PetscOptionItems *PetscOptionsObject,Tao tao)
31: {
32: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
34: PetscBool flg;
37: PetscOptionsHead(PetscOptionsObject,"Gradient Projection, Conjugate Gradient method for bound constrained optimization");
38: PetscOptionsInt("-tao_gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",NULL,gpcg->maxgpits,&gpcg->maxgpits,&flg);
39: PetscOptionsTail();
40: KSPSetFromOptions(tao->ksp);
41: TaoLineSearchSetFromOptions(tao->linesearch);
42: return(0);
43: }
45: /*------------------------------------------------------------*/
46: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
47: {
48: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
49: PetscBool isascii;
53: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
54: if (isascii) {
55: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",gpcg->total_gp_its);
56: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)gpcg->pg_ftol);
57: }
58: TaoLineSearchView(tao->linesearch,viewer);
59: return(0);
60: }
62: /* GPCGObjectiveAndGradient()
63: Compute f=0.5 * x'Hx + b'x + c
64: g=Hx + b
65: */
66: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void*tptr)
67: {
68: Tao tao = (Tao)tptr;
69: TAO_GPCG *gpcg = (TAO_GPCG*)tao->data;
71: PetscReal f1,f2;
74: MatMult(tao->hessian,X,G);
75: VecDot(G,X,&f1);
76: VecDot(gpcg->B,X,&f2);
77: VecAXPY(G,1.0,gpcg->B);
78: *f=f1/2.0 + f2 + gpcg->c;
79: return(0);
80: }
82: /* ---------------------------------------------------------- */
83: static PetscErrorCode TaoSetup_GPCG(Tao tao)
84: {
86: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
89: /* Allocate some arrays */
90: if (!tao->gradient) {
91: VecDuplicate(tao->solution, &tao->gradient);
92: }
93: if (!tao->stepdirection) {
94: VecDuplicate(tao->solution, &tao->stepdirection);
95: }
96: if (!tao->XL) {
97: VecDuplicate(tao->solution,&tao->XL);
98: VecSet(tao->XL,PETSC_NINFINITY);
99: }
100: if (!tao->XU) {
101: VecDuplicate(tao->solution,&tao->XU);
102: VecSet(tao->XU,PETSC_INFINITY);
103: }
105: VecDuplicate(tao->solution,&gpcg->B);
106: VecDuplicate(tao->solution,&gpcg->Work);
107: VecDuplicate(tao->solution,&gpcg->X_New);
108: VecDuplicate(tao->solution,&gpcg->G_New);
109: VecDuplicate(tao->solution,&gpcg->DXFree);
110: VecDuplicate(tao->solution,&gpcg->R);
111: VecDuplicate(tao->solution,&gpcg->PG);
112: /*
113: if (gpcg->ksp_type == GPCG_KSP_NASH) {
114: KSPSetType(tao->ksp,KSPNASH);
115: } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
116: KSPSetType(tao->ksp,KSPSTCG);
117: } else {
118: KSPSetType(tao->ksp,KSPGLTR);
119: }
120: if (tao->ksp->ops->setfromoptions) {
121: (*tao->ksp->ops->setfromoptions)(tao->ksp);
122: }
124: }
125: */
126: return(0);
127: }
129: static PetscErrorCode TaoSolve_GPCG(Tao tao)
130: {
131: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
132: PetscErrorCode ierr;
133: PetscInt its;
134: PetscReal actred,f,f_new,gnorm,gdx,stepsize,xtb;
135: PetscReal xtHx;
136: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
140: TaoComputeVariableBounds(tao);
141: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
142: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
144: /* Using f = .5*x'Hx + x'b + c and g=Hx + b, compute b,c */
145: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
146: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
147: VecCopy(tao->gradient, gpcg->B);
148: MatMult(tao->hessian,tao->solution,gpcg->Work);
149: VecDot(gpcg->Work, tao->solution, &xtHx);
150: VecAXPY(gpcg->B,-1.0,gpcg->Work);
151: VecDot(gpcg->B,tao->solution,&xtb);
152: gpcg->c=f-xtHx/2.0-xtb;
153: if (gpcg->Free_Local) {
154: ISDestroy(&gpcg->Free_Local);
155: }
156: VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local);
158: /* Project the gradient and calculate the norm */
159: VecCopy(tao->gradient,gpcg->G_New);
160: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,gpcg->PG);
161: VecNorm(gpcg->PG,NORM_2,&gpcg->gnorm);
162: tao->step=1.0;
163: gpcg->f = f;
165: /* Check Stopping Condition */
166: tao->reason = TAO_CONTINUE_ITERATING;
167: TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its);
168: TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step);
169: (*tao->ops->convergencetest)(tao,tao->cnvP);
171: while (tao->reason == TAO_CONTINUE_ITERATING) {
172: /* Call general purpose update function */
173: if (tao->ops->update) {
174: (*tao->ops->update)(tao, tao->niter, tao->user_update);
175: }
176: tao->ksp_its=0;
178: GPCGGradProjections(tao);
179: ISGetSize(gpcg->Free_Local,&gpcg->n_free);
181: f=gpcg->f; gnorm=gpcg->gnorm;
183: KSPReset(tao->ksp);
185: if (gpcg->n_free > 0) {
186: /* Create a reduced linear system */
187: VecDestroy(&gpcg->R);
188: VecDestroy(&gpcg->DXFree);
189: TaoVecGetSubVec(tao->gradient,gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R);
190: VecScale(gpcg->R, -1.0);
191: TaoVecGetSubVec(tao->stepdirection,gpcg->Free_Local,tao->subset_type, 0.0, &gpcg->DXFree);
192: VecSet(gpcg->DXFree,0.0);
194: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub);
196: if (tao->hessian_pre == tao->hessian) {
197: MatDestroy(&gpcg->Hsub_pre);
198: PetscObjectReference((PetscObject)gpcg->Hsub);
199: gpcg->Hsub_pre = gpcg->Hsub;
200: } else {
201: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre);
202: }
204: KSPReset(tao->ksp);
205: KSPSetOperators(tao->ksp,gpcg->Hsub,gpcg->Hsub_pre);
207: KSPSolve(tao->ksp,gpcg->R,gpcg->DXFree);
208: KSPGetIterationNumber(tao->ksp,&its);
209: tao->ksp_its+=its;
210: tao->ksp_tot_its+=its;
211: VecSet(tao->stepdirection,0.0);
212: VecISAXPY(tao->stepdirection,gpcg->Free_Local,1.0,gpcg->DXFree);
214: VecDot(tao->stepdirection,tao->gradient,&gdx);
215: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
216: f_new=f;
217: TaoLineSearchApply(tao->linesearch,tao->solution,&f_new,tao->gradient,tao->stepdirection,&stepsize,&ls_status);
219: actred = f_new - f;
221: /* Evaluate the function and gradient at the new point */
222: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU, gpcg->PG);
223: VecNorm(gpcg->PG, NORM_2, &gnorm);
224: f=f_new;
225: ISDestroy(&gpcg->Free_Local);
226: VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&gpcg->Free_Local);
227: } else {
228: actred = 0; gpcg->step=1.0;
229: /* if there were no free variables, no cg method */
230: }
232: tao->niter++;
233: gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
234: TaoLogConvergenceHistory(tao,f,gpcg->gnorm,0.0,tao->ksp_its);
235: TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step);
236: (*tao->ops->convergencetest)(tao,tao->cnvP);
237: if (tao->reason != TAO_CONTINUE_ITERATING) break;
238: } /* END MAIN LOOP */
240: return(0);
241: }
243: static PetscErrorCode GPCGGradProjections(Tao tao)
244: {
245: PetscErrorCode ierr;
246: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
247: PetscInt i;
248: PetscReal actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
249: PetscReal f_new,gdx,stepsize;
250: Vec DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work;
251: Vec X=tao->solution,G=tao->gradient;
252: TaoLineSearchConvergedReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING;
254: /*
255: The free, active, and binding variables should be already identified
256: */
258: for (i=0;i<gpcg->maxgpits;i++) {
259: if (-actred <= (gpcg->pg_ftol)*actred_max) break;
260: VecBoundGradientProjection(G,X,XL,XU,DX);
261: VecScale(DX,-1.0);
262: VecDot(DX,G,&gdx);
264: MatMult(tao->hessian,DX,Work);
265: VecDot(DX,Work,&gAg);
267: gpcg->gp_iterates++;
268: gpcg->total_gp_its++;
270: gtg=-gdx;
271: if (PetscAbsReal(gAg) == 0.0) {
272: alpha = 1.0;
273: } else {
274: alpha = PetscAbsReal(gtg/gAg);
275: }
276: TaoLineSearchSetInitialStepLength(tao->linesearch,alpha);
277: f_new=gpcg->f;
278: TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag);
280: /* Update the iterate */
281: actred = f_new - gpcg->f;
282: actred_max = PetscMax(actred_max,-(f_new - gpcg->f));
283: gpcg->f = f_new;
284: ISDestroy(&gpcg->Free_Local);
285: VecWhichInactive(XL,X,tao->gradient,XU,PETSC_TRUE,&gpcg->Free_Local);
286: }
288: gpcg->gnorm=gtg;
289: return(0);
290: } /* End gradient projections */
292: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
293: {
294: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
298: VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work);
299: VecCopy(gpcg->Work, DXL);
300: VecAXPY(DXL,-1.0,tao->gradient);
301: VecSet(DXU,0.0);
302: VecPointwiseMax(DXL,DXL,DXU);
304: VecCopy(tao->gradient,DXU);
305: VecAXPY(DXU,-1.0,gpcg->Work);
306: VecSet(gpcg->Work,0.0);
307: VecPointwiseMin(DXU,gpcg->Work,DXU);
308: return(0);
309: }
311: /*------------------------------------------------------------*/
312: /*MC
313: TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
314: conjugate-gradient based method for bound-constrained minimization
316: Options Database Keys:
317: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
318: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
320: Level: beginner
321: M*/
322: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
323: {
324: TAO_GPCG *gpcg;
328: tao->ops->setup = TaoSetup_GPCG;
329: tao->ops->solve = TaoSolve_GPCG;
330: tao->ops->view = TaoView_GPCG;
331: tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
332: tao->ops->destroy = TaoDestroy_GPCG;
333: tao->ops->computedual = TaoComputeDual_GPCG;
335: PetscNewLog(tao,&gpcg);
336: tao->data = (void*)gpcg;
338: /* Override default settings (unless already changed) */
339: if (!tao->max_it_changed) tao->max_it=500;
340: if (!tao->max_funcs_changed) tao->max_funcs = 100000;
341: #if defined(PETSC_USE_REAL_SINGLE)
342: if (!tao->gatol_changed) tao->gatol=1e-6;
343: if (!tao->grtol_changed) tao->grtol=1e-6;
344: #else
345: if (!tao->gatol_changed) tao->gatol=1e-12;
346: if (!tao->grtol_changed) tao->grtol=1e-12;
347: #endif
349: /* Initialize pointers and variables */
350: gpcg->n=0;
351: gpcg->maxgpits = 8;
352: gpcg->pg_ftol = 0.1;
354: gpcg->gp_iterates=0; /* Cumulative number */
355: gpcg->total_gp_its = 0;
357: /* Initialize pointers and variables */
358: gpcg->n_bind=0;
359: gpcg->n_free = 0;
360: gpcg->n_upper=0;
361: gpcg->n_lower=0;
362: gpcg->subset_type = TAO_SUBSET_MASK;
363: gpcg->Hsub=NULL;
364: gpcg->Hsub_pre=NULL;
366: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
367: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
368: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
369: KSPSetType(tao->ksp,KSPNASH);
371: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
372: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
373: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG);
374: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao);
375: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
376: return(0);
377: }