Actual source code: lcl.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <../src/tao/pde_constrained/impls/lcl/lcl.h>
  2: #include <../src/tao/matrix/lmvmmat.h>
  3: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
  4: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
  5: static PetscErrorCode LCLScatter(TAO_LCL*,Vec,Vec,Vec);
  6: static PetscErrorCode LCLGather(TAO_LCL*,Vec,Vec,Vec);

 10: static PetscErrorCode TaoDestroy_LCL(Tao tao)
 11: {
 12:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;

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

 52:     VecDestroy(&lclP->r);
 53:     VecDestroy(&lclP->s);

 55:     ISDestroy(&tao->state_is);
 56:     ISDestroy(&tao->design_is);

 58:     VecScatterDestroy(&lclP->state_scatter);
 59:     VecScatterDestroy(&lclP->design_scatter);
 60:   }
 61:   PetscFree(tao->data);
 62:   return(0);
 63: }

 67: static PetscErrorCode TaoSetFromOptions_LCL(PetscOptionItems *PetscOptionsObject,Tao tao)
 68: {
 69:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;

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

 94: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
 95: {
 96:   return 0;
 97: }

101: static PetscErrorCode TaoSetup_LCL(Tao tao)
102: {
103:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;
104:   PetscInt       lo, hi, nlocalstate, nlocaldesign;
106:   IS             is_state, is_design;

109:   if (!tao->state_is) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
110:   VecDuplicate(tao->solution, &tao->gradient);
111:   VecDuplicate(tao->solution, &tao->stepdirection);
112:   VecDuplicate(tao->solution, &lclP->W);
113:   VecDuplicate(tao->solution, &lclP->X0);
114:   VecDuplicate(tao->solution, &lclP->G0);
115:   VecDuplicate(tao->solution, &lclP->GL);
116:   VecDuplicate(tao->solution, &lclP->GAugL);

118:   VecDuplicate(tao->constraints, &lclP->lamda);
119:   VecDuplicate(tao->constraints, &lclP->WL);
120:   VecDuplicate(tao->constraints, &lclP->lamda0);
121:   VecDuplicate(tao->constraints, &lclP->con1);

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

125:   VecGetSize(tao->solution, &lclP->n);
126:   VecGetSize(tao->constraints, &lclP->m);

128:   VecCreate(((PetscObject)tao)->comm,&lclP->U);
129:   VecCreate(((PetscObject)tao)->comm,&lclP->V);
130:   ISGetLocalSize(tao->state_is,&nlocalstate);
131:   ISGetLocalSize(tao->design_is,&nlocaldesign);
132:   VecSetSizes(lclP->U,nlocalstate,lclP->m);
133:   VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m);
134:   VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name);
135:   VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name);
136:   VecSetFromOptions(lclP->U);
137:   VecSetFromOptions(lclP->V);
138:   VecDuplicate(lclP->U,&lclP->DU);
139:   VecDuplicate(lclP->U,&lclP->U0);
140:   VecDuplicate(lclP->U,&lclP->GU);
141:   VecDuplicate(lclP->U,&lclP->GU0);
142:   VecDuplicate(lclP->U,&lclP->GAugL_U);
143:   VecDuplicate(lclP->U,&lclP->GL_U);
144:   VecDuplicate(lclP->U,&lclP->GAugL_U0);
145:   VecDuplicate(lclP->U,&lclP->GL_U0);
146:   VecDuplicate(lclP->U,&lclP->WU);
147:   VecDuplicate(lclP->U,&lclP->r);
148:   VecDuplicate(lclP->V,&lclP->V0);
149:   VecDuplicate(lclP->V,&lclP->V1);
150:   VecDuplicate(lclP->V,&lclP->DV);
151:   VecDuplicate(lclP->V,&lclP->s);
152:   VecDuplicate(lclP->V,&lclP->GV);
153:   VecDuplicate(lclP->V,&lclP->GV0);
154:   VecDuplicate(lclP->V,&lclP->dbar);
155:   VecDuplicate(lclP->V,&lclP->GAugL_V);
156:   VecDuplicate(lclP->V,&lclP->GL_V);
157:   VecDuplicate(lclP->V,&lclP->GAugL_V0);
158:   VecDuplicate(lclP->V,&lclP->GL_V0);
159:   VecDuplicate(lclP->V,&lclP->WV);
160:   VecDuplicate(lclP->V,&lclP->g1);
161:   VecDuplicate(lclP->V,&lclP->g2);

