Actual source code: lcl.c
petsc-3.9.4 2018-09-11
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);
8: static PetscErrorCode TaoDestroy_LCL(Tao tao)
9: {
10: TAO_LCL *lclP = (TAO_LCL*)tao->data;
14: if (tao->setupcalled) {
15: MatDestroy(&lclP->R);
16: VecDestroy(&lclP->lamda);
17: VecDestroy(&lclP->lamda0);
18: VecDestroy(&lclP->WL);
19: VecDestroy(&lclP->W);
20: VecDestroy(&lclP->X0);
21: VecDestroy(&lclP->G0);
22: VecDestroy(&lclP->GL);
23: VecDestroy(&lclP->GAugL);
24: VecDestroy(&lclP->dbar);
25: VecDestroy(&lclP->U);
26: VecDestroy(&lclP->U0);
27: VecDestroy(&lclP->V);
28: VecDestroy(&lclP->V0);
29: VecDestroy(&lclP->V1);
30: VecDestroy(&lclP->GU);
31: VecDestroy(&lclP->GV);
32: VecDestroy(&lclP->GU0);
33: VecDestroy(&lclP->GV0);
34: VecDestroy(&lclP->GL_U);
35: VecDestroy(&lclP->GL_V);
36: VecDestroy(&lclP->GAugL_U);
37: VecDestroy(&lclP->GAugL_V);
38: VecDestroy(&lclP->GL_U0);
39: VecDestroy(&lclP->GL_V0);
40: VecDestroy(&lclP->GAugL_U0);
41: VecDestroy(&lclP->GAugL_V0);
42: VecDestroy(&lclP->DU);
43: VecDestroy(&lclP->DV);
44: VecDestroy(&lclP->WU);
45: VecDestroy(&lclP->WV);
46: VecDestroy(&lclP->g1);
47: VecDestroy(&lclP->g2);
48: VecDestroy(&lclP->con1);
50: VecDestroy(&lclP->r);
51: VecDestroy(&lclP->s);
53: ISDestroy(&tao->state_is);
54: ISDestroy(&tao->design_is);
56: VecScatterDestroy(&lclP->state_scatter);
57: VecScatterDestroy(&lclP->design_scatter);
58: }
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: return(0);
86: }
88: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
89: {
90: return 0;
91: }
93: static PetscErrorCode TaoSetup_LCL(Tao tao)
94: {
95: TAO_LCL *lclP = (TAO_LCL*)tao->data;
96: PetscInt lo, hi, nlocalstate, nlocaldesign;
98: IS is_state, is_design;
101: if (!tao->state_is) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
102: VecDuplicate(tao->solution, &tao->gradient);
103: VecDuplicate(tao->solution, &tao->stepdirection);
104: VecDuplicate(tao->solution, &lclP->W);
105: VecDuplicate(tao->solution, &lclP->X0);
106: VecDuplicate(tao->solution, &lclP->G0);
107: VecDuplicate(tao->solution, &lclP->GL);
108: VecDuplicate(tao->solution, &lclP->GAugL);
110: VecDuplicate(tao->constraints, &lclP->lamda);
111: VecDuplicate(tao->constraints, &lclP->WL);
112: VecDuplicate(tao->constraints, &lclP->lamda0);
113: VecDuplicate(tao->constraints, &lclP->con1);
115: VecSet(lclP->lamda,0.0);
117: VecGetSize(tao->solution, &lclP->n);
118: VecGetSize(tao->constraints, &lclP->m);
120: VecCreate(((PetscObject)tao)->comm,&lclP->U);
121: VecCreate(((PetscObject)tao)->comm,&lclP->V);
122: ISGetLocalSize(tao->state_is,&nlocalstate);
123: ISGetLocalSize(tao->design_is,&nlocaldesign);
124: VecSetSizes(lclP->U,nlocalstate,lclP->m);
125: VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m);
126: VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name);
127: VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name);
128: VecSetFromOptions(lclP->U);
129: VecSetFromOptions(lclP->V);
130: VecDuplicate(lclP->U,&lclP->DU);
131: VecDuplicate(lclP->U,&lclP->U0);
132: VecDuplicate(lclP->U,&lclP->GU);
133: VecDuplicate(lclP->U,&lclP->GU0);
134: VecDuplicate(lclP->U,&lclP->GAugL_U);
135: VecDuplicate(lclP->U,&lclP->GL_U);
136: VecDuplicate(lclP->U,&lclP->GAugL_U0);
137: VecDuplicate(lclP->U,&lclP->GL_U0);
138: VecDuplicate(lclP->U,&lclP->WU);
139: VecDuplicate(lclP->U,&lclP->r);
140: VecDuplicate(lclP->V,&lclP->V0);
141: VecDuplicate(lclP->V,&lclP->V1);
142: VecDuplicate(lclP->V,&lclP->DV);
143: VecDuplicate(lclP->V,&lclP->s);
144: VecDuplicate(lclP->V,&lclP->GV);
145: VecDuplicate(lclP->V,&lclP->GV0);
146: VecDuplicate(lclP->V,&lclP->dbar);
147: VecDuplicate(lclP->V,&lclP->GAugL_V);
148: VecDuplicate(lclP->V,&lclP->GL_V);
149: VecDuplicate(lclP->V,&lclP->GAugL_V0);
150: VecDuplicate(lclP->V,&lclP->GL_V0);
151: VecDuplicate(lclP->V,&lclP->WV);
152: VecDuplicate(lclP->V,&lclP->g1);
153: VecDuplicate(lclP->V,&lclP->g2);
155: /* create scatters for state, design subvecs */
156: VecGetOwnershipRange(lclP->U,&lo,&hi);
157: ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state);
158: VecGetOwnershipRange(lclP->V,&lo,&hi);
159: if (0) {
160: PetscInt sizeU,sizeV;
161: VecGetSize(lclP->U,&sizeU);
162: VecGetSize(lclP->V,&sizeV);
163: PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV);
164: }
165: ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design);
166: VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter);
167: VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter);
168: ISDestroy(&is_state);
169: ISDestroy(&is_design);
170: return(0);
171: }
173: static PetscErrorCode TaoSolve_LCL(Tao tao)
174: {
175: TAO_LCL *lclP = (TAO_LCL*)tao->data;
176: PetscInt phase2_iter,nlocal,its;
177: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
178: PetscReal step=1.0,f, descent, aldescent;
179: PetscReal cnorm, mnorm;
180: PetscReal adec,r2,rGL_U,rWU;
181: PetscBool set,pset,flag,pflag,symmetric;
182: PetscErrorCode ierr;
185: lclP->rho = lclP->rho0;
186: VecGetLocalSize(lclP->U,&nlocal);
187: VecGetLocalSize(lclP->V,&nlocal);
188: MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R);
189: MatLMVMAllocateVectors(lclP->R,lclP->V);
190: lclP->recompute_jacobian_flag = PETSC_TRUE;
192: /* Scatter to U,V */
193: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
195: /* Evaluate Function, Gradient, Constraints, and Jacobian */
196: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
197: TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
198: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
199: TaoComputeConstraints(tao,tao->solution, tao->constraints);
201: /* Scatter gradient to GU,GV */
202: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
204: /* Evaluate Lagrangian function and gradient */
205: /* p0 */
206: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
207: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
208: if (tao->jacobian_state_pre) {
209: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
210: } else {
211: pset = pflag = PETSC_TRUE;
212: }
213: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
214: else symmetric = PETSC_FALSE;
216: lclP->solve_type = LCL_ADJOINT2;
217: if (tao->jacobian_state_inv) {
218: if (symmetric) {
219: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
220: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
221: }
222: } else {
223: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
224: if (symmetric) {
225: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
226: } else {
227: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
228: }
229: KSPGetIterationNumber(tao->ksp,&its);
230: tao->ksp_its+=its;
231: tao->ksp_tot_its+=its;
232: }
233: VecCopy(lclP->lamda,lclP->lamda0);
234: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
236: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
237: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
239: /* Evaluate constraint norm */
240: VecNorm(tao->constraints, NORM_2, &cnorm);
241: VecNorm(lclP->GAugL, NORM_2, &mnorm);
243: /* Monitor convergence */
244: tao->reason = TAO_CONTINUE_ITERATING;
245: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
246: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
247: (*tao->ops->convergencetest)(tao,tao->cnvP);
249: while (tao->reason == TAO_CONTINUE_ITERATING) {
250: tao->ksp_its=0;
251: /* Compute a descent direction for the linearly constrained subproblem
252: minimize f(u+du, v+dv)
253: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
255: /* Store the points around the linearization */
256: VecCopy(lclP->U, lclP->U0);
257: VecCopy(lclP->V, lclP->V0);
258: VecCopy(lclP->GU,lclP->GU0);
259: VecCopy(lclP->GV,lclP->GV0);
260: VecCopy(lclP->GAugL_U,lclP->GAugL_U0);
261: VecCopy(lclP->GAugL_V,lclP->GAugL_V0);
262: VecCopy(lclP->GL_U,lclP->GL_U0);
263: VecCopy(lclP->GL_V,lclP->GL_V0);
265: lclP->aug0 = lclP->aug;
266: lclP->lgn0 = lclP->lgn;
268: /* Given the design variables, we need to project the current iterate
269: onto the linearized constraint. We choose to fix the design variables
270: and solve the linear system for the state variables. The resulting
271: point is the Newton direction */
273: /* Solve r = A\con */
274: lclP->solve_type = LCL_FORWARD1;
275: VecSet(lclP->r,0.0); /* Initial guess in CG */
277: if (tao->jacobian_state_inv) {
278: MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
279: } else {
280: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
281: KSPSolve(tao->ksp, tao->constraints, lclP->r);
282: KSPGetIterationNumber(tao->ksp,&its);
283: tao->ksp_its+=its;
284: tao->ksp_tot_its+=tao->ksp_its;
285: }
287: /* Set design step direction dv to zero */
288: VecSet(lclP->s, 0.0);
290: /*
291: Check sufficient descent for constraint merit function .5*||con||^2
292: con' Ak r >= eps1 ||r||^(2+eps2)
293: */
295: /* Compute WU= Ak' * con */
296: if (symmetric) {
297: MatMult(tao->jacobian_state,tao->constraints,lclP->WU);
298: } else {
299: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU);
300: }
301: /* Compute r * Ak' * con */
302: VecDot(lclP->r,lclP->WU,&rWU);
304: /* compute ||r||^(2+eps2) */
305: VecNorm(lclP->r,NORM_2,&r2);
306: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
307: adec = lclP->eps1 * r2;
309: if (rWU < adec) {
310: PetscInfo(tao,"Newton direction not descent for constraint, feasibility phase required\n");
311: if (lclP->verbose) {
312: PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %g -- using steepest descent\n",(double)descent);
313: }
315: PetscInfo(tao,"Using steepest descent direction instead.\n");
316: VecSet(lclP->r,0.0);
317: VecAXPY(lclP->r,-1.0,lclP->WU);
318: VecDot(lclP->r,lclP->r,&rWU);
319: VecNorm(lclP->r,NORM_2,&r2);
320: r2 = PetscPowScalar(r2,2.0+lclP->eps2);
321: VecDot(lclP->r,lclP->GAugL_U,&descent);
322: adec = lclP->eps1 * r2;
323: }
326: /*
327: Check descent for aug. lagrangian
328: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
329: GL_U = GUk - Ak'*yk
330: WU = Ak'*con
331: adec=eps1||r||^(2+eps2)
333: ==>
334: Check r'GL_U - rho*r'WU <= adec
335: */
337: VecDot(lclP->r,lclP->GL_U,&rGL_U);
338: aldescent = rGL_U - lclP->rho*rWU;
339: if (aldescent > -adec) {
340: if (lclP->verbose) {
341: PetscPrintf(PETSC_COMM_WORLD," Newton direction not descent for augmented Lagrangian: %g",(double)aldescent);
342: }
343: PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %g\n",(double)aldescent);
344: lclP->rho = (rGL_U - adec)/rWU;
345: if (lclP->rho > lclP->rhomax) {
346: lclP->rho = lclP->rhomax;
347: SETERRQ1(PETSC_COMM_WORLD,0,"rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)",(double)lclP->rho);
348: }
349: if (lclP->verbose) {
350: PetscPrintf(PETSC_COMM_WORLD," Increasing penalty parameter to %g\n",(double)lclP->rho);
351: }
352: PetscInfo1(tao," Increasing penalty parameter to %g\n",(double)lclP->rho);
353: }
355: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
356: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
357: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
359: /* We now minimize the augmented Lagrangian along the Newton direction */
360: VecScale(lclP->r,-1.0);
361: LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
362: VecScale(lclP->r,-1.0);
363: LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
364: LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);
366: lclP->recompute_jacobian_flag = PETSC_TRUE;
368: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
369: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
370: TaoLineSearchSetFromOptions(tao->linesearch);
371: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
372: if (lclP->verbose) {
373: PetscPrintf(PETSC_COMM_WORLD,"Steplength = %g\n",(double)step);
374: }
376: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
377: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
378: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
380: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
382: /* Check convergence */
383: VecNorm(lclP->GAugL, NORM_2, &mnorm);
384: VecNorm(tao->constraints, NORM_2, &cnorm);
385: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
386: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
387: (*tao->ops->convergencetest)(tao,tao->cnvP);
388: if (tao->reason != TAO_CONTINUE_ITERATING){
389: break;
390: }
392: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
393: for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++){
394: /* We now minimize the objective function starting from the fraction of
395: the Newton point accepted by applying one step of a reduced-space
396: method. The optimization problem is:
398: minimize f(u+du, v+dv)
399: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
401: In particular, we have that
402: du = -inv(A)*(Bdv + alpha g) */
404: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
405: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
407: /* Store V and constraints */
408: VecCopy(lclP->V, lclP->V1);
409: VecCopy(tao->constraints,lclP->con1);
411: /* Compute multipliers */
412: /* p1 */
413: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
414: lclP->solve_type = LCL_ADJOINT1;
415: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
416: if (tao->jacobian_state_pre) {
417: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
418: } else {
419: pset = pflag = PETSC_TRUE;
420: }
421: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
422: else symmetric = PETSC_FALSE;
424: if (tao->jacobian_state_inv) {
425: if (symmetric) {
426: MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
427: } else {
428: MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
429: }
430: } else {
431: if (symmetric) {
432: KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lamda);
433: } else {
434: KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lamda);
435: }
436: KSPGetIterationNumber(tao->ksp,&its);
437: tao->ksp_its+=its;
438: tao->ksp_tot_its+=its;
439: }
440: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
441: VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);
443: /* Compute the limited-memory quasi-newton direction */
444: if (tao->niter > 0) {
445: MatLMVMSolve(lclP->R,lclP->g1,lclP->s);
446: VecDot(lclP->s,lclP->g1,&descent);
447: if (descent <= 0) {
448: if (lclP->verbose) {
449: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %g\n",(double)descent);
450: }
451: VecCopy(lclP->g1,lclP->s);
452: }
453: } else {
454: VecCopy(lclP->g1,lclP->s);
455: }
456: VecScale(lclP->g1,-1.0);
458: /* Recover the full space direction */
459: MatMult(tao->jacobian_design,lclP->s,lclP->WU);
460: /* VecSet(lclP->r,0.0); */ /* Initial guess in CG */
461: lclP->solve_type = LCL_FORWARD2;
462: if (tao->jacobian_state_inv) {
463: MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
464: } else {
465: KSPSolve(tao->ksp, lclP->WU, lclP->r);
466: KSPGetIterationNumber(tao->ksp,&its);
467: tao->ksp_its += its;
468: tao->ksp_tot_its+=its;
469: }
471: /* We now minimize the augmented Lagrangian along the direction -r,s */
472: VecScale(lclP->r, -1.0);
473: LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
474: VecScale(lclP->r, -1.0);
475: lclP->recompute_jacobian_flag = PETSC_TRUE;
477: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
478: TaoLineSearchSetType(tao->linesearch,TAOLINESEARCHMT);
479: TaoLineSearchSetFromOptions(tao->linesearch);
480: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
481: if (lclP->verbose){
482: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength = %g\n",(double)step);
483: }
485: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
486: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
487: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
488: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
489: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
491: /* Compute the reduced gradient at the new point */
493: TaoComputeJacobianState(tao,lclP->X0,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
494: TaoComputeJacobianDesign(tao,lclP->X0,tao->jacobian_design);
496: /* p2 */
497: /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
498: if (phase2_iter==0){
499: VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
500: } else {
501: VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
502: }
504: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
505: if (tao->jacobian_state_pre) {
506: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
507: } else {
508: pset = pflag = PETSC_TRUE;
509: }
510: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
511: else symmetric = PETSC_FALSE;
513: lclP->solve_type = LCL_ADJOINT2;
514: if (tao->jacobian_state_inv) {
515: if (symmetric) {
516: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
517: } else {
518: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
519: }
520: } else {
521: if (symmetric) {
522: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
523: } else {
524: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
525: }
526: KSPGetIterationNumber(tao->ksp,&its);
527: tao->ksp_its += its;
528: tao->ksp_tot_its += its;
529: }
531: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g2);
532: VecAXPY(lclP->g2,-1.0,lclP->GV);
534: VecScale(lclP->g2,-1.0);
536: /* Update the quasi-newton approximation */
537: if (phase2_iter >= 0){
538: MatLMVMSetPrev(lclP->R,lclP->V1,lclP->g1);
539: }
540: MatLMVMUpdate(lclP->R,lclP->V,lclP->g2);
541: /* 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 */
543: }
545: VecCopy(lclP->lamda,lclP->lamda0);
547: /* Evaluate Function, Gradient, Constraints, and Jacobian */
548: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
549: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
550: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
552: TaoComputeJacobianState(tao,tao->solution,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
553: TaoComputeJacobianDesign(tao,tao->solution,tao->jacobian_design);
554: TaoComputeConstraints(tao,tao->solution, tao->constraints);
556: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
558: VecNorm(lclP->GAugL, NORM_2, &mnorm);
560: /* Evaluate constraint norm */
561: VecNorm(tao->constraints, NORM_2, &cnorm);
563: /* Monitor convergence */
564: tao->niter++;
565: TaoLogConvergenceHistory(tao,f,mnorm,cnorm,tao->ksp_its);
566: TaoMonitor(tao,tao->niter,f,mnorm,cnorm,step);
567: (*tao->ops->convergencetest)(tao,tao->cnvP);
568: }
569: MatDestroy(&lclP->R);
570: return(0);
571: }
573: /*MC
574: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
576: + -tao_lcl_eps1 - epsilon 1 tolerance
577: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
578: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
579: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
580: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
581: . -tao_lcl_verbose - Print verbose output if True
582: . -tao_lcl_tola - Tolerance for first forward solve
583: . -tao_lcl_tolb - Tolerance for first adjoint solve
584: . -tao_lcl_tolc - Tolerance for second forward solve
585: - -tao_lcl_told - Tolerance for second adjoint solve
587: Level: beginner
588: M*/
589: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
590: {
591: TAO_LCL *lclP;
593: const char *morethuente_type = TAOLINESEARCHMT;
596: tao->ops->setup = TaoSetup_LCL;
597: tao->ops->solve = TaoSolve_LCL;
598: tao->ops->view = TaoView_LCL;
599: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
600: tao->ops->destroy = TaoDestroy_LCL;
601: PetscNewLog(tao,&lclP);
602: tao->data = (void*)lclP;
604: /* Override default settings (unless already changed) */
605: if (!tao->max_it_changed) tao->max_it = 200;
606: if (!tao->catol_changed) tao->catol = 1.0e-4;
607: if (!tao->gatol_changed) tao->gttol = 1.0e-4;
608: if (!tao->grtol_changed) tao->gttol = 1.0e-4;
609: if (!tao->gttol_changed) tao->gttol = 1.0e-4;
610: lclP->rho0 = 1.0e-4;
611: lclP->rhomax=1e5;
612: lclP->eps1 = 1.0e-8;
613: lclP->eps2 = 0.0;
614: lclP->solve_type=2;
615: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
616: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
617: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
618: TaoLineSearchSetType(tao->linesearch, morethuente_type);
619: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
621: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
622: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
623: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
624: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
625: KSPSetFromOptions(tao->ksp);
626: return(0);
627: }
629: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
630: {
631: Tao tao = (Tao)ptr;
632: TAO_LCL *lclP = (TAO_LCL*)tao->data;
633: PetscBool set,pset,flag,pflag,symmetric;
634: PetscReal cdotl;
638: TaoComputeObjectiveAndGradient(tao,X,f,G);
639: LCLScatter(lclP,G,lclP->GU,lclP->GV);
640: if (lclP->recompute_jacobian_flag) {
641: TaoComputeJacobianState(tao,X,tao->jacobian_state,tao->jacobian_state_pre,tao->jacobian_state_inv);
642: TaoComputeJacobianDesign(tao,X,tao->jacobian_design);
643: }
644: TaoComputeConstraints(tao,X, tao->constraints);
645: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
646: if (tao->jacobian_state_pre) {
647: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
648: } else {
649: pset = pflag = PETSC_TRUE;
650: }
651: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
652: else symmetric = PETSC_FALSE;
654: VecDot(lclP->lamda0, tao->constraints, &cdotl);
655: lclP->lgn = *f - cdotl;
657: /* Gradient of Lagrangian GL = G - J' * lamda */
658: /* WU = A' * WL
659: WV = B' * WL */
660: if (symmetric) {
661: MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
662: } else {
663: MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
664: }
665: MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
666: VecScale(lclP->GL_U,-1.0);
667: VecScale(lclP->GL_V,-1.0);
668: VecAXPY(lclP->GL_U,1.0,lclP->GU);
669: VecAXPY(lclP->GL_V,1.0,lclP->GV);
670: LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);
672: f[0] = lclP->lgn;
673: VecCopy(lclP->GL,G);
674: return(0);
675: }
677: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
678: {
679: Tao tao = (Tao)ptr;
680: TAO_LCL *lclP = (TAO_LCL*)tao->data;
681: PetscReal con2;
682: PetscBool flag,pflag,set,pset,symmetric;
686: LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
687: LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
688: VecDot(tao->constraints,tao->constraints,&con2);
689: lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;
691: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
692: /* WU = A' * c
693: WV = B' * c */
694: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
695: if (tao->jacobian_state_pre) {
696: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
697: } else {
698: pset = pflag = PETSC_TRUE;
699: }
700: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
701: else symmetric = PETSC_FALSE;
703: if (symmetric) {
704: MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
705: } else {
706: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
707: }
709: MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
710: VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
711: VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
712: LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);
714: f[0] = lclP->aug;
715: VecCopy(lclP->GAugL,G);
716: return(0);
717: }
719: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
720: {
723: VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
724: VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
725: VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
726: VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
727: return(0);
729: }
730: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
731: {
734: VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
735: VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
736: VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
737: VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
738: return(0);
740: }