1 #include <petscksp.h>
2 #include <../src/tao/quadratic/impls/gpcg/gpcg.h> /*I "gpcg.h" I*/
3
4 static PetscErrorCode GPCGGradProjections(Tao tao);
5 static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);
6
TaoDestroy_GPCG(Tao tao)7 static PetscErrorCode TaoDestroy_GPCG(Tao tao)
8 {
9 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
10
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 }
27
TaoSetFromOptions_GPCG(Tao tao,PetscOptionItems PetscOptionsObject)28 static PetscErrorCode TaoSetFromOptions_GPCG(Tao tao, PetscOptionItems PetscOptionsObject)
29 {
30 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
31 PetscBool flg;
32
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 }
41
TaoView_GPCG(Tao tao,PetscViewer viewer)42 static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
43 {
44 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
45 PetscBool isascii;
46
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 }
56
57 /* GPCGObjectiveAndGradient()
58 Compute f=0.5 * x'Hx + b'x + c
59 g=Hx + b
60 */
GPCGObjectiveAndGradient(TaoLineSearch ls,Vec X,PetscReal * f,Vec G,void * tptr)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;
66
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 }
75
TaoSetup_GPCG(Tao tao)76 static PetscErrorCode TaoSetup_GPCG(Tao tao)
77 {
78 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
79
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));
84
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);
101
102 }
103 */
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
TaoSolve_GPCG(Tao tao)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;
114
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));
119
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));
131
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;
138
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);
144
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;
149
150 PetscCall(GPCGGradProjections(tao));
151 PetscCall(ISGetSize(gpcg->Free_Local, &gpcg->n_free));
152
153 f = gpcg->f;
154 gnorm = gpcg->gnorm;
155
156 PetscCall(KSPReset(tao->ksp));
157
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));
166
167 PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub));
168
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 }
176
177 PetscCall(KSPReset(tao->ksp));
178 PetscCall(KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre));
179
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));
186
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));
191
192 actred = f_new - f;
193
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 }
205
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 }
217
GPCGGradProjections(Tao tao)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;
227
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));
237
238 PetscCall(MatMult(tao->hessian, DX, Work));
239 PetscCall(VecDot(DX, Work, &gAg));
240
241 gpcg->gp_iterates++;
242 gpcg->total_gp_its++;
243
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));
253
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 }
261
262 gpcg->gnorm = gtg;
263 PetscFunctionReturn(PETSC_SUCCESS);
264 } /* End gradient projections */
265
TaoComputeDual_GPCG(Tao tao,Vec DXL,Vec DXU)266 static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
267 {
268 TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
269
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));
276
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 }
283
284 /*MC
285 TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
286 conjugate-gradient based method for bound-constrained minimization
287
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
291
292 Level: beginner
293 M*/
TaoCreate_GPCG(Tao tao)294 PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
295 {
296 TAO_GPCG *gpcg;
297
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;
305
306 PetscCall(PetscNew(&gpcg));
307 tao->data = (void *)gpcg;
308
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);
315
316 /* Initialize pointers and variables */
317 gpcg->n = 0;
318 gpcg->maxgpits = 8;
319 gpcg->pg_ftol = 0.1;
320
321 gpcg->gp_iterates = 0; /* Cumulative number */
322 gpcg->total_gp_its = 0;
323
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;
332
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));
337
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 }
345