163:   /* create scatters for state, design subvecs */
164:   VecGetOwnershipRange(lclP->U,&lo,&hi);
165:   ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state);
166:   VecGetOwnershipRange(lclP->V,&lo,&hi);
167:   if (0) {
168:     PetscInt sizeU,sizeV;
169:     VecGetSize(lclP->U,&sizeU);
170:     VecGetSize(lclP->V,&sizeV);
171:     PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV);
172:   }
173:   ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design);
174:   VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter);
175:   VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter);
176:   ISDestroy(&is_state);
177:   ISDestroy(&is_design);
178:   return(0);
179: }

183: static PetscErrorCode TaoSolve_LCL(Tao tao)
184: {
185:   TAO_LCL                      *lclP = (TAO_LCL*)tao->data;
186:   PetscInt                     phase2_iter,nlocal,its;
187:   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
188:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
189:   PetscReal                    step=1.0,f, descent, aldescent;
190:   PetscReal                    cnorm, mnorm;
191:   PetscReal                    adec,r2,rGL_U,rWU;
192:   PetscBool                    set,pset,flag,pflag,symmetric;
193:   PetscErrorCode               ierr;

196:   lclP->rho = lclP->rho0;
197:   VecGetLocalSize(lclP->U,&nlocal);
198:   VecGetLocalSize(lclP->V,&nlocal);
199:   MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R);
200:   MatLMVMAllocateVectors(lclP->R,lclP->V);
201:   lclP->recompute_jacobian_flag = PETSC_TRUE;

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

206:   /* Evaluate Function, Gradient, Constraints, and Jacobian */
207:   TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
208:   TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
209:   TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
210:   TaoComputeConstraints(tao,tao->solution, tao->constraints);

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

215:   /* Evaluate Lagrangian function and gradient */
216:   /* p0 */
217:   VecSet(lclP->lamda,0.0); /*  Initial guess in CG */
218:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
219:   if (tao->jacobian_state_pre) {
220:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
221:   } else {
222:     pset = pflag = PETSC_TRUE;
223:   }
224:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
225:   else symmetric = PETSC_FALSE;

227:   lclP->solve_type = LCL_ADJOINT2;
228:   if (tao->jacobian_state_inv) {
229:     if (symmetric) {
230:       MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
231:       MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
232:     }
233:   } else {
234:     KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
235:     if (symmetric) {
236:       KSPSolve(tao->ksp, lclP->GU,  lclP->lamda);
237:     } else {
238:       KSPSolveTranspose(tao->ksp, lclP->GU,  lclP->lamda);
239:     }
240:     KSPGetIterationNumber(tao->ksp,&its);
241:     tao->ksp_its+=its;
242:     tao->ksp_tot_its+=its;
243:   }
244:   VecCopy(lclP->lamda,lclP->lamda0);
245:   LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);

247:   LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
248:   LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

250:   /* Evaluate constraint norm */
251:   VecNorm(tao->constraints, NORM_2, &cnorm);
252:   VecNorm(lclP->GAugL, NORM_2, &mnorm);

254:   /* Monitor convergence */
255:   TaoMonitor(tao, tao->niter,f,mnorm,cnorm,step,&reason);

257:   while (reason == TAO_CONTINUE_ITERATING) {
258:     tao->ksp_its=0;
259:     /* Compute a descent direction for the linearly constrained subproblem
260:        minimize f(u+du, v+dv)
261:        s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */

263:     /* Store the points around the linearization */
264:     VecCopy(lclP->U, lclP->U0);
265:     VecCopy(lclP->V, lclP->V0);
266:     VecCopy(lclP->GU,lclP->GU0);
267:     VecCopy(lclP->GV,lclP->GV0);
268:     VecCopy(lclP->GAugL_U,lclP->GAugL_U0);
269:     VecCopy(lclP->GAugL_V,lclP->GAugL_V0);
270:     VecCopy(lclP->GL_U,lclP->GL_U0);
271:     VecCopy(lclP->GL_V,lclP->GL_V0);

273:     lclP->aug0 = lclP->aug;
274:     lclP->lgn0 = lclP->lgn;

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

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

285:     if (tao->jacobian_state_inv) {
286:       MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
287:     } else {
288:       KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
289:       KSPSolve(tao->ksp, tao->constraints,  lclP->r);
290:       KSPGetIterationNumber(tao->ksp,&its);
291:       tao->ksp_its+=its;
292:       tao->ksp_tot_its+=tao->ksp_its;
293:     }

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

298:     /*
299:        Check sufficient descent for constraint merit function .5*||con||^2
300:        con' Ak r >= eps1 ||r||^(2+eps2)
301:     */

303:     /* Compute WU= Ak' * con */
304:     if (symmetric)  {
305:       MatMult(tao->jacobian_state,tao->constraints,lclP->WU);
306:     } else {
307:       MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU);
308:     }
309:     /* Compute r * Ak' * con */
310:     VecDot(lclP->r,lclP->WU,&rWU);

