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 = TaoSetObjectiveAndGradientRoutine(auglag->subsolver, TaoALMMSubsolverObjectiveAndGradient_Private, (void*)auglag);CHKERRQ(ierr); 265 if (tao->bounded) { 266 /* make sure that the subsolver is a bound-constrained method */ 267 ierr = PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOCG, &is_cg);CHKERRQ(ierr); 268 ierr = PetscObjectTypeCompare((PetscObject)auglag->subsolver, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 269 if (is_cg) { 270 ierr = TaoSetType(auglag->subsolver, TAOBNCG);CHKERRQ(ierr); 271 ierr = PetscInfo(tao,"TAOCG detected for bound-constrained problem, switching to TAOBNCG instead.");CHKERRQ(ierr); 272 } 273 if (is_lmvm) { 274 ierr = TaoSetType(auglag->subsolver, TAOBQNLS);CHKERRQ(ierr); 275 ierr = PetscInfo(tao,"TAOLMVM detected for bound-constrained problem, switching to TAOBQNLS instead.");CHKERRQ(ierr); 276 } 277 /* create lower and upper bound clone vectors for subsolver */ 278 if (!auglag->PL) { 279 ierr = VecDuplicate(auglag->P, &auglag->PL);CHKERRQ(ierr); 280 } 281 if (!auglag->PU) { 282 ierr = VecDuplicate(auglag->P, &auglag->PU);CHKERRQ(ierr); 283 } 284 if (tao->ineq_constrained) { 285 /* create lower and upper bounds for slack, set lower to 0 */ 286 ierr = VecDuplicate(auglag->Ci, &SL);CHKERRQ(ierr); 287 ierr = VecSet(SL, 0.0);CHKERRQ(ierr); 288 ierr = VecDuplicate(auglag->Ci, &SU);CHKERRQ(ierr); 289 ierr = VecSet(SU, PETSC_INFINITY);CHKERRQ(ierr); 290 /* combine opt var bounds with slack bounds */ 291 ierr = TaoALMMCombinePrimal_Private(tao, tao->XL, SL, auglag->PL);CHKERRQ(ierr); 292 ierr = TaoALMMCombinePrimal_Private(tao, tao->XU, SU, auglag->PU);CHKERRQ(ierr); 293 /* destroy work vectors */ 294 ierr = VecDestroy(&SL);CHKERRQ(ierr); 295 ierr = VecDestroy(&SU);CHKERRQ(ierr); 296 } else { 297 /* no inequality constraints, just copy bounds into the subsolver */ 298 ierr = VecCopy(tao->XL, auglag->PL);CHKERRQ(ierr); 299 ierr = VecCopy(tao->XU, auglag->PU);CHKERRQ(ierr); 300 } 301 ierr = TaoSetVariableBounds(auglag->subsolver, auglag->PL, auglag->PU);CHKERRQ(ierr); 302 } 303 ierr = TaoSetUp(auglag->subsolver);CHKERRQ(ierr); 304 auglag->G = auglag->subsolver->gradient; 305 306 PetscFunctionReturn(0); 307 } 308 309 static PetscErrorCode TaoDestroy_ALMM(Tao tao) 310 { 311 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 312 PetscErrorCode ierr; 313 314 PetscFunctionBegin; 315 ierr = TaoDestroy(&auglag->subsolver);CHKERRQ(ierr); 316 if (tao->setupcalled) { 317 ierr = VecDestroy(&auglag->Xwork);CHKERRQ(ierr); /* opt work */ 318 if (tao->eq_constrained) { 319 ierr = VecDestroy(&auglag->Ye);CHKERRQ(ierr); /* equality multipliers */ 320 ierr = VecDestroy(&auglag->Cework);CHKERRQ(ierr); /* equality work vector */ 321 } 322 if (tao->ineq_constrained) { 323 ierr = VecDestroy(&auglag->Ps);CHKERRQ(ierr); /* slack vars */ 324 auglag->Parr[0] = NULL; /* clear pointer to tao->solution, will be destroyed by TaoDestroy() shell */ 325 ierr = PetscFree(auglag->Parr);CHKERRQ(ierr); /* array of primal vectors */ 326 ierr = VecDestroy(&auglag->LgradS);CHKERRQ(ierr); /* slack grad */ 327 ierr = VecDestroy(&auglag->Cizero);CHKERRQ(ierr); /* zero vector for pointwise max */ 328 ierr = VecDestroy(&auglag->Yi);CHKERRQ(ierr); /* inequality multipliers */ 329 ierr = VecDestroy(&auglag->Ciwork);CHKERRQ(ierr); /* inequality work vector */ 330 ierr = VecDestroy(&auglag->P);CHKERRQ(ierr); /* combo primal */ 331 ierr = ISDestroy(&auglag->Pis[0]);CHKERRQ(ierr); /* index set for X inside P */ 332 ierr = ISDestroy(&auglag->Pis[1]);CHKERRQ(ierr); /* index set for S inside P */ 333 ierr = PetscFree(auglag->Pis);CHKERRQ(ierr); /* array of P index sets */ 334 ierr = VecScatterDestroy(&auglag->Pscatter[0]);CHKERRQ(ierr); 335 ierr = VecScatterDestroy(&auglag->Pscatter[1]);CHKERRQ(ierr); 336 ierr = PetscFree(auglag->Pscatter);CHKERRQ(ierr); 337 if (tao->eq_constrained) { 338 ierr = VecDestroy(&auglag->Y);CHKERRQ(ierr); /* combo multipliers */ 339 ierr = PetscFree(auglag->Yarr);CHKERRQ(ierr); /* array of dual vectors */ 340 ierr = VecDestroy(&auglag->C);CHKERRQ(ierr); /* combo constraints */ 341 ierr = ISDestroy(&auglag->Yis[0]);CHKERRQ(ierr); /* index set for Ye inside Y */ 342 ierr = ISDestroy(&auglag->Yis[1]);CHKERRQ(ierr); /* index set for Yi inside Y */ 343 ierr = PetscFree(auglag->Yis);CHKERRQ(ierr); 344 ierr = VecScatterDestroy(&auglag->Yscatter[0]);CHKERRQ(ierr); 345 ierr = VecScatterDestroy(&auglag->Yscatter[1]);CHKERRQ(ierr); 346 ierr = PetscFree(auglag->Yscatter);CHKERRQ(ierr); 347 } 348 } 349 if (tao->bounded) { 350 ierr = VecDestroy(&auglag->PL);CHKERRQ(ierr); /* lower bounds for subsolver */ 351 ierr = VecDestroy(&auglag->PU);CHKERRQ(ierr); /* upper bounds for subsolver */ 352 } 353 } 354 ierr = PetscFree(tao->data);CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357 358 static PetscErrorCode TaoSetFromOptions_ALMM(PetscOptionItems *PetscOptionsObject,Tao tao) 359 { 360 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 361 PetscInt i; 362 PetscBool compatible; 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 = TaoSetFromOptions(auglag->subsolver);CHKERRQ(ierr); 379 ierr = PetscObjectTypeCompareAny((PetscObject)(auglag->subsolver), &compatible, TAOCG, TAOLMVM, TAOBNCG, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 380 if (!compatible) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Subsolver must be a first-order method (TAOCG, TABNCG, TAOLMVM, TAOBQN-family"); 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 = TaoSetOptionsPrefix(auglag->subsolver,"tao_almm_subsolver_");CHKERRQ(ierr); 494 ierr = PetscObjectIncrementTabLevel((PetscObject)auglag->subsolver,(PetscObject)tao,1);CHKERRQ(ierr); 495 496 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetType_C", TaoALMMGetType_Private);CHKERRQ(ierr); 497 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetType_C", TaoALMMSetType_Private);CHKERRQ(ierr); 498 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetSubsolver_C", TaoALMMGetSubsolver_Private);CHKERRQ(ierr); 499 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetSubsolver_C", TaoALMMSetSubsolver_Private);CHKERRQ(ierr); 500 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetMultipliers_C", TaoALMMGetMultipliers_Private);CHKERRQ(ierr); 501 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMSetMultipliers_C", TaoALMMSetMultipliers_Private);CHKERRQ(ierr); 502 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetPrimalIS_C", TaoALMMGetPrimalIS_Private);CHKERRQ(ierr); 503 ierr = PetscObjectComposeFunction((PetscObject)tao, "TaoALMMGetDualIS_C", TaoALMMGetDualIS_Private);CHKERRQ(ierr); 504 PetscFunctionReturn(0); 505 } 506 507 static PetscErrorCode TaoALMMCombinePrimal_Private(Tao tao, Vec X, Vec S, Vec P) 508 { 509 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 510 PetscErrorCode ierr; 511 512 PetscFunctionBegin; 513 if (tao->ineq_constrained) { 514 ierr = VecScatterBegin(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 515 ierr = VecScatterEnd(auglag->Pscatter[0], X, P, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 516 ierr = VecScatterBegin(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 517 ierr = VecScatterEnd(auglag->Pscatter[1], S, P, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 518 } else { 519 ierr = VecCopy(X, P);CHKERRQ(ierr); 520 } 521 PetscFunctionReturn(0); 522 } 523 524 static PetscErrorCode TaoALMMCombineDual_Private(Tao tao, Vec EQ, Vec IN, Vec Y) 525 { 526 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 527 PetscErrorCode ierr; 528 529 PetscFunctionBegin; 530 if (tao->eq_constrained) { 531 if (tao->ineq_constrained) { 532 ierr = VecScatterBegin(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 533 ierr = VecScatterEnd(auglag->Yscatter[0], EQ, Y, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 534 ierr = VecScatterBegin(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 535 ierr = VecScatterEnd(auglag->Yscatter[1], IN, Y, INSERT_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 536 } else { 537 ierr = VecCopy(EQ, Y);CHKERRQ(ierr); 538 } 539 } else { 540 ierr = VecCopy(IN, Y);CHKERRQ(ierr); 541 } 542 PetscFunctionReturn(0); 543 } 544 545 static PetscErrorCode TaoALMMSplitPrimal_Private(Tao tao, Vec P, Vec X, Vec S) 546 { 547 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 548 PetscErrorCode ierr; 549 550 PetscFunctionBegin; 551 if (tao->ineq_constrained) { 552 ierr = VecScatterBegin(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 553 ierr = VecScatterEnd(auglag->Pscatter[0], P, X, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 554 ierr = VecScatterBegin(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 555 ierr = VecScatterEnd(auglag->Pscatter[1], P, S, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 556 } else { 557 ierr = VecCopy(P, X);CHKERRQ(ierr); 558 } 559 PetscFunctionReturn(0); 560 } 561 562 /* this assumes that the latest constraints are stored in Ce and Ci, and also combined in C */ 563 static PetscErrorCode TaoALMMComputeOptimalityNorms_Private(Tao tao) 564 { 565 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 566 PetscErrorCode ierr; 567 568 PetscFunctionBegin; 569 /* if bounded, project the gradient */ 570 if (tao->bounded) { 571 ierr = VecBoundGradientProjection(auglag->LgradX, auglag->Px, tao->XL, tao->XU, auglag->LgradX);CHKERRQ(ierr); 572 } 573 if (auglag->type == TAO_ALMM_PHR) { 574 ierr = VecNorm(auglag->LgradX, NORM_INFINITY, &auglag->gnorm);CHKERRQ(ierr); 575 auglag->cenorm = 0.0; 576 if (tao->eq_constrained) { 577 ierr = VecNorm(auglag->Ce, NORM_INFINITY, &auglag->cenorm);CHKERRQ(ierr); 578 } 579 auglag->cinorm = 0.0; 580 if (tao->ineq_constrained) { 581 ierr = VecCopy(auglag->Yi, auglag->Ciwork);CHKERRQ(ierr); 582 ierr = VecScale(auglag->Ciwork, -1.0/auglag->mu);CHKERRQ(ierr); 583 ierr = VecPointwiseMax(auglag->Ciwork, auglag->Ci, auglag->Ciwork);CHKERRQ(ierr); 584 ierr = VecNorm(auglag->Ciwork, NORM_INFINITY, &auglag->cinorm);CHKERRQ(ierr); 585 } 586 auglag->cnorm_old = auglag->cnorm; 587 auglag->cnorm = PetscMax(auglag->cenorm, auglag->cinorm); 588 auglag->ytol = auglag->ytol0 * auglag->cnorm_old; 589 } else { 590 ierr = VecNorm(auglag->LgradX, NORM_2, &auglag->gnorm);CHKERRQ(ierr); 591 ierr = VecNorm(auglag->C, NORM_2, &auglag->cnorm);CHKERRQ(ierr); 592 } 593 PetscFunctionReturn(0); 594 } 595 596 static PetscErrorCode TaoALMMEvaluateIterate_Private(Tao tao, Vec P) 597 { 598 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 599 PetscErrorCode ierr; 600 601 PetscFunctionBegin; 602 /* split solution into primal and slack components */ 603 ierr = TaoALMMSplitPrimal_Private(tao, auglag->P, auglag->Px, auglag->Ps);CHKERRQ(ierr); 604 605 /* compute f, df/dx and the constraints */ 606 ierr = TaoComputeObjectiveAndGradient(tao, auglag->Px, &auglag->fval, auglag->LgradX);CHKERRQ(ierr); 607 if (tao->eq_constrained) { 608 ierr = TaoComputeEqualityConstraints(tao, auglag->Px, auglag->Ce);CHKERRQ(ierr); 609 ierr = TaoComputeJacobianEquality(tao, auglag->Px, auglag->Ae, auglag->Ae);CHKERRQ(ierr); 610 } 611 if (tao->ineq_constrained) { 612 ierr = TaoComputeInequalityConstraints(tao, auglag->Px, auglag->Ci);CHKERRQ(ierr); 613 ierr = TaoComputeJacobianInequality(tao, auglag->Px, auglag->Ai, auglag->Ai);CHKERRQ(ierr); 614 switch (auglag->type) { 615 case (TAO_ALMM_CLASSIC): 616 /* classic formulation converts inequality to equality constraints via slack variables */ 617 ierr = VecAXPY(auglag->Ci, -1.0, auglag->Ps);CHKERRQ(ierr); 618 break; 619 case (TAO_ALMM_PHR): 620 /* PHR is based on Ci <= 0 while TAO defines Ci >= 0 so we hit it with a negative sign */ 621 ierr = VecScale(auglag->Ci, -1.0);CHKERRQ(ierr); 622 ierr = MatScale(auglag->Ai, -1.0);CHKERRQ(ierr); 623 break; 624 default: 625 break; 626 } 627 } 628 /* combine constraints into one vector */ 629 ierr = TaoALMMCombineDual_Private(tao, auglag->Ce, auglag->Ci, auglag->C);CHKERRQ(ierr); 630 PetscFunctionReturn(0); 631 } 632 633 /* 634 Lphr = f + 0.5*mu*[ (Ce + Ye/mu)^T (Ce + Ye/mu) + pmin(0, Ci + Yi/mu)^T pmin(0, Ci + Yi/mu)] 635 636 dLphr/dX = dF/dX + mu*[ (Ce + Ye/mu)^T Ae + pmin(0, Ci + Yi/mu)^T Ai] 637 638 dLphr/dS = 0 639 */ 640 static PetscErrorCode TaoALMMComputePHRLagAndGradient_Private(Tao tao) 641 { 642 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 643 PetscReal eq_norm=0.0, ineq_norm=0.0; 644 PetscErrorCode ierr; 645 646 PetscFunctionBegin; 647 ierr = TaoALMMEvaluateIterate_Private(tao, auglag->P);CHKERRQ(ierr); 648 if (tao->eq_constrained) { 649 /* Ce_work = mu*(Ce + Ye/mu) */ 650 ierr = VecWAXPY(auglag->Cework, 1.0/auglag->mu, auglag->Ye, auglag->Ce);CHKERRQ(ierr); 651 ierr = VecDot(auglag->Cework, auglag->Cework, &eq_norm);CHKERRQ(ierr); /* contribution to scalar Lagrangian */ 652 ierr = VecScale(auglag->Cework, auglag->mu);CHKERRQ(ierr); 653 /* dL/dX += mu*(Ce + Ye/mu)^T Ae */ 654 ierr = MatMultTransposeAdd(auglag->Ae, auglag->Cework, auglag->LgradX, auglag->LgradX);CHKERRQ(ierr); 655 } 656 if (tao->ineq_constrained) { 657 /* Ci_work = mu * pmax(0, Ci + Yi/mu) where pmax() is pointwise max() */ 658 ierr = VecWAXPY(auglag->Ciwork, 1.0/auglag->mu, auglag->Yi, auglag->Ci);CHKERRQ(ierr); 659 ierr = VecPointwiseMax(auglag->Ciwork, auglag->Cizero, auglag->Ciwork);CHKERRQ(ierr); 660 ierr = VecDot(auglag->Ciwork, auglag->Ciwork, &ineq_norm);CHKERRQ(ierr); /* contribution to scalar Lagrangian */ 661 /* dL/dX += mu * pmax(0, Ci + Yi/mu)^T Ai */ 662 ierr = VecScale(auglag->Ciwork, auglag->mu);CHKERRQ(ierr); 663 ierr = MatMultTransposeAdd(auglag->Ai, auglag->Ciwork, auglag->LgradX, auglag->LgradX);CHKERRQ(ierr); 664 /* dL/dS = 0 because there are no slacks in PHR */ 665 ierr = VecZeroEntries(auglag->LgradS);CHKERRQ(ierr); 666 } 667 /* combine gradient together */ 668 ierr = TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);CHKERRQ(ierr); 669 /* compute L = f + 0.5 * mu * [(Ce + Ye/mu)^T (Ce + Ye/mu) + pmax(0, Ci + Yi/mu)^T pmax(0, Ci + Yi/mu)] */ 670 auglag->Lval = auglag->fval + 0.5*auglag->mu*(eq_norm + ineq_norm);CHKERRQ(ierr); 671 PetscFunctionReturn(0); 672 } 673 674 /* 675 Lc = F + Ye^TCe + Yi^T(Ci - S) + 0.5*mu*[Ce^TCe + (Ci - S)^T(Ci - S)] 676 677 dLc/dX = dF/dX + Ye^TAe + Yi^TAi + 0.5*mu*[Ce^TAe + (Ci - S)^TAi] 678 679 dLc/dS = -[Yi + mu*(Ci - S)] 680 */ 681 static PetscErrorCode TaoALMMComputeAugLagAndGradient_Private(Tao tao) 682 { 683 TAO_ALMM *auglag = (TAO_ALMM*)tao->data; 684 PetscReal yeTce=0.0, yiTcims=0.0, ceTce=0.0, cimsTcims=0.0; 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 ierr = TaoALMMEvaluateIterate_Private(tao, auglag->P);CHKERRQ(ierr); 689 if (tao->eq_constrained) { 690 /* compute scalar contributions */ 691 ierr = VecDot(auglag->Ye, auglag->Ce, &yeTce);CHKERRQ(ierr); 692 ierr = VecDot(auglag->Ce, auglag->Ce, &ceTce);CHKERRQ(ierr); 693 /* dL/dX += ye^T Ae */ 694 ierr = MatMultTransposeAdd(auglag->Ae, auglag->Ye, auglag->LgradX, auglag->LgradX);CHKERRQ(ierr); 695 /* dL/dX += mu * ce^T Ae */ 696 ierr = MatMultTranspose(auglag->Ae, auglag->Ce, auglag->Xwork);CHKERRQ(ierr); 697 ierr = VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork);CHKERRQ(ierr); 698 } 699 if (tao->ineq_constrained) { 700 /* compute scalar contributions */ 701 ierr = VecDot(auglag->Yi, auglag->Ci, &yiTcims);CHKERRQ(ierr); 702 ierr = VecDot(auglag->Ci, auglag->Ci, &cimsTcims);CHKERRQ(ierr); 703 /* dL/dX += yi^T Ai */ 704 ierr = MatMultTransposeAdd(auglag->Ai, auglag->Yi, auglag->LgradX, auglag->LgradX);CHKERRQ(ierr); 705 /* dL/dX += mu * (ci - s)^T Ai */ 706 ierr = MatMultTranspose(auglag->Ai, auglag->Ci, auglag->Xwork);CHKERRQ(ierr); 707 ierr = VecAXPY(auglag->LgradX, auglag->mu, auglag->Xwork);CHKERRQ(ierr); 708 /* dL/dS = -[yi + mu*(ci - s)] */ 709 ierr = VecWAXPY(auglag->LgradS, auglag->mu, auglag->Ci, auglag->Yi);CHKERRQ(ierr); 710 ierr = VecScale(auglag->LgradS, -1.0);CHKERRQ(ierr); 711 } 712 /* combine gradient together */ 713 ierr = TaoALMMCombinePrimal_Private(tao, auglag->LgradX, auglag->LgradS, auglag->G);CHKERRQ(ierr); 714 /* compute L = f + ye^T ce + yi^T (ci - s) + 0.5*mu*||ce||^2 + 0.5*mu*||ci - s||^2 */ 715 auglag->Lval = auglag->fval + yeTce + yiTcims + 0.5*auglag->mu*(ceTce + cimsTcims); 716 PetscFunctionReturn(0); 717 } 718 719 PetscErrorCode TaoALMMSubsolverObjectiveAndGradient_Private(Tao tao, Vec P, PetscReal *Lval, Vec G, void *ctx) 720 { 721 TAO_ALMM *auglag = (TAO_ALMM*)ctx; 722 PetscErrorCode ierr; 723 724 PetscFunctionBegin; 725 ierr = VecCopy(P, auglag->P);CHKERRQ(ierr); 726 ierr = (*auglag->sub_obj)(auglag->parent);CHKERRQ(ierr); 727 ierr = VecCopy(auglag->G, G);CHKERRQ(ierr); 728 *Lval = auglag->Lval; 729 PetscFunctionReturn(0); 730 } 731