Actual source code: lcl.c

  1: #include <../src/tao/pde_constrained/impls/lcl/lcl.h>
  2: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);
  3: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);
  4: static PetscErrorCode LCLScatter(TAO_LCL *, Vec, Vec, Vec);
  5: static PetscErrorCode LCLGather(TAO_LCL *, Vec, Vec, Vec);

  7: static PetscErrorCode TaoDestroy_LCL(Tao tao)
  8: {
  9:   TAO_LCL *lclP = (TAO_LCL *)tao->data;

 11:   if (tao->setupcalled) {
 12:     MatDestroy(&lclP->R);
 13:     VecDestroy(&lclP->lambda);
 14:     VecDestroy(&lclP->lambda0);
 15:     VecDestroy(&lclP->WL);
 16:     VecDestroy(&lclP->W);
 17:     VecDestroy(&lclP->X0);
 18:     VecDestroy(&lclP->G0);
 19:     VecDestroy(&lclP->GL);
 20:     VecDestroy(&lclP->GAugL);
 21:     VecDestroy(&lclP->dbar);
 22:     VecDestroy(&lclP->U);
 23:     VecDestroy(&lclP->U0);
 24:     VecDestroy(&lclP->V);
 25:     VecDestroy(&lclP->V0);
 26:     VecDestroy(&lclP->V1);
 27:     VecDestroy(&lclP->GU);
 28:     VecDestroy(&lclP->GV);
 29:     VecDestroy(&lclP->GU0);
 30:     VecDestroy(&lclP->GV0);
 31:     VecDestroy(&lclP->GL_U);
 32:     VecDestroy(&lclP->GL_V);
 33:     VecDestroy(&lclP->GAugL_U);
 34:     VecDestroy(&lclP->GAugL_V);
 35:     VecDestroy(&lclP->GL_U0);
 36:     VecDestroy(&lclP->GL_V0);
 37:     VecDestroy(&lclP->GAugL_U0);
 38:     VecDestroy(&lclP->GAugL_V0);
 39:     VecDestroy(&lclP->DU);
 40:     VecDestroy(&lclP->DV);
 41:     VecDestroy(&lclP->WU);
 42:     VecDestroy(&lclP->WV);
 43:     VecDestroy(&lclP->g1);
 44:     VecDestroy(&lclP->g2);
 45:     VecDestroy(&lclP->con1);

 47:     VecDestroy(&lclP->r);
 48:     VecDestroy(&lclP->s);

 50:     ISDestroy(&tao->state_is);
 51:     ISDestroy(&tao->design_is);

 53:     VecScatterDestroy(&lclP->state_scatter);
 54:     VecScatterDestroy(&lclP->design_scatter);
 55:   }
 56:   MatDestroy(&lclP->R);
 57:   KSPDestroy(&tao->ksp);
 58:   PetscFree(tao->data);
 59:   return 0;
 60: }

 62: static PetscErrorCode TaoSetFromOptions_LCL(Tao tao, PetscOptionItems *PetscOptionsObject)
 63: {
 64:   TAO_LCL *lclP = (TAO_LCL *)tao->data;

 66:   PetscOptionsHeadBegin(PetscOptionsObject, "Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
 67:   PetscOptionsReal("-tao_lcl_eps1", "epsilon 1 tolerance", "", lclP->eps1, &lclP->eps1, NULL);
 68:   PetscOptionsReal("-tao_lcl_eps2", "epsilon 2 tolerance", "", lclP->eps2, &lclP->eps2, NULL);
 69:   PetscOptionsReal("-tao_lcl_rho0", "init value for rho", "", lclP->rho0, &lclP->rho0, NULL);
 70:   PetscOptionsReal("-tao_lcl_rhomax", "max value for rho", "", lclP->rhomax, &lclP->rhomax, NULL);
 71:   lclP->phase2_niter = 1;
 72:   PetscOptionsInt("-tao_lcl_phase2_niter", "Number of phase 2 iterations in LCL algorithm", "", lclP->phase2_niter, &lclP->phase2_niter, NULL);
 73:   lclP->verbose = PETSC_FALSE;
 74:   PetscOptionsBool("-tao_lcl_verbose", "Print verbose output", "", lclP->verbose, &lclP->verbose, NULL);
 75:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
 76:   PetscOptionsReal("-tao_lcl_tola", "Tolerance for first forward solve", "", lclP->tau[0], &lclP->tau[0], NULL);
 77:   PetscOptionsReal("-tao_lcl_tolb", "Tolerance for first adjoint solve", "", lclP->tau[1], &lclP->tau[1], NULL);
 78:   PetscOptionsReal("-tao_lcl_tolc", "Tolerance for second forward solve", "", lclP->tau[2], &lclP->tau[2], NULL);
 79:   PetscOptionsReal("-tao_lcl_told", "Tolerance for second adjoint solve", "", lclP->tau[3], &lclP->tau[3], NULL);
 80:   PetscOptionsHeadEnd();
 81:   TaoLineSearchSetFromOptions(tao->linesearch);
 82:   MatSetFromOptions(lclP->R);
 83:   return 0;
 84: }

 86: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
 87: {
 88:   return 0;
 89: }

 91: static PetscErrorCode TaoSetup_LCL(Tao tao)
 92: {
 93:   TAO_LCL *lclP = (TAO_LCL *)tao->data;
 94:   PetscInt lo, hi, nlocalstate, nlocaldesign;
 95:   IS       is_state, is_design;

 98:   VecDuplicate(tao->solution, &tao->gradient);
 99:   VecDuplicate(tao->solution, &tao->stepdirection);
100:   VecDuplicate(tao->solution, &lclP->W);
101:   VecDuplicate(tao->solution, &lclP->X0);
102:   VecDuplicate(tao->solution, &lclP->G0);
103:   VecDuplicate(tao->solution, &lclP->GL);
104:   VecDuplicate(tao->solution, &lclP->GAugL);

106:   VecDuplicate(tao->constraints, &lclP->lambda);
107:   VecDuplicate(tao->constraints, &lclP->WL);
108:   VecDuplicate(tao->constraints, &lclP->lambda0);
109:   VecDuplicate(tao->constraints, &lclP->con1);

111:   VecSet(lclP->lambda, 0.0);

113:   VecGetSize(tao->solution, &lclP->n);
114:   VecGetSize(tao->constraints, &lclP->m);

116:   VecCreate(((PetscObject)tao)->comm, &lclP->U);
117:   VecCreate(((PetscObject)tao)->comm, &lclP->V);
118:   ISGetLocalSize(tao->state_is, &nlocalstate);
119:   ISGetLocalSize(tao->design_is, &nlocaldesign);
120:   VecSetSizes(lclP->U, nlocalstate, lclP->m);
121:   VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m);
122:   VecSetType(lclP->U, ((PetscObject)(tao->solution))->type_name);
123:   VecSetType(lclP->V, ((PetscObject)(tao->solution))->type_name);
124:   VecSetFromOptions(lclP->U);
125:   VecSetFromOptions(lclP->V);
126:   VecDuplicate(lclP->U, &lclP->DU);
127:   VecDuplicate(lclP->U, &lclP->U0);
128:   VecDuplicate(lclP->U, &lclP->GU);
129:   VecDuplicate(lclP->U, &lclP->GU0);
130:   VecDuplicate(lclP->U, &lclP->GAugL_U);
131:   VecDuplicate(lclP->U, &lclP->GL_U);
132:   VecDuplicate(lclP->U, &lclP->GAugL_U0);
133:   VecDuplicate(lclP->U, &lclP->GL_U0);
134:   VecDuplicate(lclP->U, &lclP->WU);
135:   VecDuplicate(lclP->U, &lclP->r);
136:   VecDuplicate(lclP->V, &lclP->V0);
137:   VecDuplicate(lclP->V, &lclP->V1);
138:   VecDuplicate(lclP->V, &lclP->DV);
139:   VecDuplicate(lclP->V, &lclP->s);
140:   VecDuplicate(lclP->V, &lclP->GV);
141:   VecDuplicate(lclP->V, &lclP->GV0);
142:   VecDuplicate(lclP->V, &lclP->dbar);
143:   VecDuplicate(lclP->V, &lclP->GAugL_V);
144:   VecDuplicate(lclP->V, &lclP->GL_V);
145:   VecDuplicate(lclP->V, &lclP->GAugL_V0);
146:   VecDuplicate(lclP->V, &lclP->GL_V0);
147:   VecDuplicate(lclP->V, &lclP->WV);
148:   VecDuplicate(lclP->V, &lclP->g1);
149:   VecDuplicate(lclP->V, &lclP->g2);

151:   /* create scatters for state, design subvecs */
152:   VecGetOwnershipRange(lclP->U, &lo, &hi);
153:   ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state);
154:   VecGetOwnershipRange(lclP->V, &lo, &hi);
155:   if (0) {
156:     PetscInt sizeU, sizeV;
157:     VecGetSize(lclP->U, &sizeU);
158:     VecGetSize(lclP->V, &sizeV);
159:     PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV);
160:   }
161:   ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design);
162:   VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter);
163:   VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter);
164:   ISDestroy(&is_state);
165:   ISDestroy(&is_design);
166:   return 0;
167: }