312:     /* compute ||r||^(2+eps2) */
313:     VecNorm(lclP->r,NORM_2,&r2);
314:     r2 = PetscPowScalar(r2,2.0+lclP->eps2);
315:     adec = lclP->eps1 * r2;

317:     if (rWU < adec) {
318:       PetscInfo(tao,"Newton direction not descent for constraint, feasibility phase required\n");
319:       if (lclP->verbose) {
320:         PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %g -- using steepest descent\n",(double)descent);
321:       }

323:       PetscInfo(tao,"Using steepest descent direction instead.\n");
324:       VecSet(lclP->r,0.0);
325:       VecAXPY(lclP->r,-1.0,lclP->WU);
326:       VecDot(lclP->r,lclP->r,&rWU);
327:       VecNorm(lclP->r,NORM_2,&r2);
328:       r2 = PetscPowScalar(r2,2.0+lclP->eps2);
329:       VecDot(lclP->r,lclP->GAugL_U,&descent);
330:       adec = lclP->eps1 * r2;
331:     }


334:     /*
335:        Check descent for aug. lagrangian
336:        r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
337:           GL_U = GUk - Ak'*yk
338:           WU   = Ak'*con
339:                                          adec=eps1||r||^(2+eps2)

341:        ==>
342:        Check r'GL_U - rho*r'WU <= adec
343:     */

345:     VecDot(lclP->r,lclP->GL_U,&rGL_U);
346:     aldescent =  rGL_U - lclP->rho*rWU;
347:     if (aldescent > -adec) {
348:       if (lclP->verbose) {
349:         PetscPrintf(PETSC_COMM_WORLD," Newton direction not descent for augmented Lagrangian: %g",(double)aldescent);
350:       }
351:       PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %g\n",(double)aldescent);
352:       lclP->rho =  (rGL_U - adec)/rWU;
353:       if (lclP->rho > lclP->rhomax) {
354:         lclP->rho = lclP->rhomax;
355:         SETERRQ1(PETSC_COMM_WORLD,0,"rho=%g > rhomax, case not implemented.  Increase rhomax (-tao_lcl_rhomax)",(double)lclP->rho);
356:       }
357:       if (lclP->verbose) {
358:         PetscPrintf(PETSC_COMM_WORLD,"  Increasing penalty parameter to %g\n",(double)lclP->rho);
359:       }
360:       PetscInfo1(tao,"  Increasing penalty parameter to %g\n",(double)lclP->rho);
361:     }

363:     LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
364:     LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
365:     LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

367:     /* We now minimize the augmented Lagrangian along the Newton direction */
368:     VecScale(lclP->r,-1.0);
369:     LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
370:     VecScale(lclP->r,-1.0);
371:     LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
372:     LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);

374:     lclP->recompute_jacobian_flag = PETSC_TRUE;

376:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
377:     TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
378:     TaoLineSearchSetFromOptions(tao->linesearch);
379:     TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
380:     if (lclP->verbose) {
381:       PetscPrintf(PETSC_COMM_WORLD,"Steplength = %g\n",(double)step);
382:     }

384:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
385:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
386:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

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

390:     /* Check convergence */
391:     VecNorm(lclP->GAugL, NORM_2, &mnorm);
392:     VecNorm(tao->constraints, NORM_2, &cnorm);
393:     TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step,&reason);
394:     if (reason != TAO_CONTINUE_ITERATING){
395:       break;
396:     }

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

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

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

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

413:       /* Store V and constraints */
414:       VecCopy(lclP->V, lclP->V1);
415:       VecCopy(tao->constraints,lclP->con1);

417:       /* Compute multipliers */
418:       /* p1 */
419:       VecSet(lclP->lamda,0.0); /*  Initial guess in CG */
420:       lclP->solve_type = LCL_ADJOINT1;
421:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
422:       if (tao->jacobian_state_pre) {
423:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
424:       } else {
425:         pset = pflag = PETSC_TRUE;
426:       }
427:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
428:       else symmetric = PETSC_FALSE;

