Actual source code: lcl.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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;

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

 49:     VecDestroy(&lclP->r);
 50:     VecDestroy(&lclP->s);

 52:     ISDestroy(&tao->state_is);
 53:     ISDestroy(&tao->design_is);

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

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

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

 89: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
 90: {
 91:   return 0;
 92: }

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

102:   if (!tao->state_is) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
103:   VecDuplicate(tao->solution, &tao->gradient);
104:   VecDuplicate(tao->solution, &tao->stepdirection);
105:   VecDuplicate(tao->solution, &lclP->W);
106:   VecDuplicate(tao->solution, &lclP->X0);
107:   VecDuplicate(tao->solution, &lclP->G0);
108:   VecDuplicate(tao->solution, &lclP->GL);
109:   VecDuplicate(tao->solution, &lclP->GAugL);

111:   VecDuplicate(tao->constraints, &lclP->lamda);
112:   VecDuplicate(tao->constraints, &lclP->WL);
113:   VecDuplicate(tao->constraints, &lclP->lamda0);
114:   VecDuplicate(tao->constraints, &lclP->con1);

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

118:   VecGetSize(tao->solution, &lclP->n);
119:   VecGetSize(tao->constraints, &lclP->m);

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

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

174: static PetscErrorCode TaoSolve_LCL(Tao tao)
175: {
176:   TAO_LCL                      *lclP = (TAO_LCL*)tao->data;
177:   PetscInt                     phase2_iter,nlocal,its;
178:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
179:   PetscReal                    step=1.0,f, descent, aldescent;
180:   PetscReal                    cnorm, mnorm;
181:   PetscReal                    adec,r2,rGL_U,rWU;
182:   PetscBool                    set,pset,flag,pflag,symmetric;
183:   PetscErrorCode               ierr;

186:   lclP->rho = lclP->rho0;
187:   VecGetLocalSize(lclP->U,&nlocal);
188:   VecGetLocalSize(lclP->V,&nlocal);
189:   MatSetSizes(lclP->R, nlocal, nlocal, lclP->n-lclP->m, lclP->n-lclP->m);
190:   MatLMVMAllocate(lclP->R,lclP->V,lclP->V);
191:   lclP->recompute_jacobian_flag = PETSC_TRUE;

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

196:   /* Evaluate Function, Gradient, Constraints, and Jacobian */
197:   TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
198:   TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
199:   TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
200:   TaoComputeConstraints(tao,tao->solution, tao->constraints);

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

205:   /* Evaluate Lagrangian function and gradient */
206:   /* p0 */
207:   VecSet(lclP->lamda,0.0); /*  Initial guess in CG */
208:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
209:   if (tao->jacobian_state_pre) {
210:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
211:   } else {
212:     pset = pflag = PETSC_TRUE;
213:   }
214:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
215:   else symmetric = PETSC_FALSE;

217:   lclP->solve_type = LCL_ADJOINT2;
218:   if (tao->jacobian_state_inv) {
219:     if (symmetric) {
220:       MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
221:       MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
222:     }
223:   } else {
224:     KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
225:     if (symmetric) {
226:       KSPSolve(tao->ksp, lclP->GU,  lclP->lamda);
227:     } else {
228:       KSPSolveTranspose(tao->ksp, lclP->GU,  lclP->lamda);
229:     }
230:     KSPGetIterationNumber(tao->ksp,&its);
231:     tao->ksp_its+=its;
232:     tao->ksp_tot_its+=its;
233:   }
234:   VecCopy(lclP->lamda,lclP->lamda0);
235:   LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);

237:   LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
238:   LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

240:   /* Evaluate constraint norm */
241:   VecNorm(tao->constraints, NORM_2, &cnorm);
242:   VecNorm(lclP->GAugL, NORM_2, &mnorm);

244:   /* Monitor convergence */
245:   tao->reason = TAO_CONTINUE_ITERATING;
246:   TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
247:   TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
248:   (*tao->ops->convergencetest)(tao,tao->cnvP);

250:   while (tao->reason == TAO_CONTINUE_ITERATING) {
251:     /* Call general purpose update function */
252:     if (tao->ops->update) {
253:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
254:     }
255:     tao->ksp_its=0;
256:     /* Compute a descent direction for the linearly constrained subproblem
257:        minimize f(u+du, v+dv)
258:        s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */

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

270:     lclP->aug0 = lclP->aug;
271:     lclP->lgn0 = lclP->lgn;

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

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

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

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

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

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

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

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

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


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

338:        ==>
339:        Check r'GL_U - rho*r'WU <= adec
340:     */

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

360:     LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
361:     LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
362:     LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);

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

371:     lclP->recompute_jacobian_flag = PETSC_TRUE;

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

381:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
382:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
383:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

545:     }

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

549:     /* Evaluate Function, Gradient, Constraints, and Jacobian */
550:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
551:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
552:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);

554:     TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
555:     TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
556:     TaoComputeConstraints(tao,tao->solution, tao->constraints);

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

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

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

565:     /* Monitor convergence */
566:     tao->niter++;
567:     TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
568:     TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
569:     (*tao->ops->convergencetest)(tao,tao->cnvP);
570:   }
571:   return(0);
572: }

574: /*MC
575:  TAOLCL - linearly constrained lagrangian method for pde-constrained optimization

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

588:   Level: beginner
589: M*/
590: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
591: {
592:   TAO_LCL        *lclP;
594:   const char     *morethuente_type = TAOLINESEARCHMT;

597:   tao->ops->setup = TaoSetup_LCL;
598:   tao->ops->solve = TaoSolve_LCL;
599:   tao->ops->view = TaoView_LCL;
600:   tao->ops->setfromoptions = TaoSetFromOptions_LCL;
601:   tao->ops->destroy = TaoDestroy_LCL;
602:   PetscNewLog(tao,&lclP);
603:   tao->data = (void*)lclP;

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

622:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
623:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
624:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
625:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
626:   KSPSetFromOptions(tao->ksp);

628:   MatCreate(((PetscObject)tao)->comm, &lclP->R);
629:   MatSetType(lclP->R, MATLMVMBFGS);
630:   return(0);
631: }

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

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

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

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

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

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

690:   LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
691:   LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
692:   VecDot(tao->constraints,tao->constraints,&con2);
693:   lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;

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

707:   if (symmetric) {
708:     MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
709:   } else {
710:     MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
711:   }

713:   MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
714:   VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
715:   VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
716:   LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);

718:   f[0] = lclP->aug;
719:   VecCopy(lclP->GAugL,G);
720:   return(0);
721: }

723: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
724: {
727:   VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
728:   VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
729:   VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
730:   VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
731:   return(0);

733: }
734: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
735: {
738:   VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
739:   VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
740:   VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
741:   VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
742:   return(0);

744: }