Actual source code: lcl.c

  1: #include "lcl.h"
  2: #include "src/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(TaoSolver tao)
 11: {
 12:   TAO_LCL *lclP = (TAO_LCL*)tao->data;
 15:   if (tao->setupcalled) {
 16:     MatDestroy(&lclP->R); 
 17:     
 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:   tao->data = PETSC_NULL;

 64:   return(0);
 65: }

 69: static PetscErrorCode TaoSetFromOptions_LCL(TaoSolver tao)
 70: {
 71:   /*TAO_LCL *lclP = (TAO_LCL*)tao->data;*/
 74:   TaoLineSearchSetFromOptions(tao->linesearch); 
 75:   return(0);
 76: }

 80: static PetscErrorCode TaoView_LCL(TaoSolver tao, PetscViewer viewer)
 81: {
 82:   /*
 83:   TAO_LCL *lclP = (TAO_LCL*)tao->data;
 86:   return(0);
 87:   */
 88:   return 0;
 89: }

 93: static PetscErrorCode TaoSetup_LCL(TaoSolver tao)
 94: {
 95:   TAO_LCL *lclP = (TAO_LCL*)tao->data;
 96:   PetscInt lo, hi, nlocalstate, nlocaldesign;
 98:   PetscBool flag;
 99:   IS is_state, is_design;
101:   /* Check for state/design IS */
102:   if (!tao->state_is) {
103:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
104:   }
105:   VecDuplicate(tao->solution, &tao->gradient); 
106:   VecDuplicate(tao->solution, &tao->stepdirection); 
107:   VecDuplicate(tao->solution, &lclP->W); 
108:   VecDuplicate(tao->solution, &lclP->X0); 
109:   VecDuplicate(tao->solution, &lclP->G0); 
110:   VecDuplicate(tao->solution, &lclP->GL); 
111:   VecDuplicate(tao->solution, &lclP->GAugL); 
112:   
113:   VecDuplicate(tao->constraints, &lclP->lamda); 
114:   VecDuplicate(tao->constraints, &lclP->WL); 
115:   VecDuplicate(tao->constraints, &lclP->lamda0); 
116:   VecDuplicate(tao->constraints, &lclP->con1); 

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

120:   VecGetSize(tao->solution, &lclP->n); 
121:   VecGetSize(tao->constraints, &lclP->m); 


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

160:   

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


179:   lclP->phase2_niter = 1;
180:   PetscOptionsInt("-tao_lcl_phase2_niter","Number of phase 2 iterations in LCL algorithm","",lclP->phase2_niter,&lclP->phase2_niter,&flag); 
181:   lclP->verbose = PETSC_FALSE;
182:   PetscOptionsBool("-tao_lcl_verbose","Print verbose output","",lclP->verbose,&lclP->verbose,&flag); 
183:   lclP->tola = 1e-4;
184:   PetscOptionsReal("-tao_lcl_tola","Tolerance for first forward solve","",lclP->tola,&lclP->tola,&flag); 
185:   lclP->tolb = 1e-4;
186:   PetscOptionsReal("-tao_lcl_tolb","Tolerance for first adjoint solve","",lclP->tolb,&lclP->tolb,&flag); 
187:   lclP->tolc = 1e-4;
188:   PetscOptionsReal("-tao_lcl_tolc","Tolerance for second forward solve","",lclP->tolc,&lclP->tolc,&flag); 
189:   lclP->told = 1e-4;
190:   PetscOptionsReal("-tao_lcl_told","Tolerance for second adjoint solve","",lclP->told,&lclP->told,&flag); 


193:   return(0);
194: }

198: static PetscErrorCode TaoSolve_LCL(TaoSolver tao)
199: {
200:   TAO_LCL *lclP = (TAO_LCL*)tao->data;
201:   PetscInt iter=0,phase2_iter,nlocal,its;
202:   TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
203:   TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
204:   PetscReal step=1.0,f, descent,con2;
205:   PetscReal cnorm, mnorm;
206:   PetscReal adec,r2,rGL_U,rWU;
207:   PetscBool set,pset,flag,pflag,symmetric;

211:   VecGetLocalSize(lclP->U,&nlocal); 
212:   VecGetLocalSize(lclP->V,&nlocal); 
213:   MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R); 
214:   MatLMVMAllocateVectors(lclP->R,lclP->V); 
215:   lclP->rho = 1.0e-4;
216:   lclP->recompute_jacobian_flag = PETSC_TRUE;
217:   
218:   /* Scatter to U,V */
219:   LCLScatter(lclP,tao->solution,lclP->U,lclP->V); 
220:   
221:   /* Evaluate Function, Gradient, Constraints, and Jacobian */
222:   TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); 
223:   TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag); 
224:   TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design); 
225:   TaoComputeConstraints(tao,tao->solution, tao->constraints); 

