Actual source code: lcl.c
petsc-3.5.4 2015-05-23
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(Tao tao)
68: {
69: TAO_LCL *lclP = (TAO_LCL*)tao->data;
70: PetscBool flg;
74: PetscOptionsHead("Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
75: PetscOptionsReal("-tao_lcl_eps1","epsilon 1 tolerance","",lclP->eps1,&lclP->eps1,&flg);
76: PetscOptionsReal("-tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,&flg);
77: PetscOptionsReal("-tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,&flg);
78: PetscOptionsReal("-tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,&flg);
79: lclP->phase2_niter = 1;
80: PetscOptionsInt("-tao_lcl_phase2_niter","Number of phase 2 iterations in LCL algorithm","",lclP->phase2_niter,&lclP->phase2_niter,&flg);
81: lclP->verbose = PETSC_FALSE;
82: PetscOptionsBool("-tao_lcl_verbose","Print verbose output","",lclP->verbose,&lclP->verbose,&flg);
83: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
84: PetscOptionsReal("-tao_lcl_tola","Tolerance for first forward solve","",lclP->tau[0],&lclP->tau[0],&flg);
85: PetscOptionsReal("-tao_lcl_tolb","Tolerance for first adjoint solve","",lclP->tau[1],&lclP->tau[1],&flg);
86: PetscOptionsReal("-tao_lcl_tolc","Tolerance for second forward solve","",lclP->tau[2],&lclP->tau[2],&flg);
87: PetscOptionsReal("-tao_lcl_told","Tolerance for second adjoint solve","",lclP->tau[3],&lclP->tau[3],&flg);
88: PetscOptionsTail();
89: TaoLineSearchSetFromOptions(tao->linesearch);
90: return(0);
91: }
95: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
96: {
97: return 0;
98: }
102: static PetscErrorCode TaoSetup_LCL(Tao tao)
103: {
104: TAO_LCL *lclP = (TAO_LCL*)tao->data;
105: PetscInt lo, hi, nlocalstate, nlocaldesign;
107: IS is_state, is_design;
110: if (!tao->state_is) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
111: VecDuplicate(tao->solution, &tao->gradient);
112: VecDuplicate(tao->solution, &tao->stepdirection);
113: VecDuplicate(tao->solution, &lclP->W);
114: VecDuplicate(tao->solution, &lclP->X0);
115: VecDuplicate(tao->solution, &lclP->G0);
116: VecDuplicate(tao->solution, &lclP->GL);
117: VecDuplicate(tao->solution, &lclP->GAugL);
119: VecDuplicate(tao->constraints, &lclP->lamda);
120: VecDuplicate(tao->constraints, &lclP->WL);
121: VecDuplicate(tao->constraints, &lclP->lamda0);
122: VecDuplicate(tao->constraints, &lclP->con1);
124: VecSet(lclP->lamda,0.0);
126: VecGetSize(tao->solution, &lclP->n);
127: VecGetSize(tao->constraints, &lclP->m);
129: VecCreate(((PetscObject)tao)->comm,&lclP->U);
130: VecCreate(((PetscObject)tao)->comm,&lclP->V);
131: ISGetLocalSize(tao->state_is,&nlocalstate);
132: ISGetLocalSize(tao->design_is,&nlocaldesign);
133: VecSetSizes(lclP->U,nlocalstate,lclP->m);
134: VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m);
135: VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name);
136: VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name);
137: VecSetFromOptions(lclP->U);
138: VecSetFromOptions(lclP->V);
139: VecDuplicate(lclP->U,&lclP->DU);
140: VecDuplicate(lclP->U,&lclP->U0);
141: VecDuplicate(lclP->U,&lclP->GU);
142: VecDuplicate(lclP->U,&lclP->GU0);
143: VecDuplicate(lclP->U,&lclP->GAugL_U);
144: VecDuplicate(lclP->U,&lclP->GL_U);
145: VecDuplicate(lclP->U,&lclP->GAugL_U0);
146: VecDuplicate(lclP->U,&lclP->GL_U0);
147: VecDuplicate(lclP->U,&lclP->WU);
148: VecDuplicate(lclP->U,&lclP->r);
149: VecDuplicate(lclP->V,&lclP->V0);
150: VecDuplicate(lclP->V,&lclP->V1);
151: VecDuplicate(lclP->V,&lclP->DV);
152: VecDuplicate(lclP->V,&lclP->s);
153: VecDuplicate(lclP->V,&lclP->GV);
154: VecDuplicate(lclP->V,&lclP->GV0);
155: VecDuplicate(lclP->V,&lclP->dbar);
156: VecDuplicate(lclP->V,&lclP->GAugL_V);
157: VecDuplicate(lclP->V,&lclP->GL_V);
158: VecDuplicate(lclP->V,&lclP->GAugL_V0);
159: VecDuplicate(lclP->V,&lclP->GL_V0);
160: VecDuplicate(lclP->V,&lclP->WV);
161: VecDuplicate(lclP->V,&lclP->g1);
162: VecDuplicate(lclP->V,&lclP->g2);
164: /* create scatters for state, design subvecs */
165: VecGetOwnershipRange(lclP->U,&lo,&hi);
166: ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state);
167: VecGetOwnershipRange(lclP->V,&lo,&hi);
168: if (0) {
169: PetscInt sizeU,sizeV;
170: VecGetSize(lclP->U,&sizeU);
171: VecGetSize(lclP->V,&sizeV);
172: PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV);
173: }
174: ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design);
175: VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter);
176: VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter);
177: ISDestroy(&is_state);
178: ISDestroy(&is_design);
179: return(0);
180: }
184: static PetscErrorCode TaoSolve_LCL(Tao tao)
185: {
186: TAO_LCL *lclP = (TAO_LCL*)tao->data;
187: PetscInt iter=0,phase2_iter,nlocal,its;
188: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
189: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
190: PetscReal step=1.0,f, descent, aldescent;
191: PetscReal cnorm, mnorm;
192: PetscReal adec,r2,rGL_U,rWU;
193: PetscBool set,pset,flag,pflag,symmetric;
194: PetscErrorCode ierr;
197: lclP->rho = lclP->rho0;
198: VecGetLocalSize(lclP->U,&nlocal);
199: VecGetLocalSize(lclP->V,&nlocal);
200: MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R);
201: MatLMVMAllocateVectors(lclP->R,lclP->V);
202: lclP->recompute_jacobian_flag = PETSC_TRUE;
204: /* Scatter to U,V */
205: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
207: /* Evaluate Function, Gradient, Constraints, and Jacobian */
208: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
209: TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
210: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
211: TaoComputeConstraints(tao,tao->solution, tao->constraints);
213: /* Scatter gradient to GU,GV */
214: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
216: /* Evaluate Lagrangian function and gradient */
217: /* p0 */
218: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
219: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
220: if (tao->jacobian_state_pre) {
221: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
222: } else {
223: pset = pflag = PETSC_TRUE;
224: }
225: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
226: else symmetric = PETSC_FALSE;
228: lclP->solve_type = LCL_ADJOINT2;
229: if (tao->jacobian_state_inv) {
230: if (symmetric) {
231: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
232: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
233: }
234: } else {
235: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
236: if (symmetric) {
237: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
238: } else {
239: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
240: }
241: KSPGetIterationNumber(tao->ksp,&its);
242: tao->ksp_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, iter,f,mnorm,cnorm,step,&reason);
257: while (reason == TAO_CONTINUE_ITERATING) {
258: /* Compute a descent direction for the linearly constrained subproblem
259: minimize f(u+du, v+dv)
260: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
262: /* Store the points around the linearization */
263: VecCopy(lclP->U, lclP->U0);
264: VecCopy(lclP->V, lclP->V0);
265: VecCopy(lclP->GU,lclP->GU0);
266: VecCopy(lclP->GV,lclP->GV0);
267: VecCopy(lclP->GAugL_U,lclP->GAugL_U0);
268: VecCopy(lclP->GAugL_V,lclP->GAugL_V0);
269: VecCopy(lclP->GL_U,lclP->GL_U0);
270: VecCopy(lclP->GL_V,lclP->GL_V0);
272: lclP->aug0 = lclP->aug;
273: lclP->lgn0 = lclP->lgn;
275: /* Given the design variables, we need to project the current iterate
276: onto the linearized constraint. We choose to fix the design variables
277: and solve the linear system for the state variables. The resulting
278: point is the Newton direction */
280: /* Solve r = A\con */
281: lclP->solve_type = LCL_FORWARD1;
282: VecSet(lclP->r,0.0); /* Initial guess in CG */
284: if (tao->jacobian_state_inv) {
285: MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
286: } else {
287: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
288: KSPSolve(tao->ksp, tao->constraints, lclP->r);
289: KSPGetIterationNumber(tao->ksp,&its);
290: tao->ksp_its += its;
291: }
293: /* Set design step direction dv to zero */
294: VecSet(lclP->s, 0.0);
296: /*
297: Check sufficient descent for constraint merit function .5*||con||^2
298: con' Ak r >= eps1 ||r||^(2+eps2)
299: */
301: /* Compute WU= Ak' * con */
302: if (symmetric) {
303: MatMult(tao->jacobian_state,tao->constraints,lclP->WU);
304: } else {
305: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU);
306: }
307: /* Compute r * Ak' * con */
308: VecDot(lclP->r,lclP->WU,&rWU);
310: /* compute ||r||^(2+eps2) */
311: VecNorm(lclP->r,NORM_2,&r2);
312: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
313: adec = lclP->eps1 * r2;
315: if (rWU < adec) {
316: PetscInfo(tao,"Newton direction not descent for constraint, feasibility phase required");
317: if (lclP->verbose) {
318: PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %g -- using steepest descent\n",(double)descent);
319: }
321: PetscInfo(tao,"Using steepest descent direction instead.\n");
322: VecSet(lclP->r,0.0);
323: VecAXPY(lclP->r,-1.0,lclP->WU);
324: VecDot(lclP->r,lclP->r,&rWU);
325: VecNorm(lclP->r,NORM_2,&r2);
326: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
327: VecDot(lclP->r,lclP->GAugL_U,&descent);
328: adec = lclP->eps1 * r2;
329: }
332: /*
333: Check descent for aug. lagrangian
334: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
335: GL_U = GUk - Ak'*yk
336: WU = Ak'*con
337: adec=eps1||r||^(2+eps2)
339: ==>
340: Check r'GL_U - rho*r'WU <= adec
341: */
343: VecDot(lclP->r,lclP->GL_U,&rGL_U);
344: aldescent = rGL_U - lclP->rho*rWU;
345: if (aldescent > -adec) {
346: if (lclP->verbose) {
347: PetscPrintf(PETSC_COMM_WORLD," Newton direction not descent for augmented Lagrangian: %g",(double)aldescent);
348: }
349: PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %g",(double)aldescent);
350: lclP->rho = (rGL_U - adec)/rWU;
351: if (lclP->rho > lclP->rhomax) {
352: lclP->rho = lclP->rhomax;
353: SETERRQ1(PETSC_COMM_WORLD,0,"rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)",(double)lclP->rho);
354: }
355: if (lclP->verbose) {
356: PetscPrintf(PETSC_COMM_WORLD," Increasing penalty parameter to %g\n",(double)lclP->rho);
357: }
358: PetscInfo1(tao," Increasing penalty parameter to %g",(double)lclP->rho);
359: }
361: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
362: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
363: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
365: /* We now minimize the augmented Lagrangian along the Newton direction */
366: VecScale(lclP->r,-1.0);
367: LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
368: VecScale(lclP->r,-1.0);
369: LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
370: LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);
372: lclP->recompute_jacobian_flag = PETSC_TRUE;
374: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
375: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
376: TaoLineSearchSetFromOptions(tao->linesearch);
377: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
378: if (lclP->verbose) {
379: PetscPrintf(PETSC_COMM_WORLD,"Steplength = %g\n",(double)step);
380: }
382: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
383: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
384: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
386: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
388: /* Check convergence */
389: VecNorm(lclP->GAugL, NORM_2, &mnorm);
390: VecNorm(tao->constraints, NORM_2, &cnorm);
391: TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);
392: if (reason != TAO_CONTINUE_ITERATING){
393: break;
394: }
396: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
397: for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++){
398: /* We now minimize the objective function starting from the fraction of
399: the Newton point accepted by applying one step of a reduced-space
400: method. The optimization problem is:
402: minimize f(u+du, v+dv)
403: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
405: In particular, we have that
406: du = -inv(A)*(Bdv + alpha g) */
408: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
409: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
411: /* Store V and constraints */
412: VecCopy(lclP->V, lclP->V1);
413: VecCopy(tao->constraints,lclP->con1);
415: /* Compute multipliers */
416: /* p1 */
417: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
418: lclP->solve_type = LCL_ADJOINT1;
419: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
420: if (tao->jacobian_state_pre) {
421: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
422: } else {
423: pset = pflag = PETSC_TRUE;
424: }
425: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
426: else symmetric = PETSC_FALSE;
428: if (tao->jacobian_state_inv) {
429: if (symmetric) {
430: MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
431: } else {
432: MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
433: }
434: } else {
435: if (symmetric) {
436: KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lamda);
437: } else {
438: KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lamda);
439: }
440: KSPGetIterationNumber(tao->ksp,&its);
441: tao->ksp_its += its;
442: }
443: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
444: VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);
446: /* Compute the limited-memory quasi-newton direction */
447: if (iter > 0) {
448: MatLMVMSolve(lclP->R,lclP->g1,lclP->s);
449: VecDot(lclP->s,lclP->g1,&descent);
450: if (descent <= 0) {
451: if (lclP->verbose) {
452: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %g\n",(double)descent);
453: }
454: VecCopy(lclP->g1,lclP->s);
455: }
456: } else {
457: VecCopy(lclP->g1,lclP->s);
458: }
459: VecScale(lclP->g1,-1.0);
461: /* Recover the full space direction */
462: MatMult(tao->jacobian_design,lclP->s,lclP->WU);
463: /* VecSet(lclP->r,0.0); */ /* Initial guess in CG */
464: lclP->solve_type = LCL_FORWARD2;
465: if (tao->jacobian_state_inv) {
466: MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
467: } else {
468: KSPSolve(tao->ksp, lclP->WU, lclP->r);
469: KSPGetIterationNumber(tao->ksp,&its);
470: tao->ksp_its += its;
471: }
473: /* We now minimize the augmented Lagrangian along the direction -r,s */
474: VecScale(lclP->r, -1.0);
475: LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
476: VecScale(lclP->r, -1.0);
477: lclP->recompute_jacobian_flag = PETSC_TRUE;
479: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
480: TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT);
481: TaoLineSearchSetFromOptions(tao->linesearch);
482: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
483: if (lclP->verbose){
484: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength = %g\n",(double)step);
485: }
487: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
488: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
489: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
490: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
491: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
493: /* Compute the reduced gradient at the new point */
495: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
496: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
498: /* p2 */
499: /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
500: if (phase2_iter==0){
501: VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
502: } else {
503: VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
504: }
506: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
507: if (tao->jacobian_state_pre) {
508: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
509: } else {
510: pset = pflag = PETSC_TRUE;
511: }
512: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
513: else symmetric = PETSC_FALSE;
515: lclP->solve_type = LCL_ADJOINT2;
516: if (tao->jacobian_state_inv) {
517: if (symmetric) {
518: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
519: } else {
520: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
521: }
522: } else {
523: if (symmetric) {
524: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
525: } else {
526: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
527: }
528: KSPGetIterationNumber(tao->ksp,&its);
529: tao->ksp_its += its;
530: }
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);
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);
554: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
555: TaoComputeConstraints(tao,tao->solution, tao->constraints);
557: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
559: VecNorm(lclP->GAugL, NORM_2, &mnorm);
561: /* Evaluate constraint norm */
562: VecNorm(tao->constraints, NORM_2, &cnorm);
564: /* Monitor convergence */
565: iter++;
566: TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);
567: }
568: MatDestroy(&lclP->R);
569: return(0);
570: }
572: /*MC
573: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
575: + -tao_lcl_eps1 - epsilon 1 tolerance
576: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,&flg);
577: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,&flg);
578: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,&flg);
579: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
580: . -tao_lcl_verbose - Print verbose output if True
581: . -tao_lcl_tola - Tolerance for first forward solve
582: . -tao_lcl_tolb - Tolerance for first adjoint solve
583: . -tao_lcl_tolc - Tolerance for second forward solve
584: - -tao_lcl_told - Tolerance for second adjoint solve
586: Level: beginner
587: M*/
588: EXTERN_C_BEGIN
591: PetscErrorCode TaoCreate_LCL(Tao tao)
592: {
593: TAO_LCL *lclP;
595: const char *morethuente_type = TAOLINESEARCHMT;
598: tao->ops->setup = TaoSetup_LCL;
599: tao->ops->solve = TaoSolve_LCL;
600: tao->ops->view = TaoView_LCL;
601: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
602: tao->ops->destroy = TaoDestroy_LCL;
603: PetscNewLog(tao,&lclP);
604: tao->data = (void*)lclP;
606: tao->max_it=200;
607: #if defined(PETSC_USE_REAL_SINGLE)
608: tao->fatol=1e-5;
609: tao->frtol=1e-5;
610: #else
611: tao->fatol=1e-8;
612: tao->frtol=1e-8;
613: #endif
614: tao->catol=1e-4;
615: tao->crtol=1e-4;
616: tao->gttol=1e-4;
617: tao->gatol=1e-4;
618: tao->grtol=1e-4;
619: lclP->rho0 = 1.0e-4;
620: lclP->rhomax=1e5;
621: lclP->eps1 = 1.0e-8;
622: lclP->eps2 = 0.0;
623: lclP->solve_type=2;
624: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
625: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
626: TaoLineSearchSetType(tao->linesearch, morethuente_type);
628: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
629: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
630: KSPSetFromOptions(tao->ksp);
631: return(0);
632: }
633: EXTERN_C_END
637: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
638: {
639: Tao tao = (Tao)ptr;
640: TAO_LCL *lclP = (TAO_LCL*)tao->data;
641: PetscBool set,pset,flag,pflag,symmetric;
642: PetscReal cdotl;
646: TaoComputeObjectiveAndGradient(tao,X,f,G);
647: LCLScatter(lclP,G,lclP->GU,lclP->GV);
648: if (lclP->recompute_jacobian_flag) {
649: TaoComputeJacobianState(tao,X,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
650: TaoComputeJacobianDesign(tao,X,tao->jacobian_design);
651: }
652: TaoComputeConstraints(tao,X, tao->constraints);
653: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
654: if (tao->jacobian_state_pre) {
655: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
656: } else {
657: pset = pflag = PETSC_TRUE;
658: }
659: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
660: else symmetric = PETSC_FALSE;
662: VecDot(lclP->lamda0, tao->constraints, &cdotl);
663: lclP->lgn = *f - cdotl;
665: /* Gradient of Lagrangian GL = G - J' * lamda */
666: /* WU = A' * WL
667: WV = B' * WL */
668: if (symmetric) {
669: MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
670: } else {
671: MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
672: }
673: MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
674: VecScale(lclP->GL_U,-1.0);
675: VecScale(lclP->GL_V,-1.0);
676: VecAXPY(lclP->GL_U,1.0,lclP->GU);
677: VecAXPY(lclP->GL_V,1.0,lclP->GV);
678: LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);
680: f[0] = lclP->lgn;
681: VecCopy(lclP->GL,G);
682: return(0);
683: }
687: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
688: {
689: Tao tao = (Tao)ptr;
690: TAO_LCL *lclP = (TAO_LCL*)tao->data;
691: PetscReal con2;
692: PetscBool flag,pflag,set,pset,symmetric;
696: LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
697: LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
698: VecDot(tao->constraints,tao->constraints,&con2);
699: lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;
701: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
702: /* WU = A' * c
703: WV = B' * c */
704: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
705: if (tao->jacobian_state_pre) {
706: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
707: } else {
708: pset = pflag = PETSC_TRUE;
709: }
710: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
711: else symmetric = PETSC_FALSE;
713: if (symmetric) {
714: MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
715: } else {
716: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
717: }
719: MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
720: VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
721: VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
722: LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);
724: f[0] = lclP->aug;
725: VecCopy(lclP->GAugL,G);
726: return(0);
727: }
731: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
732: {
735: VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
736: VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
737: VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
738: VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
739: return(0);
741: }
744: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
745: {
748: VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
749: VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
750: VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
751: VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
752: return(0);
754: }