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