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