227:   
228:   /* Scatter gradient to GU,GV */
229:   LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); 
230:   

232:   /* Evaluate Lagrangian function and gradient */
233:   /* p0 */
234:   VecSet(lclP->lamda,0.0);  /*  Initial guess in CG */
235:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); 
236:   if (tao->jacobian_state_pre) {
237:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); 
238:   } else {
239:     pset = pflag = PETSC_TRUE;
240:   }
241:   if (set && pset && flag && pflag)
242:     symmetric = PETSC_TRUE;
243:   else
244:     symmetric = PETSC_FALSE;

246:   lclP->solve_type = LCL_ADJOINT2;
247:   if (tao->jacobian_state_inv) {
248:     if (symmetric) {
249:       MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);  } else {
250:       MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda); 
251:     }
252:   } else {
253:     KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag); 
254:     if (symmetric) {
255:       KSPSolve(tao->ksp, lclP->GU,  lclP->lamda); 
256:     } else {
257:       KSPSolveTranspose(tao->ksp, lclP->GU,  lclP->lamda); 
258:     }
259:     KSPGetIterationNumber(tao->ksp,&its); 
260:     tao->ksp_its += its;
261:   }
262:     
263:   VecCopy(lclP->lamda,lclP->lamda0); 

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

267:   LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V); 
268:   LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V); 
269:   

271:   /* Evaluate constraint norm */
272:   VecNorm(tao->constraints, NORM_2, &cnorm);  
273:   VecNorm(lclP->GAugL, NORM_2, &mnorm); 
274:   

276:   /* Monitor convergence */
277:   TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason); 

279:   while (reason == TAO_CONTINUE_ITERATING) {
280:     /* Compute a descent direction for the linearly constrained subproblem
281:        minimize f(u+du, v+dv)
282:        s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */

284:     /* Store the points around the linearization */
285:     VecCopy(lclP->U, lclP->U0); 
286:     VecCopy(lclP->V, lclP->V0); 
287:     VecCopy(lclP->GU,lclP->GU0); 
288:     VecCopy(lclP->GV,lclP->GV0); 
289:     VecCopy(lclP->GAugL_U,lclP->GAugL_U0); 
290:     VecCopy(lclP->GAugL_V,lclP->GAugL_V0); 
291:     VecCopy(lclP->GL_U,lclP->GL_U0); 
292:     VecCopy(lclP->GL_V,lclP->GL_V0); 

294:     lclP->aug0 = lclP->aug;
295:     lclP->lgn0 = lclP->lgn;

297:     /* Given the design variables, we need to project the current iterate
298:        onto the linearized constraint.  We choose to fix the design variables
299:        and solve the linear system for the state variables.  The resulting
300:        point is the Newton direction */
301:     
302:     /* Solve r = A\con */
303:     lclP->solve_type = LCL_FORWARD1;
304:     VecSet(lclP->r,0.0);  /*  Initial guess in CG */
305:     
306:     if (tao->jacobian_state_inv) {
307:       MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r); 
308:     } else {
309:       KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag); 
310:       KSPSolve(tao->ksp, tao->constraints,  lclP->r); 
311:       KSPGetIterationNumber(tao->ksp,&its); 
312:       tao->ksp_its += its;
313:     }
314:       

