Actual source code: almm.c
1: #include <../src/tao/constrained/impls/almm/almm.h>
2: #include <petsctao.h>
3: #include <petsc/private/petscimpl.h>
4: #include <petsc/private/vecimpl.h>
6: static PetscErrorCode TaoALMMCombinePrimal_Private(Tao, Vec, Vec, Vec);
7: static PetscErrorCode TaoALMMCombineDual_Private(Tao, Vec, Vec, Vec);
8: static PetscErrorCode TaoALMMSplitPrimal_Private(Tao, Vec, Vec, Vec);
9: static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao);
10: static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao);
11: static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao);
13: static PetscErrorCode TaoSolve_ALMM(Tao tao)
14: {
15: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
16: TaoConvergedReason reason;
17: PetscReal updated;
19: /* reset initial multiplier/slack guess */
20: if (!tao->recycle) {
21: if (tao->ineq_constrained) {
22: VecZeroEntries(auglag->Ps);
23: TaoALMMCombinePrimal_Private(tao, auglag->Px, auglag->Ps, auglag->P);
24: VecSet(auglag->Yi, 1.0);
25: }
26: if (tao->eq_constrained) VecSet(auglag->Ye, 1.0);
27: }
29: /* compute initial nonlinear Lagrangian and its derivatives */
30: (*auglag->sub_obj)(tao);
31: TaoALMMComputeOptimalityNorms_Private(tao);
32: /* print initial step and check convergence */
33: PetscInfo(tao, "Solving with %s formulation\n", TaoALMMTypes[auglag->type]);
34: TaoLogConvergenceHistory(tao, auglag->Lval, auglag->gnorm, auglag->cnorm, tao->ksp_its);
35: TaoMonitor(tao, tao->niter, auglag->fval, auglag->gnorm, auglag->cnorm, 0.0);
36: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
37: /* set initial penalty factor and inner solver tolerance */
38: switch (auglag->type) {
39: case TAO_ALMM_CLASSIC:
40: auglag->mu = auglag->mu0;
41: break;
42: case TAO_ALMM_PHR:
43: auglag->cenorm = 0.0;
44: if (tao->eq_constrained) VecDot(auglag->Ce, auglag->Ce, &auglag->cenorm);
45: auglag->cinorm = 0.0;
46: if (tao->ineq_constrained) {
47: VecCopy(auglag->Ci, auglag->Ciwork);
48: VecScale(auglag->Ciwork, -1.0);
49: VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork);
50: VecDot(auglag->Ciwork, auglag->Ciwork, &auglag->cinorm);
51: }
52: /* determine initial penalty factor based on the balance of constraint violation and objective function value */
53: auglag->mu = PetscMax(1.e-6, PetscMin(10.0, 2.0 * PetscAbsReal(auglag->fval) / (auglag->cenorm + auglag->cinorm)));
54: break;
55: default:
56: break;
57: }
58: auglag->gtol = auglag->gtol0;
59: PetscInfo(tao, "Initial penalty: %.2g\n", (double)auglag->mu);
61: /* start aug-lag outer loop */
62: while (tao->reason == TAO_CONTINUE_ITERATING) {
63: ++tao->niter;
64: /* update subsolver tolerance */
65: PetscInfo(tao, "Subsolver tolerance: ||G|| <= %e\n", (double)auglag->gtol);
66: TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0);
67: /* solve the bound-constrained or unconstrained subproblem */
68: VecCopy(auglag->P, auglag->subsolver->solution);
69: TaoSolve(auglag->subsolver);
70: VecCopy(auglag->subsolver->solution, auglag->P);
71: TaoGetConvergedReason(auglag->subsolver, &reason);
72: tao->ksp_its += auglag->subsolver->ksp_its;
73: if (reason != TAO_CONVERGED_GATOL) PetscInfo(tao, "Subsolver failed to converge, reason: %s\n", TaoConvergedReasons[reason]);
74: /* evaluate solution and test convergence */
75: (*auglag->sub_obj)(tao);
76: TaoALMMComputeOptimalityNorms_Private(tao);
77: /* decide whether to update multipliers or not */
78: updated = 0.0;
79: if (auglag->cnorm <= auglag->ytol) {
80: PetscInfo(tao, "Multipliers updated: ||C|| <= %e\n", (double)auglag->ytol);
81: /* constraints are good, update multipliers and convergence tolerances */
82: if (tao->eq_constrained) {
83: VecAXPY(auglag->Ye, auglag->mu, auglag->Ce);
84: VecSet(auglag->Cework, auglag->ye_max);
85: VecPointwiseMin(auglag->Ye, auglag->Cework, auglag->Ye);
86: VecSet(auglag->Cework, auglag->ye_min);
87: VecPointwiseMax(auglag->Ye, auglag->Cework, auglag->Ye);
88: }
89: if (tao->ineq_constrained) {
90: VecAXPY(auglag->Yi, auglag->mu, auglag->Ci);
91: VecSet(auglag->Ciwork, auglag->yi_max);
92: VecPointwiseMin(auglag->Yi, auglag->Ciwork, auglag->Yi);
93: VecSet(auglag->Ciwork, auglag->yi_min);
94: VecPointwiseMax(auglag->Yi, auglag->Ciwork, auglag->Yi);
95: }
96: /* tolerances are updated only for non-PHR methods */
97: if (auglag->type != TAO_ALMM_PHR) {
98: auglag->ytol = PetscMax(tao->catol, auglag->ytol / PetscPowReal(auglag->mu, auglag->mu_pow_good));
99: auglag->gtol = PetscMax(tao->gatol, auglag->gtol / auglag->mu);
100: }
101: updated = 1.0;
102: } else {
103: /* constraints are bad, update penalty factor */
104: auglag->mu = PetscMin(auglag->mu_max, auglag->mu_fac * auglag->mu);
105: /* tolerances are reset only for non-PHR methods */
106: if (auglag->type != TAO_ALMM_PHR) {
107: auglag->ytol = PetscMax(tao->catol, 0.1 / PetscPowReal(auglag->mu, auglag->mu_pow_bad));
108: auglag->gtol = PetscMax(tao->gatol, 1.0 / auglag->mu);
109: }
110: PetscInfo(tao, "Penalty increased: mu = %.2g\n", (double)auglag->mu);
111: }
112: TaoLogConvergenceHistory(tao, auglag->fval, auglag->gnorm, auglag->cnorm, tao->ksp_its);
113: TaoMonitor(tao, tao->niter, auglag->fval, auglag->gnorm, auglag->cnorm, updated);
114: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
115: }
117: return 0;
118: }
120: static PetscErrorCode TaoView_ALMM(Tao tao, PetscViewer viewer)
121: {
122: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
123: PetscBool isascii;
125: PetscViewerASCIIPushTab(viewer);
126: TaoView(auglag->subsolver, viewer);
127: PetscViewerASCIIPopTab(viewer);
128: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
129: if (isascii) {
130: PetscViewerASCIIPushTab(viewer);
131: PetscViewerASCIIPrintf(viewer, "ALMM Formulation Type: %s\n", TaoALMMTypes[auglag->type]);
132: PetscViewerASCIIPopTab(viewer);
133: }
134: return 0;
135: }
137: static PetscErrorCode TaoSetUp_ALMM(Tao tao)
138: {
139: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
140: VecType vec_type;
141: Vec SL, SU;
142: PetscBool is_cg = PETSC_FALSE, is_lmvm = PETSC_FALSE;
146: TaoComputeVariableBounds(tao);
147: /* alias base vectors and create extras */
148: VecGetType(tao->solution, &vec_type);
149: auglag->Px = tao->solution;
150: if (!tao->gradient) { /* base gradient */
151: VecDuplicate(tao->solution, &tao->gradient);
152: }
153: auglag->LgradX = tao->gradient;
154: if (!auglag->Xwork) { /* opt var work vector */
155: VecDuplicate(tao->solution, &auglag->Xwork);
156: }
157: if (tao->eq_constrained) {
158: auglag->Ce = tao->constraints_equality;
159: auglag->Ae = tao->jacobian_equality;
160: if (!auglag->Ye) { /* equality multipliers */
161: VecDuplicate(auglag->Ce, &auglag->Ye);
162: }
163: if (!auglag->Cework) VecDuplicate(auglag->Ce, &auglag->Cework);
164: }
165: if (tao->ineq_constrained) {
166: auglag->Ci = tao->constraints_inequality;
167: auglag->Ai = tao->jacobian_inequality;
168: if (!auglag->Yi) { /* inequality multipliers */
169: VecDuplicate(auglag->Ci, &auglag->Yi);
170: }
171: if (!auglag->Ciwork) VecDuplicate(auglag->Ci, &auglag->Ciwork);
172: if (!auglag->Cizero) {
173: VecDuplicate(auglag->Ci, &auglag->Cizero);
174: VecZeroEntries(auglag->Cizero);
175: }
176: if (!auglag->Ps) { /* slack vars */
177: VecDuplicate(auglag->Ci, &auglag->Ps);
178: }
179: if (!auglag->LgradS) { /* slack component of Lagrangian gradient */
180: VecDuplicate(auglag->Ci, &auglag->LgradS);
181: }
182: /* create vector for combined primal space and the associated communication objects */
183: if (!auglag->P) {
184: PetscMalloc1(2, &auglag->Parr);
185: auglag->Parr[0] = auglag->Px;
186: auglag->Parr[1] = auglag->Ps;
187: VecConcatenate(2, auglag->Parr, &auglag->P, &auglag->Pis);
188: PetscMalloc1(2, &auglag->Pscatter);
189: VecScatterCreate(auglag->P, auglag->Pis[0], auglag->Px, NULL, &auglag->Pscatter[0]);
190: VecScatterCreate(auglag->P, auglag->Pis[1], auglag->Ps, NULL, &auglag->Pscatter[1]);
191: }
192: if (tao->eq_constrained) {
193: /* create vector for combined dual space and the associated communication objects */
194: if (!auglag->Y) {
195: PetscMalloc1(2, &auglag->Yarr);
196: auglag->Yarr[0] = auglag->Ye;
197: auglag->Yarr[1] = auglag->Yi;
198: VecConcatenate(2, auglag->Yarr, &auglag->Y, &auglag->Yis);
199: PetscMalloc1(2, &auglag->Yscatter);
200: VecScatterCreate(auglag->Y, auglag->Yis[0], auglag->Ye, NULL, &auglag->Yscatter[0]);
201: VecScatterCreate(auglag->Y, auglag->Yis[1], auglag->Yi, NULL, &auglag->Yscatter[1]);
202: }
203: if (!auglag->C) VecDuplicate(auglag->Y, &auglag->C);
204: } else {
205: if (!auglag->C) auglag->C = auglag->Ci;
206: if (!auglag->Y) auglag->Y = auglag->Yi;
207: }
208: } else {
209: if (!auglag->P) auglag->P = auglag->Px;
210: if (!auglag->G) auglag->G = auglag->LgradX;
211: if (!auglag->C) auglag->C = auglag->Ce;
212: if (!auglag->Y) auglag->Y = auglag->Ye;
213: }
214: /* create tao primal solution and gradient to interface with subsolver */
215: VecDuplicate(auglag->P, &auglag->Psub);
216: VecDuplicate(auglag->P, &auglag->G);
217: /* initialize parameters */
218: if (auglag->type == TAO_ALMM_PHR) {
219: auglag->mu_fac = 10.0;
220: auglag->yi_min = 0.0;
221: auglag->ytol0 = 0.5;
222: auglag->gtol0 = tao->gatol;
223: if (tao->gatol_changed && tao->catol_changed) {
224: PetscInfo(tao, "TAOALMM with PHR: different gradient and constraint tolerances are not supported, setting catol = gatol\n");
225: tao->catol = tao->gatol;
226: }
227: }
228: /* set the Lagrangian formulation type for the subsolver */
229: switch (auglag->type) {
230: case TAO_ALMM_CLASSIC:
231: auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
232: break;
233: case TAO_ALMM_PHR:
234: auglag->sub_obj = TaoALMMComputePHRLagAndGradient_Private;
235: break;
236: default:
237: break;
238: }
239: /* set up the subsolver */
240: TaoSetSolution(auglag->subsolver, auglag->Psub);
241: TaoSetObjective(auglag->subsolver, TaoALMMSubsolverObjective_Private, (void *)auglag);
242: TaoSetObjectiveAndGradient(auglag->subsolver, NULL, TaoALMMSubsolverObjectiveAndGradient_Private, (void *)auglag);
243: if (tao->bounded) {
244: /* make sure that the subsolver is a bound-constrained method */
245: PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg);
246: PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm);
247: if (is_cg) {
248: TaoSetType(auglag->subsolver, TAOBNCG);
249: PetscInfo(tao, "TAOCG detected for bound-constrained problem, switching to TAOBNCG instead.");
250: }
251: if (is_lmvm) {
252: TaoSetType(auglag->subsolver, TAOBQNLS);
253: PetscInfo(tao, "TAOLMVM detected for bound-constrained problem, switching to TAOBQNLS instead.");
254: }
255: /* create lower and upper bound clone vectors for subsolver */
256: if (!auglag->PL) VecDuplicate(auglag->P, &auglag->PL);
257: if (!auglag->PU) VecDuplicate(auglag->P, &auglag->PU);
258: if (tao->ineq_constrained) {
259: /* create lower and upper bounds for slack, set lower to 0 */
260: VecDuplicate(auglag->Ci, &SL);
261: VecSet(SL, 0.0);
262: VecDuplicate(auglag->Ci, &SU);
263: VecSet(SU, PETSC_INFINITY);
264: /* combine opt var bounds with slack bounds */
265: TaoALMMCombinePrimal_Private(tao, tao->XL, SL, auglag->PL);
266: TaoALMMCombinePrimal_Private(tao, tao->XU, SU, auglag->PU);
267: /* destroy work vectors */
268: VecDestroy(&SL);
269: VecDestroy(&SU);
270: } else {
271: /* no inequality constraints, just copy bounds into the subsolver */
272: VecCopy(tao->XL, auglag->PL);
273: VecCopy(tao->XU, auglag->PU);
274: }
275: TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU);
276: }
277: TaoSetUp(auglag->subsolver);
279: return 0;
280: }
282: static PetscErrorCode TaoDestroy_ALMM(Tao tao)
283: {
284: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
286: TaoDestroy(&auglag->subsolver);
287: VecDestroy(&auglag->Psub);
288: VecDestroy(&auglag->G);
289: if (tao->setupcalled) {
290: VecDestroy(&auglag->Xwork);
291: if (tao->eq_constrained) {
292: VecDestroy(&auglag->Ye); /* equality multipliers */
293: VecDestroy(&auglag->Cework); /* equality work vector */
294: }
295: if (tao->ineq_constrained) {
296: VecDestroy(&auglag->Ps); /* slack vars */
297: PetscFree(auglag->Parr); /* array of primal vectors */
298: VecDestroy(&auglag->LgradS); /* slack grad */
299: VecDestroy(&auglag->Cizero); /* zero vector for pointwise max */
300: VecDestroy(&auglag->Yi); /* inequality multipliers */
301: VecDestroy(&auglag->Ciwork); /* inequality work vector */
302: VecDestroy(&auglag->P); /* combo primal solution */
303: ISDestroy(&auglag->Pis[0]); /* index set for X inside P */
304: ISDestroy(&auglag->Pis[1]); /* index set for S inside P */
305: PetscFree(auglag->Pis); /* array of P index sets */
306: VecScatterDestroy(&auglag->Pscatter[0]);
307: VecScatterDestroy(&auglag->Pscatter[1]);
308: PetscFree(auglag->Pscatter);
309: if (tao->eq_constrained) {
310: VecDestroy(&auglag->Y); /* combo multipliers */
311: PetscFree(auglag->Yarr); /* array of dual vectors */
312: VecDestroy(&auglag->C); /* combo constraints */
313: ISDestroy(&auglag->Yis[0]); /* index set for Ye inside Y */
314: ISDestroy(&auglag->Yis[1]); /* index set for Yi inside Y */
315: PetscFree(auglag->Yis);
316: VecScatterDestroy(&auglag->Yscatter[0]);
317: VecScatterDestroy(&auglag->Yscatter[1]);
318: PetscFree(auglag->Yscatter);
319: }
320: }
321: if (tao->bounded) {
322: VecDestroy(&auglag->PL); /* lower bounds for subsolver */
323: VecDestroy(&auglag->PU); /* upper bounds for subsolver */
324: }
325: }
326: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", NULL);
327: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", NULL);
328: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", NULL);
329: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", NULL);
330: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", NULL);
331: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", NULL);
332: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", NULL);
333: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", NULL);
334: PetscFree(tao->data);
335: return 0;
336: }
338: static PetscErrorCode TaoSetFromOptions_ALMM(Tao tao, PetscOptionItems *PetscOptionsObject)
339: {
340: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
341: PetscInt i;
343: PetscOptionsHeadBegin(PetscOptionsObject, "Augmented Lagrangian multiplier method solves problems with general constraints by converting them into a sequence of unconstrained problems.");
344: PetscOptionsReal("-tao_almm_mu_init", "initial penalty parameter", "", auglag->mu0, &auglag->mu0, NULL);
345: PetscOptionsReal("-tao_almm_mu_factor", "increase factor for the penalty parameter", "", auglag->mu_fac, &auglag->mu_fac, NULL);
346: PetscOptionsReal("-tao_almm_mu_power_good", "exponential for penalty parameter when multiplier update is accepted", "", auglag->mu_pow_good, &auglag->mu_pow_good, NULL);
347: PetscOptionsReal("-tao_almm_mu_power_bad", "exponential for penalty parameter when multiplier update is rejected", "", auglag->mu_pow_bad, &auglag->mu_pow_bad, NULL);
348: PetscOptionsReal("-tao_almm_mu_max", "maximum safeguard for penalty parameter updates", "", auglag->mu_max, &auglag->mu_max, NULL);
349: PetscOptionsReal("-tao_almm_ye_min", "minimum safeguard for equality multiplier updates", "", auglag->ye_min, &auglag->ye_min, NULL);
350: PetscOptionsReal("-tao_almm_ye_max", "maximum safeguard for equality multipliers updates", "", auglag->ye_max, &auglag->ye_max, NULL);
351: PetscOptionsReal("-tao_almm_yi_min", "minimum safeguard for inequality multipliers updates", "", auglag->yi_min, &auglag->yi_min, NULL);
352: PetscOptionsReal("-tao_almm_yi_max", "maximum safeguard for inequality multipliers updates", "", auglag->yi_max, &auglag->yi_max, NULL);
353: PetscOptionsEnum("-tao_almm_type", "augmented Lagrangian formulation type for the subproblem", "TaoALMMType", TaoALMMTypes, (PetscEnum)auglag->type, (PetscEnum *)&auglag->type, NULL);
354: PetscOptionsHeadEnd();
355: TaoSetOptionsPrefix(auglag->subsolver, ((PetscObject)tao)->prefix);
356: TaoAppendOptionsPrefix(auglag->subsolver, "tao_almm_subsolver_");
357: TaoSetFromOptions(auglag->subsolver);
358: for (i = 0; i < tao->numbermonitors; i++) {
359: PetscObjectReference((PetscObject)tao->monitorcontext[i]);
360: TaoSetMonitor(auglag->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);
361: if (tao->monitor[i] == TaoMonitorDefault || tao->monitor[i] == TaoDefaultCMonitor || tao->monitor[i] == TaoDefaultGMonitor || tao->monitor[i] == TaoDefaultSMonitor) auglag->info = PETSC_TRUE;
362: }
363: return 0;
364: }
366: /* -------------------------------------------------------- */
368: /*MC
369: TAOALMM - Augmented Lagrangian multiplier method for solving nonlinear optimization problems with general constraints.
371: Options Database Keys:
372: + -tao_almm_mu_init <real> - initial penalty parameter (default: 10.)
373: . -tao_almm_mu_factor <real> - increase factor for the penalty parameter (default: 100.)
374: . -tao_almm_mu_max <real> - maximum safeguard for penalty parameter updates (default: 1.e20)
375: . -tao_almm_mu_power_good <real> - exponential for penalty parameter when multiplier update is accepted (default: 0.9)
376: . -tao_almm_mu_power_bad <real> - exponential for penalty parameter when multiplier update is rejected (default: 0.1)
377: . -tao_almm_ye_min <real> - minimum safeguard for equality multiplier updates (default: -1.e20)
378: . -tao_almm_ye_max <real> - maximum safeguard for equality multiplier updates (default: 1.e20)
379: . -tao_almm_yi_min <real> - minimum safeguard for inequality multiplier updates (default: -1.e20)
380: . -tao_almm_yi_max <real> - maximum safeguard for inequality multiplier updates (default: 1.e20)
381: - -tao_almm_type <phr,classic> - change formulation of the augmented Lagrangian merit function for the subproblem (default: phr)
383: Level: beginner
385: Notes:
386: This method converts a constrained problem into a sequence of unconstrained problems via the augmented
387: Lagrangian merit function. Bound constraints are pushed down to the subproblem without any modifications.
389: Two formulations are offered for the subproblem: canonical Hestenes-Powell augmented Lagrangian with slack
390: variables for inequality constraints, and a slack-less Powell-Hestenes-Rockafellar (PHR) formulation utilizing a
391: pointwise max() penalty on inequality constraints. The canonical augmented Lagrangian formulation may converge
392: faster for smaller problems but is highly susceptible to poor step lengths in the subproblem due to the positivity
393: constraint on slack variables. PHR avoids this issue by eliminating the slack variables entirely, and is highly
394: desirable for problems with a large number of inequality constraints.
396: The subproblem is solved using a nested first-order TAO solver (default: `TAOBQNLS`). The user can retrieve a
397: pointer to the subsolver via `TaoALMMGetSubsolver()` or pass command line arguments to it using the
398: "-tao_almm_subsolver_" prefix. Currently, `TAOALMM` does not support second-order methods for the subproblem.
400: .vb
401: while unconverged
402: solve argmin_x L(x) s.t. l <= x <= u
403: if ||c|| <= y_tol
404: if ||c|| <= c_tol && ||Lgrad|| <= g_tol:
405: problem converged, return solution
406: else
407: constraints sufficiently improved
408: update multipliers and tighten tolerances
409: endif
410: else
411: constraints did not improve
412: update penalty and loosen tolerances
413: endif
414: endwhile
415: .ve
417: .seealso: `TAOALMM`, `Tao`, `TaoALMMGetType()`, `TaoALMMSetType()`, `TaoALMMSetSubsolver()`, `TaoALMMGetSubsolver()`,
418: `TaoALMMGetMultipliers()`, `TaoALMMSetMultipliers()`, `TaoALMMGetPrimalIS()`, `TaoALMMGetDualIS()`
419: M*/
420: PETSC_EXTERN PetscErrorCode TaoCreate_ALMM(Tao tao)
421: {
422: TAO_ALMM *auglag;
424: PetscNew(&auglag);
426: tao->ops->destroy = TaoDestroy_ALMM;
427: tao->ops->setup = TaoSetUp_ALMM;
428: tao->ops->setfromoptions = TaoSetFromOptions_ALMM;
429: tao->ops->view = TaoView_ALMM;
430: tao->ops->solve = TaoSolve_ALMM;
432: tao->gatol = 1.e-5;
433: tao->grtol = 0.0;
434: tao->gttol = 0.0;
435: tao->catol = 1.e-5;
436: tao->crtol = 0.0;
438: tao->data = (void *)auglag;
439: auglag->parent = tao;
440: auglag->mu0 = 10.0;
441: auglag->mu = auglag->mu0;
442: auglag->mu_fac = 10.0;
443: auglag->mu_max = PETSC_INFINITY;
444: auglag->mu_pow_good = 0.9;
445: auglag->mu_pow_bad = 0.1;
446: auglag->ye_min = PETSC_NINFINITY;
447: auglag->ye_max = PETSC_INFINITY;
448: auglag->yi_min = PETSC_NINFINITY;
449: auglag->yi_max = PETSC_INFINITY;
450: auglag->ytol0 = 0.1 / PetscPowReal(auglag->mu0, auglag->mu_pow_bad);
451: auglag->ytol = auglag->ytol0;
452: auglag->gtol0 = 1.0 / auglag->mu0;
453: auglag->gtol = auglag->gtol0;
455: auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private;
456: auglag->type = TAO_ALMM_PHR;
457: auglag->info = PETSC_FALSE;
459: TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver);
460: TaoSetType(auglag->subsolver, TAOBQNLS);
461: TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0);
462: TaoSetMaximumIterations(auglag->subsolver, 1000);
463: TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000);
464: TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY);
465: PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1);
467: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private);
468: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private);
469: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private);
470: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private);
471: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private);
472: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private);
473: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private);
474: PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private);
475: return 0;
476: }
478: static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P)
479: {
480: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
482: if (tao->ineq_constrained) {
483: VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);
484: VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);
485: VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);
486: VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);
487: } else {
488: VecCopy(X, P);
489: }
490: return 0;
491: }
493: static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y)
494: {
495: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
497: if (tao->eq_constrained) {
498: if (tao->ineq_constrained) {
499: VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);
500: VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);
501: VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);
502: VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);
503: } else {
504: VecCopy(EQ, Y);
505: }
506: } else {
507: VecCopy(IN, Y);
508: }
509: return 0;
510: }
512: static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S)
513: {
514: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
516: if (tao->ineq_constrained) {
517: VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);
518: VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);
519: VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);
520: VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);
521: } else {
522: VecCopy(P, X);
523: }
524: return 0;
525: }
527: /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */
528: static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao)
529: {
530: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
532: /* if bounded, project the gradient */
533: if (tao->bounded) VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX);
534: if (auglag->type == TAO_ALMM_PHR) {
535: VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm);
536: auglag->cenorm = 0.0;
537: if (tao->eq_constrained) VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm);
538: auglag->cinorm = 0.0;
539: if (tao->ineq_constrained) {
540: VecCopy(auglag->Yi, auglag->Ciwork);
541: VecScale(auglag->Ciwork, -1.0 / auglag->mu);
542: VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork);
543: VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm);
544: }
545: auglag->cnorm_old = auglag->cnorm;
546: auglag->cnorm = PetscMax(auglag->cenorm, auglag->cinorm);
547: auglag->ytol = auglag->ytol0 * auglag->cnorm_old;
548: } else {
549: VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm);
550: VecNorm(auglag->C, NORM_2, &auglag->cnorm);
551: }
552: return 0;
553: }
555: static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao)
556: {
557: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
559: /* split solution into primal and slack components */
560: TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps);
562: /* compute f, df/dx and the constraints */
563: TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX);
564: if (tao->eq_constrained) {
565: TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce);
566: TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae);
567: }
568: if (tao->ineq_constrained) {
569: TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci);
570: TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai);
571: switch (auglag->type) {
572: case TAO_ALMM_CLASSIC:
573: /* classic formulation converts inequality to equality constraints via slack variables */
574: VecAXPY(auglag->Ci, -1.0, auglag->Ps);
575: break;
576: case TAO_ALMM_PHR:
577: /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */
578: VecScale(auglag->Ci, -1.0);
579: MatScale(auglag->Ai, -1.0);
580: break;
581: default:
582: break;
583: }
584: }
585: /* combine constraints into one vector */
586: TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C);
587: return 0;
588: }
590: /*
591: Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmin(0, Ci + Yi/mu)^T pmin(0, Ci + Yi/mu)]
593: dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmin(0, Ci + Yi/mu)^T Ai]
595: dLphr/dS = 0
596: */
597: static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao)
598: {
599: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
600: PetscReal eq_norm = 0.0, ineq_norm = 0.0;
602: TaoALMMEvaluateIterate_Private(tao);
603: if (tao->eq_constrained) {
604: /* Ce_work = mu*(Ce + Ye/mu) */
605: VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce);
606: VecDot(auglag->Cework, auglag->Cework, &eq_norm); /* contribution to scalar Lagrangian */
607: VecScale(auglag->Cework, auglag->mu);
608: /* dL/dX += mu*(Ce + Ye/mu)^T Ae */
609: MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX);
610: }
611: if (tao->ineq_constrained) {
612: /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */
613: VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci);
614: VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork);
615: VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm); /* contribution to scalar Lagrangian */
616: /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */
617: VecScale(auglag->Ciwork, auglag->mu);
618: MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX);
619: /* dL/dS = 0 because there are no slacks in PHR */
620: VecZeroEntries(auglag->LgradS);
621: }
622: /* combine gradient together */
623: TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);
624: /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */
625: auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm);
626: return 0;
627: }
629: /*
630: Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)]
632: dLc/dX = dF/dX + Ye^TAe + Yi^TAi + 0.5*mu*[Ce^TAe + (Ci - S)^TAi]
634: dLc/dS = -[Yi + mu*(Ci - S)]
635: */
636: static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao)
637: {
638: TAO_ALMM *auglag = (TAO_ALMM *)tao->data;
639: PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0;
641: TaoALMMEvaluateIterate_Private(tao);
642: if (tao->eq_constrained) {
643: /* compute scalar contributions */
644: VecDot(auglag->Ye, auglag->Ce, &yeTce);
645: VecDot(auglag->Ce, auglag->Ce, &ceTce);
646: /* dL/dX += ye^T Ae */
647: MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX);
648: /* dL/dX += 0.5 * mu * ce^T Ae */
649: MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork);
650: VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork);
651: }
652: if (tao->ineq_constrained) {
653: /* compute scalar contributions */
654: VecDot(auglag->Yi, auglag->Ci, &yiTcims);
655: VecDot(auglag->Ci, auglag->Ci, &cimsTcims);
656: /* dL/dX += yi^T Ai */
657: MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX);
658: /* dL/dX += 0.5 * mu * (ci - s)^T Ai */
659: MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork);
660: VecAXPY(auglag->LgradX, 0.5 * auglag->mu, auglag->Xwork);
661: /* dL/dS = -[yi + mu*(ci - s)] */
662: VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi);
663: VecScale(auglag->LgradS, -1.0);
664: }
665: /* combine gradient together */
666: TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);
667: /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */
668: auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims);
669: return 0;
670: }
672: PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, void *ctx)
673: {
674: TAO_ALMM *auglag = (TAO_ALMM *)ctx;
676: VecCopy(P, auglag->P);
677: (*auglag->sub_obj)(auglag->parent);
678: *Lval = auglag->Lval;
679: return 0;
680: }
682: PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx)
683: {
684: TAO_ALMM *auglag = (TAO_ALMM *)ctx;
686: VecCopy(P, auglag->P);
687: (*auglag->sub_obj)(auglag->parent);
688: VecCopy(auglag->G, G);
689: *Lval = auglag->Lval;
690: return 0;
691: }