169: static PetscErrorCode TaoSolve_LCL(Tao tao)
170: {
171:   TAO_LCL                     *lclP = (TAO_LCL *)tao->data;
172:   PetscInt                     phase2_iter, nlocal, its;
173:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
174:   PetscReal                    step      = 1.0, f, descent, aldescent;
175:   PetscReal                    cnorm, mnorm;
176:   PetscReal                    adec, r2, rGL_U, rWU;
177:   PetscBool                    set, pset, flag, pflag, symmetric;

179:   lclP->rho = lclP->rho0;
180:   VecGetLocalSize(lclP->U, &nlocal);
181:   VecGetLocalSize(lclP->V, &nlocal);
182:   MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m);
183:   MatLMVMAllocate(lclP->R, lclP->V, lclP->V);
184:   lclP->recompute_jacobian_flag = PETSC_TRUE;

186:   /* Scatter to U,V */
187:   LCLScatter(lclP, tao->solution, lclP->U, lclP->V);

189:   /* Evaluate Function, Gradient, Constraints, and Jacobian */
190:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
191:   TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
192:   TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design);
193:   TaoComputeConstraints(tao, tao->solution, tao->constraints);

195:   /* Scatter gradient to GU,GV */
196:   LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);

198:   /* Evaluate Lagrangian function and gradient */
199:   /* p0 */
200:   VecSet(lclP->lambda, 0.0); /*  Initial guess in CG */
201:   MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
202:   if (tao->jacobian_state_pre) {
203:     MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
204:   } else {
205:     pset = pflag = PETSC_TRUE;
206:   }
207:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
208:   else symmetric = PETSC_FALSE;

