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