430:       if (tao->jacobian_state_inv) {
431:         if (symmetric) {
432:           MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
433:         } else {
434:           MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
435:         }
436:       } else {
437:         if (symmetric) {
438:           KSPSolve(tao->ksp, lclP->GAugL_U,  lclP->lamda);
439:         } else {
440:           KSPSolveTranspose(tao->ksp, lclP->GAugL_U,  lclP->lamda);
441:         }
442:         KSPGetIterationNumber(tao->ksp,&its);
443:         tao->ksp_its+=its;
444:         tao->ksp_tot_its+=its;
445:       }
446:       MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
447:       VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);

449:       /* Compute the limited-memory quasi-newton direction */
450:       if (tao->niter > 0) {
451:         MatLMVMSolve(lclP->R,lclP->g1,lclP->s);
452:         VecDot(lclP->s,lclP->g1,&descent);
453:         if (descent <= 0) {
454:           if (lclP->verbose) {
455:             PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %g\n",(double)descent);
456:           }
457:           VecCopy(lclP->g1,lclP->s);
458:         }
459:       } else {
460:         VecCopy(lclP->g1,lclP->s);
461:       }
462:       VecScale(lclP->g1,-1.0);

464:       /* Recover the full space direction */
465:       MatMult(tao->jacobian_design,lclP->s,lclP->WU);
466:       /* VecSet(lclP->r,0.0); */ /*  Initial guess in CG */
467:       lclP->solve_type = LCL_FORWARD2;
468:       if (tao->jacobian_state_inv) {
469:         MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
470:       } else {
471:         KSPSolve(tao->ksp, lclP->WU, lclP->r);
472:         KSPGetIterationNumber(tao->ksp,&its);
473:         tao->ksp_its += its;
474:         tao->ksp_tot_its+=its;
475:       }

477:       /* We now minimize the augmented Lagrangian along the direction -r,s */
478:       VecScale(lclP->r, -1.0);
479:       LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
480:       VecScale(lclP->r, -1.0);
481:       lclP->recompute_jacobian_flag = PETSC_TRUE;

483:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
484:       TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT);
485:       TaoLineSearchSetFromOptions(tao->linesearch);
486:       TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
487:       if (lclP->verbose){
488:         PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength =  %g\n",(double)step);
489:       }

491:       LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
492:       LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
493:       LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
494:       TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
495:       LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

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

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

502:       /* p2 */
503:       /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
504:       if (phase2_iter==0){
505:         VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
506:       } else {
507:         VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
508:       }

510:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
511:       if (tao->jacobian_state_pre) {
512:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
513:       } else {
514:         pset = pflag = PETSC_TRUE;
515:       }
516:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
517:       else symmetric = PETSC_FALSE;

519:       lclP->solve_type = LCL_ADJOINT2;
520:       if (tao->jacobian_state_inv) {
521:         if (symmetric) {
522:           MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
523:         } else {
524:           MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
525:         }
526:       } else {
527:         if (symmetric) {
528:           KSPSolve(tao->ksp, lclP->GU,  lclP->lamda);
529:         } else {
530:           KSPSolveTranspose(tao->ksp, lclP->GU,  lclP->lamda);
531:         }
532:         KSPGetIterationNumber(tao->ksp,&its);
533:         tao->ksp_its += its;
534:         tao->ksp_tot_its += its;
535:       }

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

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

542:       /* Update the quasi-newton approximation */
543:       if (phase2_iter >= 0){
544:         MatLMVMSetPrev(lclP->R,lclP->V1,lclP->g1);
545:       }
546:       MatLMVMUpdate(lclP->R,lclP->V,lclP->g2);
547:       /* 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 */

549:     }

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

553:     /* Evaluate Function, Gradient, Constraints, and Jacobian */
554:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
555:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
556:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

558:     TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
559:     TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
560:     TaoComputeConstraints(tao,tao->solution, tao->constraints);

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

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

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

569:     /* Monitor convergence */
570:     tao->niter++;
571:     TaoMonitor(tao, tao->niter,f,mnorm,cnorm,step,&reason);
572:   }
573:   MatDestroy(&lclP->R);
574:   return(0);
575: }

