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: static PetscErrorCode TaoDestroy_GPCG(Tao tao)
8: {
9: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
11: /* Free allocated memory in GPCG structure */
12: PetscFunctionBegin;
13: PetscCall(VecDestroy(&gpcg->B));
14: PetscCall(VecDestroy(&gpcg->Work));
15: PetscCall(VecDestroy(&gpcg->X_New));
16: PetscCall(VecDestroy(&gpcg->G_New));
17: PetscCall(VecDestroy(&gpcg->DXFree));
18: PetscCall(VecDestroy(&gpcg->R));
19: PetscCall(VecDestroy(&gpcg->PG));
20: PetscCall(MatDestroy(&gpcg->Hsub));
21: PetscCall(MatDestroy(&gpcg->Hsub_pre));
22: PetscCall(ISDestroy(&gpcg->Free_Local));
23: PetscCall(KSPDestroy(&tao->ksp));
24: PetscCall(PetscFree(tao->data));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode TaoSetFromOptions_GPCG(Tao tao, PetscOptionItems PetscOptionsObject)
29: {
30: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
31: PetscBool flg;
33: PetscFunctionBegin;
34: PetscOptionsHeadBegin(PetscOptionsObject, "Gradient Projection, Conjugate Gradient method for bound constrained optimization");
35: PetscCall(PetscOptionsInt("-tao_gpcg_maxpgits", "maximum number of gradient projections per GPCG iterate", NULL, gpcg->maxgpits, &gpcg->maxgpits, &flg));
36: PetscOptionsHeadEnd();
37: PetscCall(KSPSetFromOptions(tao->ksp));
38: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
43: {
44: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
45: PetscBool isascii;
47: PetscFunctionBegin;
48: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
49: if (isascii) {
50: PetscCall(PetscViewerASCIIPrintf(viewer, "Total PG its: %" PetscInt_FMT ",", gpcg->total_gp_its));
51: PetscCall(PetscViewerASCIIPrintf(viewer, "PG tolerance: %g \n", (double)gpcg->pg_ftol));
52: }
53: PetscCall(TaoLineSearchView(tao->linesearch, viewer));
54: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
68: PetscCall(MatMult(tao->hessian, X, G));
69: PetscCall(VecDot(G, X, &f1));
70: PetscCall(VecDot(gpcg->B, X, &f2));
71: PetscCall(VecAXPY(G, 1.0, gpcg->B));
72: *f = f1 / 2.0 + f2 + gpcg->c;
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode TaoSetup_GPCG(Tao tao)
77: {
78: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
80: PetscFunctionBegin;
81: /* Allocate some arrays */
82: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
83: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
85: PetscCall(VecDuplicate(tao->solution, &gpcg->B));
86: PetscCall(VecDuplicate(tao->solution, &gpcg->Work));
87: PetscCall(VecDuplicate(tao->solution, &gpcg->X_New));
88: PetscCall(VecDuplicate(tao->solution, &gpcg->G_New));
89: PetscCall(VecDuplicate(tao->solution, &gpcg->DXFree));
90: PetscCall(VecDuplicate(tao->solution, &gpcg->R));
91: PetscCall(VecDuplicate(tao->solution, &gpcg->PG));
92: /*
93: if (gpcg->ksp_type == GPCG_KSP_NASH) {
94: PetscCall(KSPSetType(tao->ksp,KSPNASH));
95: } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
96: PetscCall(KSPSetType(tao->ksp,KSPSTCG));
97: } else {
98: PetscCall(KSPSetType(tao->ksp,KSPGLTR));
99: }
100: if (tao->ksp->ops->setfromoptions) (*tao->ksp->ops->setfromoptions)(tao->ksp);
102: }
103: */
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode TaoSolve_GPCG(Tao tao)
108: {
109: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
110: PetscInt its;
111: PetscReal actred, f, f_new, gnorm, gdx, stepsize, xtb;
112: PetscReal xtHx;
113: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
115: PetscFunctionBegin;
116: PetscCall(TaoComputeVariableBounds(tao));
117: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
118: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
120: /* Using f = .5*x'Hx + x'b + c and g=Hx + b, compute b,c */
121: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
122: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
123: PetscCall(VecCopy(tao->gradient, gpcg->B));
124: PetscCall(MatMult(tao->hessian, tao->solution, gpcg->Work));
125: PetscCall(VecDot(gpcg->Work, tao->solution, &xtHx));
126: PetscCall(VecAXPY(gpcg->B, -1.0, gpcg->Work));
127: PetscCall(VecDot(gpcg->B, tao->solution, &xtb));
128: gpcg->c = f - xtHx / 2.0 - xtb;
129: if (gpcg->Free_Local) PetscCall(ISDestroy(&gpcg->Free_Local));
130: PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));
132: /* Project the gradient and calculate the norm */
133: PetscCall(VecCopy(tao->gradient, gpcg->G_New));
134: PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
135: PetscCall(VecNorm(gpcg->PG, NORM_2, &gpcg->gnorm));
136: tao->step = 1.0;
137: gpcg->f = f;
139: /* Check Stopping Condition */
140: tao->reason = TAO_CONTINUE_ITERATING;
141: PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
142: PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
143: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
145: while (tao->reason == TAO_CONTINUE_ITERATING) {
146: /* Call general purpose update function */
147: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
148: tao->ksp_its = 0;
150: PetscCall(GPCGGradProjections(tao));
151: PetscCall(ISGetSize(gpcg->Free_Local, &gpcg->n_free));
153: f = gpcg->f;
154: gnorm = gpcg->gnorm;
156: PetscCall(KSPReset(tao->ksp));
158: if (gpcg->n_free > 0) {
159: /* Create a reduced linear system */
160: PetscCall(VecDestroy(&gpcg->R));
161: PetscCall(VecDestroy(&gpcg->DXFree));
162: PetscCall(TaoVecGetSubVec(tao->gradient, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R));
163: PetscCall(VecScale(gpcg->R, -1.0));
164: PetscCall(TaoVecGetSubVec(tao->stepdirection, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->DXFree));
165: PetscCall(VecSet(gpcg->DXFree, 0.0));
167: PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub));
169: if (tao->hessian_pre == tao->hessian) {
170: PetscCall(MatDestroy(&gpcg->Hsub_pre));
171: PetscCall(PetscObjectReference((PetscObject)gpcg->Hsub));
172: gpcg->Hsub_pre = gpcg->Hsub;
173: } else {
174: PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre));
175: }
177: PetscCall(KSPReset(tao->ksp));
178: PetscCall(KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre));
180: PetscCall(KSPSolve(tao->ksp, gpcg->R, gpcg->DXFree));
181: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
182: tao->ksp_its += its;
183: tao->ksp_tot_its += its;
184: PetscCall(VecSet(tao->stepdirection, 0.0));
185: PetscCall(VecISAXPY(tao->stepdirection, gpcg->Free_Local, 1.0, gpcg->DXFree));
187: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
188: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
189: f_new = f;
190: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &stepsize, &ls_status));
192: actred = f_new - f;
194: /* Evaluate the function and gradient at the new point */
195: PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
196: PetscCall(VecNorm(gpcg->PG, NORM_2, &gnorm));
197: f = f_new;
198: PetscCall(ISDestroy(&gpcg->Free_Local));
199: PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));
200: } else {
201: actred = 0;
202: gpcg->step = 1.0;
203: /* if there were no free variables, no cg method */
204: }
206: tao->niter++;
207: gpcg->f = f;
208: gpcg->gnorm = gnorm;
209: gpcg->actred = actred;
210: PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
211: PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
212: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
213: if (tao->reason != TAO_CONTINUE_ITERATING) break;
214: } /* END MAIN LOOP */
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode GPCGGradProjections(Tao tao)
219: {
220: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
221: PetscInt i;
222: PetscReal actred = -1.0, actred_max = 0.0, gAg, gtg = gpcg->gnorm, alpha;
223: PetscReal f_new, gdx, stepsize;
224: Vec DX = tao->stepdirection, XL = tao->XL, XU = tao->XU, Work = gpcg->Work;
225: Vec X = tao->solution, G = tao->gradient;
226: TaoLineSearchConvergedReason lsflag = TAOLINESEARCH_CONTINUE_ITERATING;
228: /*
229: The free, active, and binding variables should be already identified
230: */
231: PetscFunctionBegin;
232: for (i = 0; i < gpcg->maxgpits; i++) {
233: if (-actred <= (gpcg->pg_ftol) * actred_max) break;
234: PetscCall(VecBoundGradientProjection(G, X, XL, XU, DX));
235: PetscCall(VecScale(DX, -1.0));
236: PetscCall(VecDot(DX, G, &gdx));
238: PetscCall(MatMult(tao->hessian, DX, Work));
239: PetscCall(VecDot(DX, Work, &gAg));
241: gpcg->gp_iterates++;
242: gpcg->total_gp_its++;
244: gtg = -gdx;
245: if (PetscAbsReal(gAg) == 0.0) {
246: alpha = 1.0;
247: } else {
248: alpha = PetscAbsReal(gtg / gAg);
249: }
250: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, alpha));
251: f_new = gpcg->f;
252: PetscCall(TaoLineSearchApply(tao->linesearch, X, &f_new, G, DX, &stepsize, &lsflag));
254: /* Update the iterate */
255: actred = f_new - gpcg->f;
256: actred_max = PetscMax(actred_max, -(f_new - gpcg->f));
257: gpcg->f = f_new;
258: PetscCall(ISDestroy(&gpcg->Free_Local));
259: PetscCall(VecWhichInactive(XL, X, tao->gradient, XU, PETSC_TRUE, &gpcg->Free_Local));
260: }
262: gpcg->gnorm = gtg;
263: PetscFunctionReturn(PETSC_SUCCESS);
264: } /* End gradient projections */
266: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
267: {
268: TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
270: PetscFunctionBegin;
271: PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work));
272: PetscCall(VecCopy(gpcg->Work, DXL));
273: PetscCall(VecAXPY(DXL, -1.0, tao->gradient));
274: PetscCall(VecSet(DXU, 0.0));
275: PetscCall(VecPointwiseMax(DXL, DXL, DXU));
277: PetscCall(VecCopy(tao->gradient, DXU));
278: PetscCall(VecAXPY(DXU, -1.0, gpcg->Work));
279: PetscCall(VecSet(gpcg->Work, 0.0));
280: PetscCall(VecPointwiseMin(DXU, gpcg->Work, DXU));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*MC
285: TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
286: conjugate-gradient based method for bound-constrained minimization
288: Options Database Keys:
289: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
290: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets
292: Level: beginner
293: M*/
294: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
295: {
296: TAO_GPCG *gpcg;
298: PetscFunctionBegin;
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: PetscCall(PetscNew(&gpcg));
307: tao->data = (void *)gpcg;
309: /* Override default settings (unless already changed) */
310: PetscCall(TaoParametersInitialize(tao));
311: PetscObjectParameterSetDefault(tao, max_it, 500);
312: PetscObjectParameterSetDefault(tao, max_funcs, 100000);
313: PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
314: PetscObjectParameterSetDefault(tao, grtol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
316: /* Initialize pointers and variables */
317: gpcg->n = 0;
318: gpcg->maxgpits = 8;
319: gpcg->pg_ftol = 0.1;
321: gpcg->gp_iterates = 0; /* Cumulative number */
322: gpcg->total_gp_its = 0;
324: /* Initialize pointers and variables */
325: gpcg->n_bind = 0;
326: gpcg->n_free = 0;
327: gpcg->n_upper = 0;
328: gpcg->n_lower = 0;
329: gpcg->subset_type = TAO_SUBSET_MASK;
330: gpcg->Hsub = NULL;
331: gpcg->Hsub_pre = NULL;
333: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
334: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
335: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
336: PetscCall(KSPSetType(tao->ksp, KSPNASH));
338: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
339: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
340: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG));
341: PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao));
342: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }