Actual source code: lcl.c
petsc-3.7.5 2017-01-01
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: }