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