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->lamda);
 14:     VecDestroy(&lclP->lamda0);
 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:   PetscFree(tao->data);
 58:   return 0;
 59: }

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

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

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

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

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

105:   VecDuplicate(tao->constraints, &lclP->lamda);
106:   VecDuplicate(tao->constraints, &lclP->WL);
107:   VecDuplicate(tao->constraints, &lclP->lamda0);
108:   VecDuplicate(tao->constraints, &lclP->con1);

110:   VecSet(lclP->lamda,0.0);

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

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

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

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

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

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

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

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

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

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

229:   LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
230:   LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

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

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

242:   while (tao->reason == TAO_CONTINUE_ITERATING) {
243:     /* Call general purpose update function */
244:     if (tao->ops->update) {
245:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
246:     }
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) {
309:         PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %g -- using steepest descent\n",(double)descent);
310:       }

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

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

329:        ==>
330:        Check r'GL_U - rho*r'WU <= adec
331:     */

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

351:     LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
352:     LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
353:     LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

355:     /* We now minimize the augmented Lagrangian along the Newton direction */
356:     VecScale(lclP->r,-1.0);
357:     LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
358:     VecScale(lclP->r,-1.0);
359:     LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
360:     LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);

362:     lclP->recompute_jacobian_flag = PETSC_TRUE;

364:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
365:     TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
366:     TaoLineSearchSetFromOptions(tao->linesearch);
367:     TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
368:     if (lclP->verbose) {
369:       PetscPrintf(PETSC_COMM_WORLD,"Steplength = %g\n",(double)step);
370:     }

372:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
373:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
374:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

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

378:     /* Check convergence */
379:     VecNorm(lclP->GAugL, NORM_2, &mnorm);
380:     VecNorm(tao->constraints, NORM_2, &cnorm);
381:     TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
382:     TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
383:     (*tao->ops->convergencetest)(tao,tao->cnvP);
384:     if (tao->reason != TAO_CONTINUE_ITERATING) {
385:       break;
386:     }

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

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

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

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

403:       /* Store V and constraints */
404:       VecCopy(lclP->V, lclP->V1);
405:       VecCopy(tao->constraints,lclP->con1);

407:       /* Compute multipliers */
408:       /* p1 */
409:       VecSet(lclP->lamda,0.0); /*  Initial guess in CG */
410:       lclP->solve_type = LCL_ADJOINT1;
411:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
412:       if (tao->jacobian_state_pre) {
413:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
414:       } else {
415:         pset = pflag = PETSC_TRUE;
416:       }
417:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
418:       else symmetric = PETSC_FALSE;

420:       if (tao->jacobian_state_inv) {
421:         if (symmetric) {
422:           MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
423:         } else {
424:           MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
425:         }
426:       } else {
427:         if (symmetric) {
428:           KSPSolve(tao->ksp, lclP->GAugL_U,  lclP->lamda);
429:         } else {
430:           KSPSolveTranspose(tao->ksp, lclP->GAugL_U,  lclP->lamda);
431:         }
432:         KSPGetIterationNumber(tao->ksp,&its);
433:         tao->ksp_its+=its;
434:         tao->ksp_tot_its+=its;
435:       }
436:       MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
437:       VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);

439:       /* Compute the limited-memory quasi-newton direction */
440:       if (tao->niter > 0) {
441:         MatSolve(lclP->R,lclP->g1,lclP->s);
442:         VecDot(lclP->s,lclP->g1,&descent);
443:         if (descent <= 0) {
444:           if (lclP->verbose) {
445:             PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %g\n",(double)descent);
446:           }
447:           VecCopy(lclP->g1,lclP->s);
448:         }
449:       } else {
450:         VecCopy(lclP->g1,lclP->s);
451:       }
452:       VecScale(lclP->g1,-1.0);

454:       /* Recover the full space direction */
455:       MatMult(tao->jacobian_design,lclP->s,lclP->WU);
456:       /* VecSet(lclP->r,0.0); */ /*  Initial guess in CG */
457:       lclP->solve_type = LCL_FORWARD2;
458:       if (tao->jacobian_state_inv) {
459:         MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
460:       } else {
461:         KSPSolve(tao->ksp, lclP->WU, lclP->r);
462:         KSPGetIterationNumber(tao->ksp,&its);
463:         tao->ksp_its += its;
464:         tao->ksp_tot_its+=its;
465:       }

467:       /* We now minimize the augmented Lagrangian along the direction -r,s */
468:       VecScale(lclP->r, -1.0);
469:       LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
470:       VecScale(lclP->r, -1.0);
471:       lclP->recompute_jacobian_flag = PETSC_TRUE;

473:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
474:       TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT);
475:       TaoLineSearchSetFromOptions(tao->linesearch);
476:       TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
477:       if (lclP->verbose) {
478:         PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength =  %g\n",(double)step);
479:       }

481:       LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
482:       LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
483:       LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
484:       TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
485:       LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

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

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

492:       /* p2 */
493:       /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
494:       if (phase2_iter==0) {
495:         VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
496:       } else {
497:         VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
498:       }

500:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
501:       if (tao->jacobian_state_pre) {
502:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
503:       } else {
504:         pset = pflag = PETSC_TRUE;
505:       }
506:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
507:       else symmetric = PETSC_FALSE;

509:       lclP->solve_type = LCL_ADJOINT2;
510:       if (tao->jacobian_state_inv) {
511:         if (symmetric) {
512:           MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
513:         } else {
514:           MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
515:         }
516:       } else {
517:         if (symmetric) {
518:           KSPSolve(tao->ksp, lclP->GU,  lclP->lamda);
519:         } else {
520:           KSPSolveTranspose(tao->ksp, lclP->GU,  lclP->lamda);
521:         }
522:         KSPGetIterationNumber(tao->ksp,&its);
523:         tao->ksp_its += its;
524:         tao->ksp_tot_its += its;
525:       }

