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