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 /* -------------------------------------------------------- */ 415 416 /*MC 417 TAOALMM - Augmented Lagrangian multiplier method for solving nonlinear optimization problems with general constraints. 418 419 Options Database Keys: 420 + -tao_almm_mu_init <real> - initial penalty parameter (default: 10.) 421 . -tao_almm_mu_factor <real> - increase factor for the penalty parameter (default: 100.) 422 . -tao_almm_mu_max <real> - maximum safeguard for penalty parameter updates (default: 1.e20) 423 . -tao_almm_mu_power_good <real> - exponential for penalty parameter when multiplier update is accepted (default: 0.9) 424 . -tao_almm_mu_power_bad <real> - exponential for penalty parameter when multiplier update is rejected (default: 0.1) 425 . -tao_almm_ye_min <real> - minimum safeguard for equality multiplier updates (default: -1.e20) 426 . -tao_almm_ye_max <real> - maximum safeguard for equality multiplier updates (default: 1.e20) 427 . -tao_almm_yi_min <real> - minimum safeguard for inequality multiplier updates (default: -1.e20) 428 . -tao_almm_yi_max <real> - maximum safeguard for inequality multiplier updates (default: 1.e20) 429 - -tao_almm_type <phr,classic> - change formulation of the augmented Lagrangian merit function for the subproblem (default: phr) 430 431 Level: beginner 432 433 Notes: 434 This method converts a constrained problem into a sequence of unconstrained problems via the augmented 435 Lagrangian merit function. Bound constraints are pushed down to the subproblem without any modifications. 436 437 Two formulations are offered for the subproblem: canonical Hestenes-Powell augmented Lagrangian with slack 438 variables for inequality constraints, and a slack-less Powell-Hestenes-Rockafellar (PHR) formulation utilizing a 439 pointwise max() penalty on inequality constraints. The canonical augmented Lagrangian formulation may converge 440 faster for smaller problems but is highly susceptible to poor step lengths in the subproblem due to the positivity 441 constraint on slack variables. PHR avoids this issue by eliminating the slack variables entirely, and is highly 442 desirable for problems with a large number of inequality constraints. 443 444 The subproblem is solved using a nested first-order TAO solver (default: `TAOBQNLS`). The user can retrieve a 445 pointer to the subsolver via `TaoALMMGetSubsolver()` or pass command line arguments to it using the 446 "-tao_almm_subsolver_" prefix. Currently, `TAOALMM` does not support second-order methods for the subproblem. 447 448 .vb 449 while unconverged 450 solve argmin_x L(x) s.t. l <= x <= u 451 if ||c|| <= y_tol 452 if ||c|| <= c_tol && ||Lgrad|| <= g_tol: 453 problem converged, return solution 454 else 455 constraints sufficiently improved 456 update multipliers and tighten tolerances 457 endif 458 else 459 constraints did not improve 460 update penalty and loosen tolerances 461 endif 462 endwhile 463 .ve 464 465 .seealso: `TAOALMM`, `Tao`, `TaoALMMGetType()`, `TaoALMMSetType()`, `TaoALMMSetSubsolver()`, `TaoALMMGetSubsolver()`, 466 `TaoALMMGetMultipliers()`, `TaoALMMSetMultipliers()`, `TaoALMMGetPrimalIS()`, `TaoALMMGetDualIS()` 467 M*/ 468 PETSC_EXTERN PetscErrorCode TaoCreate_ALMM(Tao tao) 469 { 470 TAO_ALMM *auglag; 471 472 PetscFunctionBegin; 473 PetscCall(PetscNew(&auglag)); 474 475 tao->ops->destroy = TaoDestroy_ALMM; 476 tao->ops->setup = TaoSetUp_ALMM; 477 tao->ops->setfromoptions = TaoSetFromOptions_ALMM; 478 tao->ops->view = TaoView_ALMM; 479 tao->ops->solve = TaoSolve_ALMM; 480 481 PetscCall(TaoParametersInitialize(tao)); 482 PetscObjectParameterSetDefault(tao, gatol, 1.e-5); 483 PetscObjectParameterSetDefault(tao, grtol, 0.0); 484 PetscObjectParameterSetDefault(tao, gttol, 0.0); 485 PetscObjectParameterSetDefault(tao, catol, 1.e-5); 486 PetscObjectParameterSetDefault(tao, crtol, 0.0); 487 488 tao->data = (void *)auglag; 489 auglag->parent = tao; 490 auglag->mu0 = 10.0; 491 auglag->mu = auglag->mu0; 492 auglag->mu_fac = 10.0; 493 auglag->mu_max = PETSC_INFINITY; 494 auglag->mu_pow_good = 0.9; 495 auglag->mu_pow_bad = 0.1; 496 auglag->ye_min = PETSC_NINFINITY; 497 auglag->ye_max = PETSC_INFINITY; 498 auglag->yi_min = PETSC_NINFINITY; 499 auglag->yi_max = PETSC_INFINITY; 500 auglag->ytol0 = 1.0 / PetscPowReal(auglag->mu0, auglag->mu_pow_bad); 501 auglag->ytol = auglag->ytol0; 502 auglag->gtol0 = 1.0 / auglag->mu0; 503 auglag->gtol = auglag->gtol0; 504 505 auglag->sub_obj = TaoALMMComputeAugLagAndGradient_Private; 506 auglag->type = TAO_ALMM_PHR; 507 auglag->info = PETSC_FALSE; 508 509 PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &auglag->subsolver)); 510 PetscCall(TaoSetType(auglag->subsolver, TAOBQNLS)); 511 PetscCall(TaoSetTolerances(auglag->subsolver, auglag->gtol, 0.0, 0.0)); 512 PetscCall(TaoSetMaximumIterations(auglag->subsolver, 1000)); 513 PetscCall(TaoSetMaximumFunctionEvaluations(auglag->subsolver, 10000)); 514 PetscCall(TaoSetFunctionLowerBound(auglag->subsolver, PETSC_NINFINITY)); 515 PetscCall(PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver, (PetscObject)tao, 1)); 516 517 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private)); 518 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private)); 519 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private)); 520 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private)); 521 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private)); 522 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private)); 523 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private)); 524 PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private)); 525 PetscFunctionReturn(PETSC_SUCCESS); 526 } 527 528 static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P) 529 { 530 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 531 532 PetscFunctionBegin; 533 if (tao->ineq_constrained) { 534 PetscCall(VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE)); 535 PetscCall(VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE)); 536 PetscCall(VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE)); 537 PetscCall(VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE)); 538 } else { 539 PetscCall(VecCopy(X, P)); 540 } 541 PetscFunctionReturn(PETSC_SUCCESS); 542 } 543 544 static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y) 545 { 546 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 547 548 PetscFunctionBegin; 549 if (tao->eq_constrained) { 550 if (tao->ineq_constrained) { 551 PetscCall(VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE)); 552 PetscCall(VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE)); 553 PetscCall(VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE)); 554 PetscCall(VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE)); 555 } else { 556 PetscCall(VecCopy(EQ, Y)); 557 } 558 } else { 559 PetscCall(VecCopy(IN, Y)); 560 } 561 PetscFunctionReturn(PETSC_SUCCESS); 562 } 563 564 static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S) 565 { 566 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 567 568 PetscFunctionBegin; 569 if (tao->ineq_constrained) { 570 PetscCall(VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD)); 571 PetscCall(VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD)); 572 PetscCall(VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD)); 573 PetscCall(VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD)); 574 } else { 575 PetscCall(VecCopy(P, X)); 576 } 577 PetscFunctionReturn(PETSC_SUCCESS); 578 } 579 580 /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */ 581 static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao) 582 { 583 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 584 585 PetscFunctionBegin; 586 /* if bounded, project the gradient */ 587 if (tao->bounded) PetscCall(VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX)); 588 if (auglag->type == TAO_ALMM_PHR) { 589 PetscCall(VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm)); 590 auglag->cenorm = 0.0; 591 if (tao->eq_constrained) PetscCall(VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm)); 592 auglag->cinorm = 0.0; 593 if (tao->ineq_constrained) { 594 PetscCall(VecCopy(auglag->Yi, auglag->Ciwork)); 595 PetscCall(VecScale(auglag->Ciwork, -1.0 / auglag->mu)); 596 PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork)); 597 PetscCall(VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm)); 598 } 599 auglag->cnorm_old = auglag->cnorm; 600 auglag->cnorm = PetscMax(auglag->cenorm, auglag->cinorm); 601 auglag->ytol = auglag->ytol0 * auglag->cnorm_old; 602 } else { 603 PetscCall(VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm)); 604 PetscCall(VecNorm(auglag->C, NORM_2, &auglag->cnorm)); 605 } 606 PetscFunctionReturn(PETSC_SUCCESS); 607 } 608 609 static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao) 610 { 611 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 612 613 PetscFunctionBegin; 614 /* split solution into primal and slack components */ 615 PetscCall(TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps)); 616 617 /* compute f, df/dx and the constraints */ 618 PetscCall(TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX)); 619 if (tao->eq_constrained) { 620 PetscCall(TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce)); 621 PetscCall(TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae)); 622 } 623 if (tao->ineq_constrained) { 624 PetscCall(TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci)); 625 PetscCall(TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai)); 626 switch (auglag->type) { 627 case TAO_ALMM_CLASSIC: 628 /* classic formulation converts inequality to equality constraints via slack variables */ 629 PetscCall(VecAXPY(auglag->Ci, -1.0, auglag->Ps)); 630 break; 631 case TAO_ALMM_PHR: 632 /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */ 633 PetscCall(VecScale(auglag->Ci, -1.0)); 634 PetscCall(MatScale(auglag->Ai, -1.0)); 635 break; 636 default: 637 break; 638 } 639 } 640 /* combine constraints into one vector */ 641 PetscCall(TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C)); 642 PetscFunctionReturn(PETSC_SUCCESS); 643 } 644 645 /* 646 Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] 647 648 dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmax(0, Ci + Yi/mu)^T Ai] 649 650 dLphr/dS = 0 651 */ 652 static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao) 653 { 654 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 655 PetscReal eq_norm = 0.0, ineq_norm = 0.0; 656 657 PetscFunctionBegin; 658 PetscCall(TaoALMMEvaluateIterate_Private(tao)); 659 if (tao->eq_constrained) { 660 /* Ce_work = mu*(Ce + Ye/mu) */ 661 PetscCall(VecWAXPY(auglag->Cework, 1.0 / auglag->mu, auglag->Ye, auglag->Ce)); 662 PetscCall(VecDot(auglag->Cework, auglag->Cework, &eq_norm)); /* contribution to scalar Lagrangian */ 663 PetscCall(VecScale(auglag->Cework, auglag->mu)); 664 /* dL/dX += mu*(Ce + Ye/mu)^T Ae */ 665 PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX)); 666 } 667 if (tao->ineq_constrained) { 668 /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */ 669 PetscCall(VecWAXPY(auglag->Ciwork, 1.0 / auglag->mu, auglag->Yi, auglag->Ci)); 670 PetscCall(VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork)); 671 PetscCall(VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm)); /* contribution to scalar Lagrangian */ 672 /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */ 673 PetscCall(VecScale(auglag->Ciwork, auglag->mu)); 674 PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX)); 675 /* dL/dS = 0 because there are no slacks in PHR */ 676 PetscCall(VecZeroEntries(auglag->LgradS)); 677 } 678 /* combine gradient together */ 679 PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G)); 680 /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */ 681 auglag->Lval = auglag->fval + 0.5 * auglag->mu * (eq_norm + ineq_norm); 682 PetscFunctionReturn(PETSC_SUCCESS); 683 } 684 685 /* 686 Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)] 687 688 dLc/dX = dF/dX + Ye^TAe + Yi^TAi + mu*[Ce^TAe + (Ci - S)^TAi] 689 690 dLc/dS = -[Yi + mu*(Ci - S)] 691 */ 692 static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao) 693 { 694 TAO_ALMM *auglag = (TAO_ALMM *)tao->data; 695 PetscReal yeTce = 0.0, yiTcims = 0.0, ceTce = 0.0, cimsTcims = 0.0; 696 697 PetscFunctionBegin; 698 PetscCall(TaoALMMEvaluateIterate_Private(tao)); 699 if (tao->eq_constrained) { 700 /* compute scalar contributions */ 701 PetscCall(VecDot(auglag->Ye, auglag->Ce, &yeTce)); 702 PetscCall(VecDot(auglag->Ce, auglag->Ce, &ceTce)); 703 /* dL/dX += ye^T Ae */ 704 PetscCall(MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX)); 705 /* dL/dX += mu * ce^T Ae */ 706 PetscCall(MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork)); 707 PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork)); 708 } 709 if (tao->ineq_constrained) { 710 /* compute scalar contributions */ 711 PetscCall(VecDot(auglag->Yi, auglag->Ci, &yiTcims)); 712 PetscCall(VecDot(auglag->Ci, auglag->Ci, &cimsTcims)); 713 /* dL/dX += yi^T Ai */ 714 PetscCall(MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX)); 715 /* dL/dX += mu * (ci - s)^T Ai */ 716 PetscCall(MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork)); 717 PetscCall(VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork)); 718 /* dL/dS = -[yi + mu*(ci - s)] */ 719 PetscCall(VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi)); 720 PetscCall(VecScale(auglag->LgradS, -1.0)); 721 } 722 /* combine gradient together */ 723 PetscCall(TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G)); 724 /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */ 725 auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5 * auglag->mu * (ceTce + cimsTcims); 726 PetscFunctionReturn(PETSC_SUCCESS); 727 } 728 729 PetscErrorCode TaoALMMSubsolverObjective_Private(Tao tao, Vec P, PetscReal *Lval, void *ctx) 730 { 731 TAO_ALMM *auglag = (TAO_ALMM *)ctx; 732 733 PetscFunctionBegin; 734 PetscCall(VecCopy(P, auglag->P)); 735 PetscCall((*auglag->sub_obj)(auglag->parent)); 736 *Lval = auglag->Lval; 737 PetscFunctionReturn(PETSC_SUCCESS); 738 } 739 740 PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx) 741 { 742 TAO_ALMM *auglag = (TAO_ALMM *)ctx; 743 744 PetscFunctionBegin; 745 PetscCall(VecCopy(P, auglag->P)); 746 PetscCall((*auglag->sub_obj)(auglag->parent)); 747 PetscCall(VecCopy(auglag->G, G)); 748 *Lval = auglag->Lval; 749 PetscFunctionReturn(PETSC_SUCCESS); 750 } 751