577: /*MC
578:  TAOLCL - linearly constrained lagrangian method for pde-constrained optimization

580: + -tao_lcl_eps1 - epsilon 1 tolerance
581: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
582: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
583: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
584: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
585: . -tao_lcl_verbose - Print verbose output if True
586: . -tao_lcl_tola - Tolerance for first forward solve
587: . -tao_lcl_tolb - Tolerance for first adjoint solve
588: . -tao_lcl_tolc - Tolerance for second forward solve
589: - -tao_lcl_told - Tolerance for second adjoint solve

591:   Level: beginner
592: M*/
595: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
596: {
597:   TAO_LCL        *lclP;
599:   const char     *morethuente_type = TAOLINESEARCHMT;

602:   tao->ops->setup = TaoSetup_LCL;
603:   tao->ops->solve = TaoSolve_LCL;
604:   tao->ops->view = TaoView_LCL;
605:   tao->ops->setfromoptions = TaoSetFromOptions_LCL;
606:   tao->ops->destroy = TaoDestroy_LCL;
607:   PetscNewLog(tao,&lclP);
608:   tao->data = (void*)lclP;

610:   /* Override default settings (unless already changed) */
611:   if (!tao->max_it_changed) tao->max_it = 200;
612:   if (!tao->catol_changed) tao->catol = 1.0e-4;
613:   if (!tao->gatol_changed) tao->gttol = 1.0e-4;
614:   if (!tao->grtol_changed) tao->gttol = 1.0e-4;
615:   if (!tao->gttol_changed) tao->gttol = 1.0e-4;
616:   lclP->rho0 = 1.0e-4;
617:   lclP->rhomax=1e5;
618:   lclP->eps1 = 1.0e-8;
619:   lclP->eps2 = 0.0;
620:   lclP->solve_type=2;
621:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
622:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
623:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
624:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

626:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
627:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
628:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
629:   KSPSetFromOptions(tao->ksp);
630:   return(0);
631: }

635: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
636: {
637:   Tao            tao = (Tao)ptr;
638:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;
639:   PetscBool      set,pset,flag,pflag,symmetric;
640:   PetscReal      cdotl;

644:   TaoComputeObjectiveAndGradient(tao,X,f,G);
645:   LCLScatter(lclP,G,lclP->GU,lclP->GV);
646:   if (lclP->recompute_jacobian_flag) {
647:     TaoComputeJacobianState(tao,X,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
648:     TaoComputeJacobianDesign(tao,X,tao->jacobian_design);
649:   }
650:   TaoComputeConstraints(tao,X, tao->constraints);
651:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
652:   if (tao->jacobian_state_pre) {
653:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
654:   } else {
655:     pset = pflag = PETSC_TRUE;
656:   }
657:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
658:   else symmetric = PETSC_FALSE;

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

663:   /* Gradient of Lagrangian GL = G - J' * lamda */
664:   /*      WU = A' * WL
665:           WV = B' * WL */
666:   if (symmetric) {
667:     MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
668:   } else {
669:     MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
670:   }
671:   MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
672:   VecScale(lclP->GL_U,-1.0);
673:   VecScale(lclP->GL_V,-1.0);
674:   VecAXPY(lclP->GL_U,1.0,lclP->GU);
675:   VecAXPY(lclP->GL_V,1.0,lclP->GV);
676:   LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);

678:   f[0] = lclP->lgn;
679:   VecCopy(lclP->GL,G);
680:   return(0);
681: }

685: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
686: {
687:   Tao            tao = (Tao)ptr;
688:   TAO_LCL        *lclP = (TAO_LCL*)tao->data;
689:   PetscReal      con2;
690:   PetscBool      flag,pflag,set,pset,symmetric;

694:   LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
695:   LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
696:   VecDot(tao->constraints,tao->constraints,&con2);
697:   lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;

699:   /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
700:   /*      WU = A' * c
701:           WV = B' * c */
702:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
703:   if (tao->jacobian_state_pre) {
704:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
705:   } else {
706:     pset = pflag = PETSC_TRUE;
707:   }
708:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
709:   else symmetric = PETSC_FALSE;

711:   if (symmetric) {
712:     MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
713:   } else {
714:     MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
715:   }

717:   MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
718:   VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
719:   VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
720:   LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);

722:   f[0] = lclP->aug;
723:   VecCopy(lclP->GAugL,G);
724:   return(0);
725: }

729: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
730: {
733:   VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
734:   VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
735:   VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
736:   VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
737:   return(0);

739: }
742: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
743: {
746:   VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
747:   VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
748:   VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
749:   VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
750:   return(0);

752: }