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