316:     VecNorm(lclP->r,NORM_2,&cnorm); 
317:     VecSet(lclP->s, 0.0); 

319:     /* Make sure the Newton direction is a descent direction for the merit function */
320:     if (symmetric)  {
321:       MatMult(tao->jacobian_state,tao->constraints,lclP->WU); 
322:     } else {
323:       MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU); 
324:     }
325:     
326:     VecDot(lclP->r,lclP->WU,&descent); 
327:     if (descent <= 0) {
328:       PetscInfo1(tao,"Newton direction not descent: %G\n",descent); 
329:       if (lclP->verbose) {
330:         PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent: %G\n",descent); 
331:       }
332:       PetscInfo(tao,"Using steepest descent direction instead.\n"); 
333:       VecSet(lclP->r,0.0); 
334:       VecAXPY(lclP->r,-1.0,lclP->WU); 
335:     }

337:     /* Check descent for aug. lagrangian */
338:     VecDot(lclP->r,lclP->r,&r2); 
339:     VecDot(lclP->r,lclP->GAugL_U,&descent);  
340:     adec = 1e-8 * r2;
341:     if (descent <= adec) {
342:       PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %G",descent); 
343:       VecDot(tao->constraints,tao->constraints,&con2); 
344:       VecDot(lclP->r,lclP->GL_U,&rGL_U);
345:       VecDot(lclP->r,lclP->WU,&rWU);
346:       lclP->rho = 2 * (adec - rGL_U) / rWU;
347:       PetscInfo1(tao,"  Increasing penalty parameter to %G",lclP->rho);
348:       if (lclP->verbose) {
349:         PetscPrintf(PETSC_COMM_WORLD,"  Increasing penalty parameter to %G\n",lclP->rho); 
350:       }
351:     }

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

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

364:     lclP->recompute_jacobian_flag = PETSC_TRUE;

366:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 
367:     TaoLineSearchSetType(tao->linesearch, TAOLINESEARCH_MT); 
368:     TaoLineSearchSetFromOptions(tao->linesearch); 
369:     TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason); 
370:     if (lclP->verbose) {
371:       PetscPrintf(PETSC_COMM_WORLD,"Steplength = %G\n",step); 
372:     }
373:     
374:     LCLScatter(lclP,tao->solution,lclP->U,lclP->V); 
375:     TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); 
376:     LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); 

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

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


389:     /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
390:     for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++){
391:       /* We now minimize the objective function starting from the fraction of
392:          the Newton point accepted by applying one step of a reduced-space
393:          method.  The optimization problem is:
394:        
395:          minimize f(u+du, v+dv)
396:          s. t.    A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)

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

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

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

408:       /* Compute multipliers */
409:       /* p1 */
410:       VecSet(lclP->lamda,0.0);  /*  Initial guess in CG */
411:       lclP->solve_type = LCL_ADJOINT1;
412:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); 
413:       if (tao->jacobian_state_pre) {
414:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); 
415:       } else {
416:         pset = pflag = PETSC_TRUE;
417:       }
418:       if (set && pset && flag && pflag)
419:         symmetric = PETSC_TRUE;
420:       else
421:         symmetric = PETSC_FALSE;
422:       
423:       if (tao->jacobian_state_inv) {
424:         if (symmetric) {
425:           MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda); 
426:         } else {
427:           MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda); 
428:         }
429:       } else {
430:         if (symmetric) {
431:           KSPSolve(tao->ksp, lclP->GAugL_U,  lclP->lamda); 
432:         } else {
433:           KSPSolveTranspose(tao->ksp, lclP->GAugL_U,  lclP->lamda); 
434:         }
435:         KSPGetIterationNumber(tao->ksp,&its); 
436:         tao->ksp_its += its;
437:       }
438:       MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1); 
439:       VecAXPY(lclP->g1,-1.0,lclP->GAugL_V); 
440:     
441:       /* Compute the limited-memory quasi-newton direction */
442:       if (iter > 0) {
443:         MatLMVMSolve(lclP->R,lclP->g1,lclP->s); 
444:         VecDot(lclP->s,lclP->g1,&descent); 
445:         if (descent < 0e-8) {
446:           if (lclP->verbose) {
447:             PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %G\n",descent); 
448:           }
449:           VecCopy(lclP->g1,lclP->s); 
450:         }
451:       } else {
452:         VecCopy(lclP->g1,lclP->s); 
453:       }
454:       VecScale(lclP->g1,-1.0); 


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

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

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

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

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

