Actual source code: gpcg.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/kspimpl.h>
2: #include <../src/tao/bound/impls/gpcg/gpcg.h> /*I "gpcg.h" I*/
5: static PetscErrorCode GPCGGradProjections(Tao tao);
6: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
8: /*------------------------------------------------------------*/
11: static PetscErrorCode TaoDestroy_GPCG(Tao tao)
12: {
13: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
16: /* Free allocated memory in GPCG structure */
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 = NULL;
30: return(0);
31: }
33: /*------------------------------------------------------------*/
36: static PetscErrorCode TaoSetFromOptions_GPCG(PetscOptions *PetscOptionsObject,Tao tao)
37: {
38: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
40: PetscBool flg;
43: PetscOptionsHead(PetscOptionsObject,"Gradient Projection, Conjugate Gradient method for bound constrained optimization");
44: ierr=PetscOptionsInt("-tao_gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",NULL,gpcg->maxgpits,&gpcg->maxgpits,&flg);
45: PetscOptionsTail();
46: KSPSetFromOptions(tao->ksp);
47: TaoLineSearchSetFromOptions(tao->linesearch);
48: return(0);
49: }
51: /*------------------------------------------------------------*/
54: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
55: {
56: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
57: PetscBool isascii;
61: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
62: if (isascii) {
63: PetscViewerASCIIPushTab(viewer);
64: PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",gpcg->total_gp_its);
65: PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)gpcg->pg_ftol);
66: PetscViewerASCIIPopTab(viewer);
67: }
68: TaoLineSearchView(tao->linesearch,viewer);
69: return(0);
70: }
72: /* GPCGObjectiveAndGradient()
73: Compute f=0.5 * x'Hx + b'x + c
74: g=Hx + b
75: */
78: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void*tptr)
79: {
80: Tao tao = (Tao)tptr;
81: TAO_GPCG *gpcg = (TAO_GPCG*)tao->data;
83: PetscReal f1,f2;
86: MatMult(tao->hessian,X,G);
87: VecDot(G,X,&f1);
88: VecDot(gpcg->B,X,&f2);
89: VecAXPY(G,1.0,gpcg->B);
90: *f=f1/2.0 + f2 + gpcg->c;
91: return(0);
92: }
94: /* ---------------------------------------------------------- */
97: static PetscErrorCode TaoSetup_GPCG(Tao tao)
98: {
100: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
103: /* Allocate some arrays */
104: if (!tao->gradient) {
105: VecDuplicate(tao->solution, &tao->gradient);
106: }
107: if (!tao->stepdirection) {
108: VecDuplicate(tao->solution, &tao->stepdirection);
109: }
110: if (!tao->XL) {
111: VecDuplicate(tao->solution,&tao->XL);
112: VecSet(tao->XL,PETSC_NINFINITY);
113: }
114: if (!tao->XU) {
115: VecDuplicate(tao->solution,&tao->XU);
116: VecSet(tao->XU,PETSC_INFINITY);
117: }
119: ierr=VecDuplicate(tao->solution,&gpcg->B);
120: ierr=VecDuplicate(tao->solution,&gpcg->Work);
121: ierr=VecDuplicate(tao->solution,&gpcg->X_New);
122: ierr=VecDuplicate(tao->solution,&gpcg->G_New);
123: ierr=VecDuplicate(tao->solution,&gpcg->DXFree);
124: ierr=VecDuplicate(tao->solution,&gpcg->R);
125: ierr=VecDuplicate(tao->solution,&gpcg->PG);
126: /*
127: if (gpcg->ksp_type == GPCG_KSP_NASH) {
128: KSPSetType(tao->ksp,KSPNASH);
129: } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
130: KSPSetType(tao->ksp,KSPSTCG);
131: } else {
132: KSPSetType(tao->ksp,KSPGLTR);
133: }
134: if (tao->ksp->ops->setfromoptions) {
135: (*tao->ksp->ops->setfromoptions)(tao->ksp);
136: }
138: }
139: */
140: return(0);
141: }
145: static PetscErrorCode TaoSolve_GPCG(Tao tao)
146: {
147: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
148: PetscErrorCode ierr;
149: PetscInt its;
150: PetscReal actred,f,f_new,gnorm,gdx,stepsize,xtb;
151: PetscReal xtHx;
152: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
153: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
156: gpcg->Hsub=NULL;
157: gpcg->Hsub_pre=NULL;
159: TaoComputeVariableBounds(tao);
160: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
161: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
163: /* Using f = .5*x'Hx + x'b + c and g=Hx + b, compute b,c */
164: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
165: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
166: VecCopy(tao->gradient, gpcg->B);
167: MatMult(tao->hessian,tao->solution,gpcg->Work);
168: VecDot(gpcg->Work, tao->solution, &xtHx);
169: VecAXPY(gpcg->B,-1.0,gpcg->Work);
170: VecDot(gpcg->B,tao->solution,&xtb);
171: gpcg->c=f-xtHx/2.0-xtb;
172: if (gpcg->Free_Local) {
173: ISDestroy(&gpcg->Free_Local);
174: }
175: VecWhichBetween(tao->XL,tao->solution,tao->XU,&gpcg->Free_Local);
177: /* Project the gradient and calculate the norm */
178: VecCopy(tao->gradient,gpcg->G_New);
179: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,gpcg->PG);
180: VecNorm(gpcg->PG,NORM_2,&gpcg->gnorm);
181: tao->step=1.0;
182: gpcg->f = f;
184: /* Check Stopping Condition */
185: ierr=TaoMonitor(tao,tao->niter,f,gpcg->gnorm,0.0,tao->step,&reason);
187: while (reason == TAO_CONTINUE_ITERATING){
188: tao->ksp_its=0;
190: GPCGGradProjections(tao);
191: ISGetSize(gpcg->Free_Local,&gpcg->n_free);
193: f=gpcg->f; gnorm=gpcg->gnorm;
195: KSPReset(tao->ksp);
197: if (gpcg->n_free > 0){
198: /* Create a reduced linear system */
199: VecDestroy(&gpcg->R);
200: VecDestroy(&gpcg->DXFree);
201: TaoVecGetSubVec(tao->gradient,gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R);
202: VecScale(gpcg->R, -1.0);
203: TaoVecGetSubVec(tao->stepdirection,gpcg->Free_Local,tao->subset_type, 0.0, &gpcg->DXFree);
204: VecSet(gpcg->DXFree,0.0);
206: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub);
208: if (tao->hessian_pre == tao->hessian) {
209: MatDestroy(&gpcg->Hsub_pre);
210: PetscObjectReference((PetscObject)gpcg->Hsub);
211: gpcg->Hsub_pre = gpcg->Hsub;
212: } else {
213: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre);
214: }
216: KSPReset(tao->ksp);
217: KSPSetOperators(tao->ksp,gpcg->Hsub,gpcg->Hsub_pre);
219: KSPSolve(tao->ksp,gpcg->R,gpcg->DXFree);
220: KSPGetIterationNumber(tao->ksp,&its);
221: tao->ksp_its+=its;
222: tao->ksp_tot_its+=its;
223: VecSet(tao->stepdirection,0.0);
224: VecISAXPY(tao->stepdirection,gpcg->Free_Local,1.0,gpcg->DXFree);
226: VecDot(tao->stepdirection,tao->gradient,&gdx);
227: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
228: f_new=f;
229: TaoLineSearchApply(tao->linesearch,tao->solution,&f_new,tao->gradient,tao->stepdirection,&stepsize,&ls_status);
231: actred = f_new - f;
233: /* Evaluate the function and gradient at the new point */
234: VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU, gpcg->PG);
235: VecNorm(gpcg->PG, NORM_2, &gnorm);
236: f=f_new;
237: ISDestroy(&gpcg->Free_Local);
238: VecWhichBetween(tao->XL,tao->solution,tao->XU,&gpcg->Free_Local);
239: } else {
240: actred = 0; gpcg->step=1.0;
241: /* if there were no free variables, no cg method */
242: }
244: tao->niter++;
245: TaoMonitor(tao,tao->niter,f,gnorm,0.0,gpcg->step,&reason);
246: gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
247: if (reason!=TAO_CONTINUE_ITERATING) break;
248: } /* END MAIN LOOP */
250: return(0);
251: }
255: static PetscErrorCode GPCGGradProjections(Tao tao)
256: {
257: PetscErrorCode ierr;
258: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
259: PetscInt i;
260: PetscReal actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
261: PetscReal f_new,gdx,stepsize;
262: Vec DX=tao->stepdirection,XL=tao->XL,XU=tao->XU,Work=gpcg->Work;
263: Vec X=tao->solution,G=tao->gradient;
264: TaoLineSearchConvergedReason lsflag=TAOLINESEARCH_CONTINUE_ITERATING;
266: /*
267: The free, active, and binding variables should be already identified
268: */
270: for (i=0;i<gpcg->maxgpits;i++){
271: if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
272: VecBoundGradientProjection(G,X,XL,XU,DX);
273: VecScale(DX,-1.0);
274: VecDot(DX,G,&gdx);
276: MatMult(tao->hessian,DX,Work);
277: VecDot(DX,Work,&gAg);
279: gpcg->gp_iterates++;
280: gpcg->total_gp_its++;
282: gtg=-gdx;
283: alpha = PetscAbsReal(gtg/gAg);
284: TaoLineSearchSetInitialStepLength(tao->linesearch,alpha);
285: f_new=gpcg->f;
286: TaoLineSearchApply(tao->linesearch,X,&f_new,G,DX,&stepsize,&lsflag);
288: /* Update the iterate */
289: actred = f_new - gpcg->f;
290: actred_max = PetscMax(actred_max,-(f_new - gpcg->f));
291: gpcg->f = f_new;
292: ISDestroy(&gpcg->Free_Local);
293: VecWhichBetween(XL,X,XU,&gpcg->Free_Local);
294: }
296: gpcg->gnorm=gtg;
297: return(0);
298: } /* End gradient projections */
302: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
303: {
304: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
308: VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work);
309: VecCopy(gpcg->Work, DXL);
310: VecAXPY(DXL,-1.0,tao->gradient);
311: VecSet(DXU,0.0);
312: VecPointwiseMax(DXL,DXL,DXU);
314: VecCopy(tao->gradient,DXU);
315: VecAXPY(DXU,-1.0,gpcg->Work);
316: VecSet(gpcg->Work,0.0);
317: VecPointwiseMin(DXU,gpcg->Work,DXU);
318: return(0);
319: }
321: /*------------------------------------------------------------*/
322: /*MC
323: TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
324: conjugate-gradient based method for bound-constrained minimization
326: Options Database Keys:
327: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
328: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
330: Level: beginner
331: M*/
334: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
335: {
336: TAO_GPCG *gpcg;
340: tao->ops->setup = TaoSetup_GPCG;
341: tao->ops->solve = TaoSolve_GPCG;
342: tao->ops->view = TaoView_GPCG;
343: tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
344: tao->ops->destroy = TaoDestroy_GPCG;
345: tao->ops->computedual = TaoComputeDual_GPCG;
347: PetscNewLog(tao,&gpcg);
348: tao->data = (void*)gpcg;
350: /* Override default settings (unless already changed) */
351: if (!tao->max_it_changed) tao->max_it=500;
352: if (!tao->max_funcs_changed) tao->max_funcs = 100000;
353: #if defined(PETSC_USE_REAL_SINGLE)
354: if (!tao->fatol_changed) tao->fatol=1e-6;
355: if (!tao->frtol_changed) tao->frtol=1e-6;
356: if (!tao->gatol_changed) tao->grtol=1e-6;
357: if (!tao->grtol_changed) tao->grtol=1e-6;
358: #else
359: if (!tao->fatol_changed) tao->fatol=1e-12;
360: if (!tao->frtol_changed) tao->frtol=1e-12;
361: if (!tao->gatol_changed) tao->grtol=1e-12;
362: if (!tao->grtol_changed) tao->grtol=1e-12;
363: #endif
365: /* Initialize pointers and variables */
366: gpcg->n=0;
367: gpcg->maxgpits = 8;
368: gpcg->pg_ftol = 0.1;
370: gpcg->gp_iterates=0; /* Cumulative number */
371: gpcg->total_gp_its = 0;
373: /* Initialize pointers and variables */
374: gpcg->n_bind=0;
375: gpcg->n_free = 0;
376: gpcg->n_upper=0;
377: gpcg->n_lower=0;
378: gpcg->subset_type = TAO_SUBSET_MASK;
379: /* gpcg->ksp_type = GPCG_KSP_STCG; */
381: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
382: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
383: KSPSetType(tao->ksp,KSPNASH);
385: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
386: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG);
387: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao);
388: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
389: return(0);
390: }