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