527:       MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g2);
528:       VecAXPY(lclP->g2,-1.0,lclP->GV);

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

532:       /* Update the quasi-newton approximation */
533:       MatLMVMUpdate(lclP->R,lclP->V,lclP->g2);
534:       /* 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 */

536:     }

538:     VecCopy(lclP->lamda,lclP->lamda0);

540:     /* Evaluate Function, Gradient, Constraints, and Jacobian */
541:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
542:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
543:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

545:     TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
546:     TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
547:     TaoComputeConstraints(tao,tao->solution, tao->constraints);

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

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

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

556:     /* Monitor convergence */
557:     tao->niter++;
558:     TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
559:     TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
560:     (*tao->ops->convergencetest)(tao,tao->cnvP);
561:   }
562:   return 0;
563: }

565: /*MC
566:  TAOLCL - linearly constrained lagrangian method for pde-constrained optimization

568: + -tao_lcl_eps1 - epsilon 1 tolerance
569: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
570: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
571: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
572: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
573: . -tao_lcl_verbose - Print verbose output if True
574: . -tao_lcl_tola - Tolerance for first forward solve
575: . -tao_lcl_tolb - Tolerance for first adjoint solve
576: . -tao_lcl_tolc - Tolerance for second forward solve
577: - -tao_lcl_told - Tolerance for second adjoint solve

579:   Level: beginner
580: M*/
581: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
582: {
583:   TAO_LCL        *lclP;
584:   const char     *morethuente_type = TAOLINESEARCHMT;

586:   tao->ops->setup = TaoSetup_LCL;
587:   tao->ops->solve = TaoSolve_LCL;
588:   tao->ops->view = TaoView_LCL;
589:   tao->ops->setfromoptions = TaoSetFromOptions_LCL;
590:   tao->ops->destroy = TaoDestroy_LCL;
591:   PetscNewLog(tao,&lclP);
592:   tao->data = (void*)lclP;

594:   /* Override default settings (unless already changed) */
595:   if (!tao->max_it_changed) tao->max_it = 200;
596:   if (!tao->catol_changed) tao->catol = 1.0e-4;
597:   if (!tao->gatol_changed) tao->gttol = 1.0e-4;
598:   if (!tao->grtol_changed) tao->gttol = 1.0e-4;
599:   if (!tao->gttol_changed) tao->gttol = 1.0e-4;
600:   lclP->rho0 = 1.0e-4;
601:   lclP->rhomax=1e5;
602:   lclP->eps1 = 1.0e-8;
603:   lclP->eps2 = 0.0;
604:   lclP->solve_type=2;
605:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
606:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
607:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
608:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
609:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

611:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
612:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
613:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
614:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
615:   KSPSetFromOptions(tao->ksp);

617:   MatCreate(((PetscObject)tao)->comm, &lclP->R);
618:   MatSetType(lclP->R, MATLMVMBFGS);
619:   return 0;
620: }

622: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
623: {
624:   Tao            tao = (Tao)ptr;
625:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;
626:   PetscBool      set,pset,flag,pflag,symmetric;
627:   PetscReal      cdotl;

629:   TaoComputeObjectiveAndGradient(tao,X,f,G);
630:   LCLScatter(lclP,G,lclP->GU,lclP->GV);
631:   if (lclP->recompute_jacobian_flag) {
632:     TaoComputeJacobianState(tao,X,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
633:     TaoComputeJacobianDesign(tao,X,tao->jacobian_design);
634:   }
635:   TaoComputeConstraints(tao,X, tao->constraints);
636:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
637:   if (tao->jacobian_state_pre) {
638:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
639:   } else {
640:     pset = pflag = PETSC_TRUE;
641:   }
642:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
643:   else symmetric = PETSC_FALSE;

645:   VecDot(lclP->lamda0, tao->constraints, &cdotl);
646:   lclP->lgn = *f - cdotl;

648:   /* Gradient of Lagrangian GL = G - J' * lamda */
649:   /*      WU = A' * WL
650:           WV = B' * WL */
651:   if (symmetric) {
652:     MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
653:   } else {
654:     MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
655:   }
656:   MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
657:   VecScale(lclP->GL_U,-1.0);
658:   VecScale(lclP->GL_V,-1.0);
659:   VecAXPY(lclP->GL_U,1.0,lclP->GU);
660:   VecAXPY(lclP->GL_V,1.0,lclP->GV);
661:   LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);

663:   f[0] = lclP->lgn;
664:   VecCopy(lclP->GL,G);
665:   return 0;
666: }

668: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
669: {
670:   Tao            tao = (Tao)ptr;
671:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;
672:   PetscReal      con2;
673:   PetscBool      flag,pflag,set,pset,symmetric;

675:   LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
676:   LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
677:   VecDot(tao->constraints,tao->constraints,&con2);
678:   lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;

680:   /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
681:   /*      WU = A' * c
682:           WV = B' * c */
683:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
684:   if (tao->jacobian_state_pre) {
685:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
686:   } else {
687:     pset = pflag = PETSC_TRUE;
688:   }
689:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
690:   else symmetric = PETSC_FALSE;

692:   if (symmetric) {
693:     MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
694:   } else {
695:     MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
696:   }

698:   MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
699:   VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
700:   VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
701:   LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);

703:   f[0] = lclP->aug;
704:   VecCopy(lclP->GAugL,G);
705:   return 0;
706: }

708: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
709: {
710:   VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
711:   VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
712:   VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
713:   VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
714:   return 0;

716: }
717: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
718: {
719:   VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
720:   VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
721:   VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
722:   VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
723:   return 0;

725: }