210:   lclP->solve_type = LCL_ADJOINT2;
211:   if (tao->jacobian_state_inv) {
212:     if (symmetric) {
213:       MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
214:     } else {
215:       MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
216:     }
217:   } else {
218:     KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
219:     if (symmetric) {
220:       KSPSolve(tao->ksp, lclP->GU, lclP->lambda);
221:     } else {
222:       KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda);
223:     }
224:     KSPGetIterationNumber(tao->ksp, &its);
225:     tao->ksp_its += its;
226:     tao->ksp_tot_its += its;
227:   }
228:   VecCopy(lclP->lambda, lclP->lambda0);
229:   LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);

231:   LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
232:   LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);

234:   /* Evaluate constraint norm */
235:   VecNorm(tao->constraints, NORM_2, &cnorm);
236:   VecNorm(lclP->GAugL, NORM_2, &mnorm);

238:   /* Monitor convergence */
239:   tao->reason = TAO_CONTINUE_ITERATING;
240:   TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
241:   TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
242:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

244:   while (tao->reason == TAO_CONTINUE_ITERATING) {
245:     /* Call general purpose update function */
246:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
247:     tao->ksp_its = 0;
248:     /* Compute a descent direction for the linearly constrained subproblem
249:        minimize f(u+du, v+dv)
250:        s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */

252:     /* Store the points around the linearization */
253:     VecCopy(lclP->U, lclP->U0);
254:     VecCopy(lclP->V, lclP->V0);
255:     VecCopy(lclP->GU, lclP->GU0);
256:     VecCopy(lclP->GV, lclP->GV0);
257:     VecCopy(lclP->GAugL_U, lclP->GAugL_U0);
258:     VecCopy(lclP->GAugL_V, lclP->GAugL_V0);
259:     VecCopy(lclP->GL_U, lclP->GL_U0);
260:     VecCopy(lclP->GL_V, lclP->GL_V0);

262:     lclP->aug0 = lclP->aug;
263:     lclP->lgn0 = lclP->lgn;

265:     /* Given the design variables, we need to project the current iterate
266:        onto the linearized constraint.  We choose to fix the design variables
267:        and solve the linear system for the state variables.  The resulting
268:        point is the Newton direction */

270:     /* Solve r = A\con */
271:     lclP->solve_type = LCL_FORWARD1;
272:     VecSet(lclP->r, 0.0); /*  Initial guess in CG */

274:     if (tao->jacobian_state_inv) {
275:       MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
276:     } else {
277:       KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
278:       KSPSolve(tao->ksp, tao->constraints, lclP->r);
279:       KSPGetIterationNumber(tao->ksp, &its);
280:       tao->ksp_its += its;
281:       tao->ksp_tot_its += tao->ksp_its;
282:     }

284:     /* Set design step direction dv to zero */
285:     VecSet(lclP->s, 0.0);

287:     /*
288:        Check sufficient descent for constraint merit function .5*||con||^2
289:        con' Ak r >= eps1 ||r||^(2+eps2)
290:     */

292:     /* Compute WU= Ak' * con */
293:     if (symmetric) {
294:       MatMult(tao->jacobian_state, tao->constraints, lclP->WU);
295:     } else {
296:       MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->WU);
297:     }
298:     /* Compute r * Ak' * con */
299:     VecDot(lclP->r, lclP->WU, &rWU);

