Actual source code: lcl.c
1: #include <../src/tao/pde_constrained/impls/lcl/lcl.h>
2: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
3: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
4: static PetscErrorCode LCLScatter(TAO_LCL*,Vec,Vec,Vec);
5: static PetscErrorCode LCLGather(TAO_LCL*,Vec,Vec,Vec);
7: static PetscErrorCode TaoDestroy_LCL(Tao tao)
8: {
9: TAO_LCL *lclP = (TAO_LCL*)tao->data;
11: if (tao->setupcalled) {
12: MatDestroy(&lclP->R);
13: VecDestroy(&lclP->lamda);
14: VecDestroy(&lclP->lamda0);
15: VecDestroy(&lclP->WL);
16: VecDestroy(&lclP->W);
17: VecDestroy(&lclP->X0);
18: VecDestroy(&lclP->G0);
19: VecDestroy(&lclP->GL);
20: VecDestroy(&lclP->GAugL);
21: VecDestroy(&lclP->dbar);
22: VecDestroy(&lclP->U);
23: VecDestroy(&lclP->U0);
24: VecDestroy(&lclP->V);
25: VecDestroy(&lclP->V0);
26: VecDestroy(&lclP->V1);
27: VecDestroy(&lclP->GU);
28: VecDestroy(&lclP->GV);
29: VecDestroy(&lclP->GU0);
30: VecDestroy(&lclP->GV0);
31: VecDestroy(&lclP->GL_U);
32: VecDestroy(&lclP->GL_V);
33: VecDestroy(&lclP->GAugL_U);
34: VecDestroy(&lclP->GAugL_V);
35: VecDestroy(&lclP->GL_U0);
36: VecDestroy(&lclP->GL_V0);
37: VecDestroy(&lclP->GAugL_U0);
38: VecDestroy(&lclP->GAugL_V0);
39: VecDestroy(&lclP->DU);
40: VecDestroy(&lclP->DV);
41: VecDestroy(&lclP->WU);
42: VecDestroy(&lclP->WV);
43: VecDestroy(&lclP->g1);
44: VecDestroy(&lclP->g2);
45: VecDestroy(&lclP->con1);
47: VecDestroy(&lclP->r);
48: VecDestroy(&lclP->s);
50: ISDestroy(&tao->state_is);
51: ISDestroy(&tao->design_is);
53: VecScatterDestroy(&lclP->state_scatter);
54: VecScatterDestroy(&lclP->design_scatter);
55: }
56: MatDestroy(&lclP->R);
57: PetscFree(tao->data);
58: return 0;
59: }
61: static PetscErrorCode TaoSetFromOptions_LCL(PetscOptionItems *PetscOptionsObject,Tao tao)
62: {
63: TAO_LCL *lclP = (TAO_LCL*)tao->data;
65: PetscOptionsHead(PetscOptionsObject,"Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
66: PetscOptionsReal("-tao_lcl_eps1","epsilon 1 tolerance","",lclP->eps1,&lclP->eps1,NULL);
67: PetscOptionsReal("-tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
68: PetscOptionsReal("-tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
69: PetscOptionsReal("-tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
70: lclP->phase2_niter = 1;
71: PetscOptionsInt("-tao_lcl_phase2_niter","Number of phase 2 iterations in LCL algorithm","",lclP->phase2_niter,&lclP->phase2_niter,NULL);
72: lclP->verbose = PETSC_FALSE;
73: PetscOptionsBool("-tao_lcl_verbose","Print verbose output","",lclP->verbose,&lclP->verbose,NULL);
74: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
75: PetscOptionsReal("-tao_lcl_tola","Tolerance for first forward solve","",lclP->tau[0],&lclP->tau[0],NULL);
76: PetscOptionsReal("-tao_lcl_tolb","Tolerance for first adjoint solve","",lclP->tau[1],&lclP->tau[1],NULL);
77: PetscOptionsReal("-tao_lcl_tolc","Tolerance for second forward solve","",lclP->tau[2],&lclP->tau[2],NULL);
78: PetscOptionsReal("-tao_lcl_told","Tolerance for second adjoint solve","",lclP->tau[3],&lclP->tau[3],NULL);
79: PetscOptionsTail();
80: TaoLineSearchSetFromOptions(tao->linesearch);
81: MatSetFromOptions(lclP->R);
82: return 0;
83: }
85: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
86: {
87: return 0;
88: }
90: static PetscErrorCode TaoSetup_LCL(Tao tao)
91: {
92: TAO_LCL *lclP = (TAO_LCL*)tao->data;
93: PetscInt lo, hi, nlocalstate, nlocaldesign;
94: IS is_state, is_design;
97: VecDuplicate(tao->solution, &tao->gradient);
98: VecDuplicate(tao->solution, &tao->stepdirection);
99: VecDuplicate(tao->solution, &lclP->W);
100: VecDuplicate(tao->solution, &lclP->X0);
101: VecDuplicate(tao->solution, &lclP->G0);
102: VecDuplicate(tao->solution, &lclP->GL);
103: VecDuplicate(tao->solution, &lclP->GAugL);
105: VecDuplicate(tao->constraints, &lclP->lamda);
106: VecDuplicate(tao->constraints, &lclP->WL);
107: VecDuplicate(tao->constraints, &lclP->lamda0);
108: VecDuplicate(tao->constraints, &lclP->con1);
110: VecSet(lclP->lamda,0.0);
112: VecGetSize(tao->solution, &lclP->n);
113: VecGetSize(tao->constraints, &lclP->m);
115: VecCreate(((PetscObject)tao)->comm,&lclP->U);
116: VecCreate(((PetscObject)tao)->comm,&lclP->V);
117: ISGetLocalSize(tao->state_is,&nlocalstate);
118: ISGetLocalSize(tao->design_is,&nlocaldesign);
119: VecSetSizes(lclP->U,nlocalstate,lclP->m);
120: VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m);
121: VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name);
122: VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name);
123: VecSetFromOptions(lclP->U);
124: VecSetFromOptions(lclP->V);
125: VecDuplicate(lclP->U,&lclP->DU);
126: VecDuplicate(lclP->U,&lclP->U0);
127: VecDuplicate(lclP->U,&lclP->GU);
128: VecDuplicate(lclP->U,&lclP->GU0);
129: VecDuplicate(lclP->U,&lclP->GAugL_U);
130: VecDuplicate(lclP->U,&lclP->GL_U);
131: VecDuplicate(lclP->U,&lclP->GAugL_U0);
132: VecDuplicate(lclP->U,&lclP->GL_U0);
133: VecDuplicate(lclP->U,&lclP->WU);
134: VecDuplicate(lclP->U,&lclP->r);
135: VecDuplicate(lclP->V,&lclP->V0);
136: VecDuplicate(lclP->V,&lclP->V1);
137: VecDuplicate(lclP->V,&lclP->DV);
138: VecDuplicate(lclP->V,&lclP->s);
139: VecDuplicate(lclP->V,&lclP->GV);
140: VecDuplicate(lclP->V,&lclP->GV0);
141: VecDuplicate(lclP->V,&lclP->dbar);
142: VecDuplicate(lclP->V,&lclP->GAugL_V);
143: VecDuplicate(lclP->V,&lclP->GL_V);
144: VecDuplicate(lclP->V,&lclP->GAugL_V0);
145: VecDuplicate(lclP->V,&lclP->GL_V0);
146: VecDuplicate(lclP->V,&lclP->WV);
147: VecDuplicate(lclP->V,&lclP->g1);
148: VecDuplicate(lclP->V,&lclP->g2);
150: /* create scatters for state, design subvecs */
151: VecGetOwnershipRange(lclP->U,&lo,&hi);
152: ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state);
153: VecGetOwnershipRange(lclP->V,&lo,&hi);
154: if (0) {
155: PetscInt sizeU,sizeV;
156: VecGetSize(lclP->U,&sizeU);
157: VecGetSize(lclP->V,&sizeV);
158: PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV);
159: }
160: ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design);
161: VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter);
162: VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter);
163: ISDestroy(&is_state);
164: ISDestroy(&is_design);
165: return 0;
166: }
168: static PetscErrorCode TaoSolve_LCL(Tao tao)
169: {
170: TAO_LCL *lclP = (TAO_LCL*)tao->data;
171: PetscInt phase2_iter,nlocal,its;
172: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
173: PetscReal step=1.0,f, descent, aldescent;
174: PetscReal cnorm, mnorm;
175: PetscReal adec,r2,rGL_U,rWU;
176: PetscBool set,pset,flag,pflag,symmetric;
178: lclP->rho = lclP->rho0;
179: VecGetLocalSize(lclP->U,&nlocal);
180: VecGetLocalSize(lclP->V,&nlocal);
181: MatSetSizes(lclP->R, nlocal, nlocal, lclP->n-lclP->m, lclP->n-lclP->m);
182: MatLMVMAllocate(lclP->R,lclP->V,lclP->V);
183: lclP->recompute_jacobian_flag = PETSC_TRUE;
185: /* Scatter to U,V */
186: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
188: /* Evaluate Function, Gradient, Constraints, and Jacobian */
189: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
190: TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
191: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
192: TaoComputeConstraints(tao,tao->solution, tao->constraints);
194: /* Scatter gradient to GU,GV */
195: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
197: /* Evaluate Lagrangian function and gradient */
198: /* p0 */
199: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
200: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
201: if (tao->jacobian_state_pre) {
202: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
203: } else {
204: pset = pflag = PETSC_TRUE;
205: }
206: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
207: else symmetric = PETSC_FALSE;
209: lclP->solve_type = LCL_ADJOINT2;
210: if (tao->jacobian_state_inv) {
211: if (symmetric) {
212: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
213: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
214: }
215: } else {
216: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
217: if (symmetric) {
218: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
219: } else {
220: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
221: }
222: KSPGetIterationNumber(tao->ksp,&its);
223: tao->ksp_its+=its;
224: tao->ksp_tot_its+=its;
225: }
226: VecCopy(lclP->lamda,lclP->lamda0);
227: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
229: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
230: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
232: /* Evaluate constraint norm */
233: VecNorm(tao->constraints, NORM_2, &cnorm);
234: VecNorm(lclP->GAugL, NORM_2, &mnorm);
236: /* Monitor convergence */
237: tao->reason = TAO_CONTINUE_ITERATING;
238: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
239: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
240: (*tao->ops->convergencetest)(tao,tao->cnvP);
242: while (tao->reason == TAO_CONTINUE_ITERATING) {
243: /* Call general purpose update function */
244: if (tao->ops->update) {
245: (*tao->ops->update)(tao, tao->niter, tao->user_update);
246: }
247: tao->ksp_its=0;
248: /* Compute a descent direction for the linearly constrained subproblem
249: minimize f(u+du, v+dv)
250: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
252: /* Store the points around the linearization */
253: VecCopy(lclP->U, lclP->U0);
254: VecCopy(lclP->V, lclP->V0);
255: VecCopy(lclP->GU,lclP->GU0);
256: VecCopy(lclP->GV,lclP->GV0);
257: VecCopy(lclP->GAugL_U,lclP->GAugL_U0);
258: VecCopy(lclP->GAugL_V,lclP->GAugL_V0);
259: VecCopy(lclP->GL_U,lclP->GL_U0);
260: VecCopy(lclP->GL_V,lclP->GL_V0);
262: lclP->aug0 = lclP->aug;
263: lclP->lgn0 = lclP->lgn;
265: /* Given the design variables, we need to project the current iterate
266: onto the linearized constraint. We choose to fix the design variables
267: and solve the linear system for the state variables. The resulting
268: point is the Newton direction */
270: /* Solve r = A\con */
271: lclP->solve_type = LCL_FORWARD1;
272: VecSet(lclP->r,0.0); /* Initial guess in CG */
274: if (tao->jacobian_state_inv) {
275: MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
276: } else {
277: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
278: KSPSolve(tao->ksp, tao->constraints, lclP->r);
279: KSPGetIterationNumber(tao->ksp,&its);
280: tao->ksp_its+=its;
281: tao->ksp_tot_its+=tao->ksp_its;
282: }
284: /* Set design step direction dv to zero */
285: VecSet(lclP->s, 0.0);
287: /*
288: Check sufficient descent for constraint merit function .5*||con||^2
289: con' Ak r >= eps1 ||r||^(2+eps2)
290: */
292: /* Compute WU= Ak' * con */
293: if (symmetric) {
294: MatMult(tao->jacobian_state,tao->constraints,lclP->WU);
295: } else {
296: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU);
297: }
298: /* Compute r * Ak' * con */
299: VecDot(lclP->r,lclP->WU,&rWU);
301: /* compute ||r||^(2+eps2) */
302: VecNorm(lclP->r,NORM_2,&r2);
303: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
304: adec = lclP->eps1 * r2;
306: if (rWU < adec) {
307: PetscInfo(tao,"Newton direction not descent for constraint, feasibility phase required\n");
308: if (lclP->verbose) {
309: PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %g -- using steepest descent\n",(double)descent);
310: }
312: PetscInfo(tao,"Using steepest descent direction instead.\n");
313: VecSet(lclP->r,0.0);
314: VecAXPY(lclP->r,-1.0,lclP->WU);
315: VecDot(lclP->r,lclP->r,&rWU);
316: VecNorm(lclP->r,NORM_2,&r2);
317: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
318: VecDot(lclP->r,lclP->GAugL_U,&descent);
319: adec = lclP->eps1 * r2;
320: }
322: /*
323: Check descent for aug. lagrangian
324: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
325: GL_U = GUk - Ak'*yk
326: WU = Ak'*con
327: adec=eps1||r||^(2+eps2)
329: ==>
330: Check r'GL_U - rho*r'WU <= adec
331: */
333: VecDot(lclP->r,lclP->GL_U,&rGL_U);
334: aldescent = rGL_U - lclP->rho*rWU;
335: if (aldescent > -adec) {
336: if (lclP->verbose) {
337: PetscPrintf(PETSC_COMM_WORLD," Newton direction not descent for augmented Lagrangian: %g",(double)aldescent);
338: }
339: PetscInfo(tao,"Newton direction not descent for augmented Lagrangian: %g\n",(double)aldescent);
340: lclP->rho = (rGL_U - adec)/rWU;
341: if (lclP->rho > lclP->rhomax) {
342: lclP->rho = lclP->rhomax;
343: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)",(double)lclP->rho);
344: }
345: if (lclP->verbose) {
346: PetscPrintf(PETSC_COMM_WORLD," Increasing penalty parameter to %g\n",(double)lclP->rho);
347: }
348: PetscInfo(tao," Increasing penalty parameter to %g\n",(double)lclP->rho);
349: }
351: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
352: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
353: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
355: /* We now minimize the augmented Lagrangian along the Newton direction */
356: VecScale(lclP->r,-1.0);
357: LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
358: VecScale(lclP->r,-1.0);
359: LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
360: LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);
362: lclP->recompute_jacobian_flag = PETSC_TRUE;
364: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
365: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
366: TaoLineSearchSetFromOptions(tao->linesearch);
367: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
368: if (lclP->verbose) {
369: PetscPrintf(PETSC_COMM_WORLD,"Steplength = %g\n",(double)step);
370: }
372: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
373: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
374: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
376: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
378: /* Check convergence */
379: VecNorm(lclP->GAugL, NORM_2, &mnorm);
380: VecNorm(tao->constraints, NORM_2, &cnorm);
381: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
382: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
383: (*tao->ops->convergencetest)(tao,tao->cnvP);
384: if (tao->reason != TAO_CONTINUE_ITERATING) {
385: break;
386: }
388: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
389: for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++) {
390: /* We now minimize the objective function starting from the fraction of
391: the Newton point accepted by applying one step of a reduced-space
392: method. The optimization problem is:
394: minimize f(u+du, v+dv)
395: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
397: In particular, we have that
398: du = -inv(A)*(Bdv + alpha g) */
400: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
401: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
403: /* Store V and constraints */
404: VecCopy(lclP->V, lclP->V1);
405: VecCopy(tao->constraints,lclP->con1);
407: /* Compute multipliers */
408: /* p1 */
409: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
410: lclP->solve_type = LCL_ADJOINT1;
411: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
412: if (tao->jacobian_state_pre) {
413: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
414: } else {
415: pset = pflag = PETSC_TRUE;
416: }
417: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
418: else symmetric = PETSC_FALSE;
420: if (tao->jacobian_state_inv) {
421: if (symmetric) {
422: MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
423: } else {
424: MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
425: }
426: } else {
427: if (symmetric) {
428: KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lamda);
429: } else {
430: KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lamda);
431: }
432: KSPGetIterationNumber(tao->ksp,&its);
433: tao->ksp_its+=its;
434: tao->ksp_tot_its+=its;
435: }
436: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
437: VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);
439: /* Compute the limited-memory quasi-newton direction */
440: if (tao->niter > 0) {
441: MatSolve(lclP->R,lclP->g1,lclP->s);
442: VecDot(lclP->s,lclP->g1,&descent);
443: if (descent <= 0) {
444: if (lclP->verbose) {
445: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %g\n",(double)descent);
446: }
447: VecCopy(lclP->g1,lclP->s);
448: }
449: } else {
450: VecCopy(lclP->g1,lclP->s);
451: }
452: VecScale(lclP->g1,-1.0);
454: /* Recover the full space direction */
455: MatMult(tao->jacobian_design,lclP->s,lclP->WU);
456: /* VecSet(lclP->r,0.0); */ /* Initial guess in CG */
457: lclP->solve_type = LCL_FORWARD2;
458: if (tao->jacobian_state_inv) {
459: MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
460: } else {
461: KSPSolve(tao->ksp, lclP->WU, lclP->r);
462: KSPGetIterationNumber(tao->ksp,&its);
463: tao->ksp_its += its;
464: tao->ksp_tot_its+=its;
465: }
467: /* We now minimize the augmented Lagrangian along the direction -r,s */
468: VecScale(lclP->r, -1.0);
469: LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
470: VecScale(lclP->r, -1.0);
471: lclP->recompute_jacobian_flag = PETSC_TRUE;
473: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
474: TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT);
475: TaoLineSearchSetFromOptions(tao->linesearch);
476: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
477: if (lclP->verbose) {
478: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength = %g\n",(double)step);
479: }
481: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
482: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
483: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
484: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
485: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
487: /* Compute the reduced gradient at the new point */
489: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
490: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
492: /* p2 */
493: /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
494: if (phase2_iter==0) {
495: VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
496: } else {
497: VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
498: }
500: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
501: if (tao->jacobian_state_pre) {
502: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
503: } else {
504: pset = pflag = PETSC_TRUE;
505: }
506: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
507: else symmetric = PETSC_FALSE;
509: lclP->solve_type = LCL_ADJOINT2;
510: if (tao->jacobian_state_inv) {
511: if (symmetric) {
512: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
513: } else {
514: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
515: }
516: } else {
517: if (symmetric) {
518: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
519: } else {
520: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
521: }
522: KSPGetIterationNumber(tao->ksp,&its);
523: tao->ksp_its += its;
524: tao->ksp_tot_its += its;
525: }
527: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g2);
528: VecAXPY(lclP->g2,-1.0,lclP->GV);
530: VecScale(lclP->g2,-1.0);
532: /* Update the quasi-newton approximation */
533: MatLMVMUpdate(lclP->R,lclP->V,lclP->g2);
534: /* Use "-tao_ls_type gpcg -tao_ls_ftol 0 -tao_lmm_broyden_phi 0.0 -tao_lmm_scale_type scalar" to obtain agreement with Matlab code */
536: }
538: VecCopy(lclP->lamda,lclP->lamda0);
540: /* Evaluate Function, Gradient, Constraints, and Jacobian */
541: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
542: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
543: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
545: TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
546: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
547: TaoComputeConstraints(tao,tao->solution, tao->constraints);
549: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
551: VecNorm(lclP->GAugL, NORM_2, &mnorm);
553: /* Evaluate constraint norm */
554: VecNorm(tao->constraints, NORM_2, &cnorm);
556: /* Monitor convergence */
557: tao->niter++;
558: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
559: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
560: (*tao->ops->convergencetest)(tao,tao->cnvP);
561: }
562: return 0;
563: }
565: /*MC
566: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
568: + -tao_lcl_eps1 - epsilon 1 tolerance
569: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
570: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
571: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
572: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
573: . -tao_lcl_verbose - Print verbose output if True
574: . -tao_lcl_tola - Tolerance for first forward solve
575: . -tao_lcl_tolb - Tolerance for first adjoint solve
576: . -tao_lcl_tolc - Tolerance for second forward solve
577: - -tao_lcl_told - Tolerance for second adjoint solve
579: Level: beginner
580: M*/
581: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
582: {
583: TAO_LCL *lclP;
584: const char *morethuente_type = TAOLINESEARCHMT;
586: tao->ops->setup = TaoSetup_LCL;
587: tao->ops->solve = TaoSolve_LCL;
588: tao->ops->view = TaoView_LCL;
589: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
590: tao->ops->destroy = TaoDestroy_LCL;
591: PetscNewLog(tao,&lclP);
592: tao->data = (void*)lclP;
594: /* Override default settings (unless already changed) */
595: if (!tao->max_it_changed) tao->max_it = 200;
596: if (!tao->catol_changed) tao->catol = 1.0e-4;
597: if (!tao->gatol_changed) tao->gttol = 1.0e-4;
598: if (!tao->grtol_changed) tao->gttol = 1.0e-4;
599: if (!tao->gttol_changed) tao->gttol = 1.0e-4;
600: lclP->rho0 = 1.0e-4;
601: lclP->rhomax=1e5;
602: lclP->eps1 = 1.0e-8;
603: lclP->eps2 = 0.0;
604: lclP->solve_type=2;
605: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
606: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
607: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
608: TaoLineSearchSetType(tao->linesearch, morethuente_type);
609: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
611: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
612: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
613: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
614: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
615: KSPSetFromOptions(tao->ksp);
617: MatCreate(((PetscObject)tao)->comm, &lclP->R);
618: MatSetType(lclP->R, MATLMVMBFGS);
619: return 0;
620: }
622: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
623: {
624: Tao tao = (Tao)ptr;
625: TAO_LCL *lclP = (TAO_LCL*)tao->data;
626: PetscBool set,pset,flag,pflag,symmetric;
627: PetscReal cdotl;
629: TaoComputeObjectiveAndGradient(tao,X,f,G);
630: LCLScatter(lclP,G,lclP->GU,lclP->GV);
631: if (lclP->recompute_jacobian_flag) {
632: TaoComputeJacobianState(tao,X,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
633: TaoComputeJacobianDesign(tao,X,tao->jacobian_design);
634: }
635: TaoComputeConstraints(tao,X, tao->constraints);
636: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
637: if (tao->jacobian_state_pre) {
638: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
639: } else {
640: pset = pflag = PETSC_TRUE;
641: }
642: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
643: else symmetric = PETSC_FALSE;
645: VecDot(lclP->lamda0, tao->constraints, &cdotl);
646: lclP->lgn = *f - cdotl;
648: /* Gradient of Lagrangian GL = G - J' * lamda */
649: /* WU = A' * WL
650: WV = B' * WL */
651: if (symmetric) {
652: MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
653: } else {
654: MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
655: }
656: MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
657: VecScale(lclP->GL_U,-1.0);
658: VecScale(lclP->GL_V,-1.0);
659: VecAXPY(lclP->GL_U,1.0,lclP->GU);
660: VecAXPY(lclP->GL_V,1.0,lclP->GV);
661: LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);
663: f[0] = lclP->lgn;
664: VecCopy(lclP->GL,G);
665: return 0;
666: }
668: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
669: {
670: Tao tao = (Tao)ptr;
671: TAO_LCL *lclP = (TAO_LCL*)tao->data;
672: PetscReal con2;
673: PetscBool flag,pflag,set,pset,symmetric;
675: LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
676: LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
677: VecDot(tao->constraints,tao->constraints,&con2);
678: lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;
680: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
681: /* WU = A' * c
682: WV = B' * c */
683: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
684: if (tao->jacobian_state_pre) {
685: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
686: } else {
687: pset = pflag = PETSC_TRUE;
688: }
689: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
690: else symmetric = PETSC_FALSE;
692: if (symmetric) {
693: MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
694: } else {
695: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
696: }
698: MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
699: VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
700: VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
701: LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);
703: f[0] = lclP->aug;
704: VecCopy(lclP->GAugL,G);
705: return 0;
706: }
708: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
709: {
710: VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
711: VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
712: VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
713: VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
714: return 0;
716: }
717: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
718: {
719: VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
720: VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
721: VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
722: VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
723: return 0;
725: }