491:       TaoComputeJacobianState(tao,lclP->X0,&tao->jacobian_state,&tao->jacobian_state_pre,&tao->jacobian_state_inv,&lclP->statematflag); 
492:       TaoComputeJacobianDesign(tao,lclP->X0,&tao->jacobian_design); 

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

503:       MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); 
504:       if (tao->jacobian_state_pre) {
505:         MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); 
506:       } else {
507:         pset = pflag = PETSC_TRUE;
508:       }
509:       if (set && pset && flag && pflag)
510:         symmetric = PETSC_TRUE;
511:       else
512:         symmetric = PETSC_FALSE;

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

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

537:       /* Update the quasi-newton approximation */
538:       if (phase2_iter >= 0){
539:         MatLMVMSetPrev(lclP->R,lclP->V1,lclP->g1);
540:       }
541:       MatLMVMUpdate(lclP->R,lclP->V,lclP->g2); 
542:       /* 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 */

544:     }

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

553:     TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag); 
554:     TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design); 
555:     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); 
564:   
565:     /* Monitor convergence */
566:     iter++;
567:     TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason); 
568:     
569:   }

571:    MatDestroy(&lclP->R); 

573:    return(0);
574: }

579: PetscErrorCode TaoCreate_LCL(TaoSolver tao)
580: {
581:   TAO_LCL *lclP;
583:   const char *morethuente_type = TAOLINESEARCH_MT;
585:   
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:   

593:   PetscNewLog(tao,TAO_LCL,&lclP); 
594:   tao->data = (void*)lclP;


597:   tao->max_it=200;
598:   tao->fatol=1e-8;
599:   tao->frtol=1e-8;
600:   tao->catol=1e-4;
601:   tao->crtol=1e-4;
602:   tao->gttol=1e-4;
603:   tao->gatol=1e-4;
604:   tao->grtol=1e-4;
605:   
606:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
607:   TaoLineSearchSetType(tao->linesearch, morethuente_type); 

609:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao); 
610:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp); 
611:   KSPSetFromOptions(tao->ksp); 


614:   
615:   return(0);

617: }


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

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

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

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

668:   f[0] = lclP->lgn;
669:   VecCopy(lclP->GL,G);

671:   return(0);
672: }

676: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
677: {
678:   TaoSolver tao = (TaoSolver)ptr;
679:   TAO_LCL *lclP = (TAO_LCL*)tao->data;
680:   PetscReal con2;
681:   PetscBool flag,pflag,set,pset,symmetric;

685:   LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao); 
686:   LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V); 
687:   VecDot(tao->constraints,tao->constraints,&con2); 
688:   lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;
689:   
690:   /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
691:   /*      WU = A' * c
692:           WV = B' * c */
693:   MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); 
694:   if (tao->jacobian_state_pre) {
695:     MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); 
696:   } else {
697:     pset = pflag = PETSC_TRUE;
698:   }
699:   if (set && pset && flag && pflag)
700:     symmetric = PETSC_TRUE;
701:   else
702:     symmetric = PETSC_FALSE;

704:   if (symmetric) {
705:     MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);  
706:   } else {
707:     MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U); 
708:   }

710:   MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V); 
711:   VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U); 
712:   VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V); 
713:   LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL); 

715:   f[0] = lclP->aug;
716:   VecCopy(lclP->GAugL,G); 

718:   return(0);
719: }

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);
732:   
733: }
736: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
737: {
740:   VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD); 
741:   VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD); 
742:   VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD); 
743:   VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD); 
744:   return(0);
745:   
746: }