301:     /* compute ||r||^(2+eps2) */
302:     VecNorm(lclP->r, NORM_2, &r2);
303:     r2   = PetscPowScalar(r2, 2.0 + lclP->eps2);
304:     adec = lclP->eps1 * r2;

306:     if (rWU < adec) {
307:       PetscInfo(tao, "Newton direction not descent for constraint, feasibility phase required\n");
308:       if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent);

310:       PetscInfo(tao, "Using steepest descent direction instead.\n");
311:       VecSet(lclP->r, 0.0);
312:       VecAXPY(lclP->r, -1.0, lclP->WU);
313:       VecDot(lclP->r, lclP->r, &rWU);
314:       VecNorm(lclP->r, NORM_2, &r2);
315:       r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
316:       VecDot(lclP->r, lclP->GAugL_U, &descent);
317:       adec = lclP->eps1 * r2;
318:     }

320:     /*
321:        Check descent for aug. lagrangian
322:        r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
323:           GL_U = GUk - Ak'*yk
324:           WU   = Ak'*con
325:                                          adec=eps1||r||^(2+eps2)

327:        ==>
328:        Check r'GL_U - rho*r'WU <= adec
329:     */

331:     VecDot(lclP->r, lclP->GL_U, &rGL_U);
332:     aldescent = rGL_U - lclP->rho * rWU;
333:     if (aldescent > -adec) {
334:       if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent);
335:       PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent);
336:       lclP->rho = (rGL_U - adec) / rWU;
337:       if (lclP->rho > lclP->rhomax) {
338:         lclP->rho = lclP->rhomax;
339:         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "rho=%g > rhomax, case not implemented.  Increase rhomax (-tao_lcl_rhomax)", (double)lclP->rho);
340:       }
341:       if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "  Increasing penalty parameter to %g\n", (double)lclP->rho);
342:       PetscInfo(tao, "  Increasing penalty parameter to %g\n", (double)lclP->rho);
343:     }

