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: KSPDestroy(&tao->ksp);
24: PetscFree(tao->data);
25: return 0;
26: }
28: /*------------------------------------------------------------*/
29: static PetscErrorCode TaoSetFromOptions_GPCG(Tao tao, PetscOptionItems *PetscOptionsObject)
30: {
31: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
32: PetscBool flg;
34: PetscOptionsHeadBegin(PetscOptionsObject, "Gradient Projection, Conjugate Gradient method for bound constrained optimization");
35: PetscOptionsInt("-tao_gpcg_maxpgits", "maximum number of gradient projections per GPCG iterate", NULL, gpcg->maxgpits, &gpcg->maxgpits, &flg);
36: PetscOptionsHeadEnd();
37: KSPSetFromOptions(tao->ksp);
38: TaoLineSearchSetFromOptions(tao->linesearch);
39: return 0;
40: }
42: /*------------------------------------------------------------*/
43: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
44: {
45: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
46: PetscBool isascii;
48: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
49: if (isascii) {
50: PetscViewerASCIIPrintf(viewer, "Total PG its: %" PetscInt_FMT ",", gpcg->total_gp_its);
51: PetscViewerASCIIPrintf(viewer, "PG tolerance: %g \n", (double)gpcg->pg_ftol);
52: }
53: TaoLineSearchView(tao->linesearch, viewer);
54: return 0;
55: }
57: /* GPCGObjectiveAndGradient()
58: Compute f=0.5 * x'Hx + b'x + c
59: g=Hx + b
60: */
61: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *tptr)
62: {
63: Tao tao = (Tao)tptr;
64: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
65: PetscReal f1, f2;
67: MatMult(tao->hessian, X, G);
68: VecDot(G, X, &f1);
69: VecDot(gpcg->B, X, &f2);
70: VecAXPY(G, 1.0, gpcg->B);
71: *f = f1 / 2.0 + f2 + gpcg->c;
72: return 0;
73: }
75: /* ---------------------------------------------------------- */
76: static PetscErrorCode TaoSetup_GPCG(Tao tao)
77: {
78: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
80: /* Allocate some arrays */
81: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
82: if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
84: VecDuplicate(tao->solution, &gpcg->B);
85: VecDuplicate(tao->solution, &gpcg->Work);
86: VecDuplicate(tao->solution, &gpcg->X_New);
87: VecDuplicate(tao->solution, &gpcg->G_New);
88: VecDuplicate(tao->solution, &gpcg->DXFree);
89: VecDuplicate(tao->solution, &gpcg->R);
90: VecDuplicate(tao->solution, &gpcg->PG);
91: /*
92: if (gpcg->ksp_type == GPCG_KSP_NASH) {
93: KSPSetType(tao->ksp,KSPNASH);
94: } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
95: KSPSetType(tao->ksp,KSPSTCG);
96: } else {
97: KSPSetType(tao->ksp,KSPGLTR);
98: }
99: if (tao->ksp->ops->setfromoptions) {
100: (*tao->ksp->ops->setfromoptions)(tao->ksp);
101: }
103: }
104: */
105: return 0;
106: }
108: static PetscErrorCode TaoSolve_GPCG(Tao tao)
109: {
110: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
111: PetscInt its;
112: PetscReal actred, f, f_new, gnorm, gdx, stepsize, xtb;
113: PetscReal xtHx;
114: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
117: TaoComputeVariableBounds(tao);
118: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
119: TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU);
121: /* Using f = .5*x'Hx + x'b + c and g=Hx + b, compute b,c */
122: TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
123: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
124: VecCopy(tao->gradient, gpcg->B);
125: MatMult(tao->hessian, tao->solution, gpcg->Work);
126: VecDot(gpcg->Work, tao->solution, &xtHx);
127: VecAXPY(gpcg->B, -1.0, gpcg->Work);
128: VecDot(gpcg->B, tao->solution, &xtb);
129: gpcg->c = f - xtHx / 2.0 - xtb;
130: if (gpcg->Free_Local) ISDestroy(&gpcg->Free_Local);
131: VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local);
133: /* Project the gradient and calculate the norm */
134: VecCopy(tao->gradient, gpcg->G_New);
135: VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG);
136: VecNorm(gpcg->PG, NORM_2, &gpcg->gnorm);
137: tao->step = 1.0;
138: gpcg->f = f;
140: /* Check Stopping Condition */
141: tao->reason = TAO_CONTINUE_ITERATING;
142: TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its);
143: TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step);
144: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
146: while (tao->reason == TAO_CONTINUE_ITERATING) {
147: /* Call general purpose update function */
148: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
149: tao->ksp_its = 0;
151: GPCGGradProjections(tao);
152: ISGetSize(gpcg->Free_Local, &gpcg->n_free);
154: f = gpcg->f;
155: gnorm = gpcg->gnorm;
157: KSPReset(tao->ksp);
159: if (gpcg->n_free > 0) {
160: /* Create a reduced linear system */
161: VecDestroy(&gpcg->R);
162: VecDestroy(&gpcg->DXFree);
163: TaoVecGetSubVec(tao->gradient, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R);
164: VecScale(gpcg->R, -1.0);
165: TaoVecGetSubVec(tao->stepdirection, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->DXFree);
166: VecSet(gpcg->DXFree, 0.0);
168: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub);
170: if (tao->hessian_pre == tao->hessian) {
171: MatDestroy(&gpcg->Hsub_pre);
172: PetscObjectReference((PetscObject)gpcg->Hsub);
173: gpcg->Hsub_pre = gpcg->Hsub;
174: } else {
175: TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre);
176: }
178: KSPReset(tao->ksp);
179: KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre);
181: KSPSolve(tao->ksp, gpcg->R, gpcg->DXFree);
182: KSPGetIterationNumber(tao->ksp, &its);
183: tao->ksp_its += its;
184: tao->ksp_tot_its += its;
185: VecSet(tao->stepdirection, 0.0);
186: VecISAXPY(tao->stepdirection, gpcg->Free_Local, 1.0, gpcg->DXFree);
188: VecDot(tao->stepdirection, tao->gradient, &gdx);
189: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
190: f_new = f;
191: TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &stepsize, &ls_status);
193: actred = f_new - f;
195: /* Evaluate the function and gradient at the new point */
196: VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG);
197: VecNorm(gpcg->PG, NORM_2, &gnorm);
198: f = f_new;
199: ISDestroy(&gpcg->Free_Local);
200: VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local);
201: } else {
202: actred = 0;
203: gpcg->step = 1.0;
204: /* if there were no free variables, no cg method */
205: }
207: tao->niter++;
208: gpcg->f = f;
209: gpcg->gnorm = gnorm;
210: gpcg->actred = actred;
211: TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its);
212: TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step);
213: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
214: if (tao->reason != TAO_CONTINUE_ITERATING) break;
215: } /* END MAIN LOOP */
217: return 0;
218: }
220: static PetscErrorCode GPCGGradProjections(Tao tao)
221: {
222: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
223: PetscInt i;
224: PetscReal actred = -1.0, actred_max = 0.0, gAg, gtg = gpcg->gnorm, alpha;
225: PetscReal f_new, gdx, stepsize;
226: Vec DX = tao->stepdirection, XL = tao->XL, XU = tao->XU, Work = gpcg->Work;
227: Vec X = tao->solution, G = tao->gradient;
228: TaoLineSearchConvergedReason lsflag = TAOLINESEARCH_CONTINUE_ITERATING;
230: /*
231: The free, active, and binding variables should be already identified
232: */
233: for (i = 0; i < gpcg->maxgpits; i++) {
234: if (-actred <= (gpcg->pg_ftol) * actred_max) break;
235: VecBoundGradientProjection(G, X, XL, XU, DX);
236: VecScale(DX, -1.0);
237: VecDot(DX, G, &gdx);
239: MatMult(tao->hessian, DX, Work);
240: VecDot(DX, Work, &gAg);
242: gpcg->gp_iterates++;
243: gpcg->total_gp_its++;
245: gtg = -gdx;
246: if (PetscAbsReal(gAg) == 0.0) {
247: alpha = 1.0;
248: } else {
249: alpha = PetscAbsReal(gtg / gAg);
250: }
251: TaoLineSearchSetInitialStepLength(tao->linesearch, alpha);
252: f_new = gpcg->f;
253: TaoLineSearchApply(tao->linesearch, X, &f_new, G, DX, &stepsize, &lsflag);
255: /* Update the iterate */
256: actred = f_new - gpcg->f;
257: actred_max = PetscMax(actred_max, -(f_new - gpcg->f));
258: gpcg->f = f_new;
259: ISDestroy(&gpcg->Free_Local);
260: VecWhichInactive(XL, X, tao->gradient, XU, PETSC_TRUE, &gpcg->Free_Local);
261: }
263: gpcg->gnorm = gtg;
264: return 0;
265: } /* End gradient projections */
267: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
268: {
269: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
271: VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work);
272: VecCopy(gpcg->Work, DXL);
273: VecAXPY(DXL, -1.0, tao->gradient);
274: VecSet(DXU, 0.0);
275: VecPointwiseMax(DXL, DXL, DXU);
277: VecCopy(tao->gradient, DXU);
278: VecAXPY(DXU, -1.0, gpcg->Work);
279: VecSet(gpcg->Work, 0.0);
280: VecPointwiseMin(DXU, gpcg->Work, DXU);
281: return 0;
282: }
284: /*------------------------------------------------------------*/
285: /*MC
286: TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
287: conjugate-gradient based method for bound-constrained minimization
289: Options Database Keys:
290: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
291: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
293: Level: beginner
294: M*/
295: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
296: {
297: TAO_GPCG *gpcg;
299: tao->ops->setup = TaoSetup_GPCG;
300: tao->ops->solve = TaoSolve_GPCG;
301: tao->ops->view = TaoView_GPCG;
302: tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
303: tao->ops->destroy = TaoDestroy_GPCG;
304: tao->ops->computedual = TaoComputeDual_GPCG;
306: PetscNew(&gpcg);
307: tao->data = (void *)gpcg;
309: /* Override default settings (unless already changed) */
310: if (!tao->max_it_changed) tao->max_it = 500;
311: if (!tao->max_funcs_changed) tao->max_funcs = 100000;
312: #if defined(PETSC_USE_REAL_SINGLE)
313: if (!tao->gatol_changed) tao->gatol = 1e-6;
314: if (!tao->grtol_changed) tao->grtol = 1e-6;
315: #else
316: if (!tao->gatol_changed) tao->gatol = 1e-12;
317: if (!tao->grtol_changed) tao->grtol = 1e-12;
318: #endif
320: /* Initialize pointers and variables */
321: gpcg->n = 0;
322: gpcg->maxgpits = 8;
323: gpcg->pg_ftol = 0.1;
325: gpcg->gp_iterates = 0; /* Cumulative number */
326: gpcg->total_gp_its = 0;
328: /* Initialize pointers and variables */
329: gpcg->n_bind = 0;
330: gpcg->n_free = 0;
331: gpcg->n_upper = 0;
332: gpcg->n_lower = 0;
333: gpcg->subset_type = TAO_SUBSET_MASK;
334: gpcg->Hsub = NULL;
335: gpcg->Hsub_pre = NULL;
337: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
338: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
339: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
340: KSPSetType(tao->ksp, KSPNASH);
342: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
343: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
344: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG);
345: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao);
346: TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
347: return 0;
348: }