345:     LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);
346:     LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
347:     LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);

349:     /* We now minimize the augmented Lagrangian along the Newton direction */
350:     VecScale(lclP->r, -1.0);
351:     LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection);
352:     VecScale(lclP->r, -1.0);
353:     LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
354:     LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0);

356:     lclP->recompute_jacobian_flag = PETSC_TRUE;

358:     TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
359:     TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
360:     TaoLineSearchSetFromOptions(tao->linesearch);
361:     TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
362:     if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step);

364:     LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
365:     TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
366:     LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);

368:     LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);

370:     /* Check convergence */
371:     VecNorm(lclP->GAugL, NORM_2, &mnorm);
372:     VecNorm(tao->constraints, NORM_2, &cnorm);
373:     TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
374:     TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
375:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
376:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

378:     /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
379:     for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
380:       /* We now minimize the objective function starting from the fraction of
381:          the Newton point accepted by applying one step of a reduced-space
382:          method.  The optimization problem is:

384:          minimize f(u+du, v+dv)
385:          s. t.    A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)

387:          In particular, we have that
388:          du = -inv(A)*(Bdv + alpha g) */

390:       TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
391:       TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design);

393:       /* Store V and constraints */
394:       VecCopy(lclP->V, lclP->V1);
395:       VecCopy(tao->constraints, lclP->con1);

397:       /* Compute multipliers */
398:       /* p1 */
399:       VecSet(lclP->lambda, 0.0); /*  Initial guess in CG */
400:       lclP->solve_type = LCL_ADJOINT1;
401:       MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
402:       if (tao->jacobian_state_pre) {
403:         MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
404:       } else {
405:         pset = pflag = PETSC_TRUE;
406:       }
407:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
408:       else symmetric = PETSC_FALSE;

410:       if (tao->jacobian_state_inv) {
411:         if (symmetric) {
412:           MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda);
413:         } else {
414:           MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda);
415:         }
416:       } else {
417:         if (symmetric) {
418:           KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda);
419:         } else {
420:           KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda);
421:         }
422:         KSPGetIterationNumber(tao->ksp, &its);
423:         tao->ksp_its += its;
424:         tao->ksp_tot_its += its;
425:       }
426:       MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1);
427:       VecAXPY(lclP->g1, -1.0, lclP->GAugL_V);

429:       /* Compute the limited-memory quasi-newton direction */
430:       if (tao->niter > 0) {
431:         MatSolve(lclP->R, lclP->g1, lclP->s);
432:         VecDot(lclP->s, lclP->g1, &descent);
433:         if (descent <= 0) {
434:           if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent);
435:           VecCopy(lclP->g1, lclP->s);
436:         }
437:       } else {
438:         VecCopy(lclP->g1, lclP->s);
439:       }
440:       VecScale(lclP->g1, -1.0);

442:       /* Recover the full space direction */
443:       MatMult(tao->jacobian_design, lclP->s, lclP->WU);
444:       /* VecSet(lclP->r,0.0); */ /*  Initial guess in CG */
445:       lclP->solve_type = LCL_FORWARD2;
446:       if (tao->jacobian_state_inv) {
447:         MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r);
448:       } else {
449:         KSPSolve(tao->ksp, lclP->WU, lclP->r);
450:         KSPGetIterationNumber(tao->ksp, &its);
451:         tao->ksp_its += its;
452:         tao->ksp_tot_its += its;
453:       }

455:       /* We now minimize the augmented Lagrangian along the direction -r,s */
456:       VecScale(lclP->r, -1.0);
457:       LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection);
458:       VecScale(lclP->r, -1.0);
459:       lclP->recompute_jacobian_flag = PETSC_TRUE;

461:       TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
462:       TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
463:       TaoLineSearchSetFromOptions(tao->linesearch);
464:       TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
465:       if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength =  %g\n", (double)step);

467:       LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
468:       LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
469:       LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);
470:       TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
471:       LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);

473:       /* Compute the reduced gradient at the new point */

475:       TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
476:       TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design);

478:       /* p2 */
479:       /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
480:       if (phase2_iter == 0) {
481:         VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0);
482:       } else {
483:         VecAXPY(lclP->lambda, -lclP->rho, tao->constraints);
484:       }

486:       MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
487:       if (tao->jacobian_state_pre) {
488:         MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
489:       } else {
490:         pset = pflag = PETSC_TRUE;
491:       }
492:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
493:       else symmetric = PETSC_FALSE;

495:       lclP->solve_type = LCL_ADJOINT2;
496:       if (tao->jacobian_state_inv) {
497:         if (symmetric) {
498:           MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
499:         } else {
500:           MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
501:         }
502:       } else {
503:         if (symmetric) {
504:           KSPSolve(tao->ksp, lclP->GU, lclP->lambda);
505:         } else {
506:           KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda);
507:         }
508:         KSPGetIterationNumber(tao->ksp, &its);
509:         tao->ksp_its += its;
510:         tao->ksp_tot_its += its;
511:       }

513:       MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2);
514:       VecAXPY(lclP->g2, -1.0, lclP->GV);

516:       VecScale(lclP->g2, -1.0);

518:       /* Update the quasi-newton approximation */
519:       MatLMVMUpdate(lclP->R, lclP->V, lclP->g2);
520:       /* Use "-tao_ls_type gpcg -tao_ls_ftol 0 -tao_lmm_broyden_phi 0.0 -tao_lmm_scale_type scalar" to obtain agreement with MATLAB code */
521:     }

523:     VecCopy(lclP->lambda, lclP->lambda0);

525:     /* Evaluate Function, Gradient, Constraints, and Jacobian */
526:     TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
527:     LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
528:     LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);

530:     TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
531:     TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design);
532:     TaoComputeConstraints(tao, tao->solution, tao->constraints);

534:     LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);

536:     VecNorm(lclP->GAugL, NORM_2, &mnorm);

538:     /* Evaluate constraint norm */
539:     VecNorm(tao->constraints, NORM_2, &cnorm);

541:     /* Monitor convergence */
542:     tao->niter++;
543:     TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
544:     TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
545:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
546:   }
547:   return 0;
548: }

550: /*MC
551:  TAOLCL - linearly constrained lagrangian method for pde-constrained optimization

553: + -tao_lcl_eps1 - epsilon 1 tolerance
554: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
555: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
556: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
557: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
558: . -tao_lcl_verbose - Print verbose output if True
559: . -tao_lcl_tola - Tolerance for first forward solve
560: . -tao_lcl_tolb - Tolerance for first adjoint solve
561: . -tao_lcl_tolc - Tolerance for second forward solve
562: - -tao_lcl_told - Tolerance for second adjoint solve

564:   Level: beginner
565: M*/
566: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
567: {
568:   TAO_LCL    *lclP;
569:   const char *morethuente_type = TAOLINESEARCHMT;

571:   tao->ops->setup          = TaoSetup_LCL;
572:   tao->ops->solve          = TaoSolve_LCL;
573:   tao->ops->view           = TaoView_LCL;
574:   tao->ops->setfromoptions = TaoSetFromOptions_LCL;
575:   tao->ops->destroy        = TaoDestroy_LCL;
576:   PetscNew(&lclP);
577:   tao->data = (void *)lclP;

579:   /* Override default settings (unless already changed) */
580:   if (!tao->max_it_changed) tao->max_it = 200;
581:   if (!tao->catol_changed) tao->catol = 1.0e-4;
582:   if (!tao->gatol_changed) tao->gttol = 1.0e-4;
583:   if (!tao->grtol_changed) tao->gttol = 1.0e-4;
584:   if (!tao->gttol_changed) tao->gttol = 1.0e-4;
585:   lclP->rho0       = 1.0e-4;
586:   lclP->rhomax     = 1e5;
587:   lclP->eps1       = 1.0e-8;
588:   lclP->eps2       = 0.0;
589:   lclP->solve_type = 2;
590:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
591:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
592:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
593:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
594:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);

596:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao);
597:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
598:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
599:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
600:   KSPSetFromOptions(tao->ksp);

602:   MatCreate(((PetscObject)tao)->comm, &lclP->R);
603:   MatSetType(lclP->R, MATLMVMBFGS);
604:   return 0;
605: }

607: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
608: {
609:   Tao       tao  = (Tao)ptr;
610:   TAO_LCL  *lclP = (TAO_LCL *)tao->data;
611:   PetscBool set, pset, flag, pflag, symmetric;
612:   PetscReal cdotl;

614:   TaoComputeObjectiveAndGradient(tao, X, f, G);
615:   LCLScatter(lclP, G, lclP->GU, lclP->GV);
616:   if (lclP->recompute_jacobian_flag) {
617:     TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
618:     TaoComputeJacobianDesign(tao, X, tao->jacobian_design);
619:   }
620:   TaoComputeConstraints(tao, X, tao->constraints);
621:   MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
622:   if (tao->jacobian_state_pre) {
623:     MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
624:   } else {
625:     pset = pflag = PETSC_TRUE;
626:   }
627:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
628:   else symmetric = PETSC_FALSE;

630:   VecDot(lclP->lambda0, tao->constraints, &cdotl);
631:   lclP->lgn = *f - cdotl;

633:   /* Gradient of Lagrangian GL = G - J' * lambda */
634:   /*      WU = A' * WL
635:           WV = B' * WL */
636:   if (symmetric) {
637:     MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U);
638:   } else {
639:     MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U);
640:   }
641:   MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V);
642:   VecScale(lclP->GL_U, -1.0);
643:   VecScale(lclP->GL_V, -1.0);
644:   VecAXPY(lclP->GL_U, 1.0, lclP->GU);
645:   VecAXPY(lclP->GL_V, 1.0, lclP->GV);
646:   LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL);

648:   f[0] = lclP->lgn;
649:   VecCopy(lclP->GL, G);
650:   return 0;
651: }

653: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
654: {
655:   Tao       tao  = (Tao)ptr;
656:   TAO_LCL  *lclP = (TAO_LCL *)tao->data;
657:   PetscReal con2;
658:   PetscBool flag, pflag, set, pset, symmetric;

660:   LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao);
661:   LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V);
662:   VecDot(tao->constraints, tao->constraints, &con2);
663:   lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;

665:   /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
666:   /*      WU = A' * c
667:           WV = B' * c */
668:   MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
669:   if (tao->jacobian_state_pre) {
670:     MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
671:   } else {
672:     pset = pflag = PETSC_TRUE;
673:   }
674:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
675:   else symmetric = PETSC_FALSE;

677:   if (symmetric) {
678:     MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U);
679:   } else {
680:     MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U);
681:   }

683:   MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V);
684:   VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U);
685:   VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V);
686:   LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL);

688:   f[0] = lclP->aug;
689:   VecCopy(lclP->GAugL, G);
690:   return 0;
691: }

693: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
694: {
695:   VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
696:   VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
697:   VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
698:   VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
699:   return 0;
700: }
701: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
702: {
703:   VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
704:   VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
705:   VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
706:   VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
707:   return 0;
708: }