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->lambda);
14: VecDestroy(&lclP->lambda0);
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: KSPDestroy(&tao->ksp);
58: PetscFree(tao->data);
59: return 0;
60: }
62: static PetscErrorCode TaoSetFromOptions_LCL(Tao tao, PetscOptionItems *PetscOptionsObject)
63: {
64: TAO_LCL *lclP = (TAO_LCL *)tao->data;
66: PetscOptionsHeadBegin(PetscOptionsObject, "Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
67: PetscOptionsReal("-tao_lcl_eps1", "epsilon 1 tolerance", "", lclP->eps1, &lclP->eps1, NULL);
68: PetscOptionsReal("-tao_lcl_eps2", "epsilon 2 tolerance", "", lclP->eps2, &lclP->eps2, NULL);
69: PetscOptionsReal("-tao_lcl_rho0", "init value for rho", "", lclP->rho0, &lclP->rho0, NULL);
70: PetscOptionsReal("-tao_lcl_rhomax", "max value for rho", "", lclP->rhomax, &lclP->rhomax, NULL);
71: lclP->phase2_niter = 1;
72: PetscOptionsInt("-tao_lcl_phase2_niter", "Number of phase 2 iterations in LCL algorithm", "", lclP->phase2_niter, &lclP->phase2_niter, NULL);
73: lclP->verbose = PETSC_FALSE;
74: PetscOptionsBool("-tao_lcl_verbose", "Print verbose output", "", lclP->verbose, &lclP->verbose, NULL);
75: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
76: PetscOptionsReal("-tao_lcl_tola", "Tolerance for first forward solve", "", lclP->tau[0], &lclP->tau[0], NULL);
77: PetscOptionsReal("-tao_lcl_tolb", "Tolerance for first adjoint solve", "", lclP->tau[1], &lclP->tau[1], NULL);
78: PetscOptionsReal("-tao_lcl_tolc", "Tolerance for second forward solve", "", lclP->tau[2], &lclP->tau[2], NULL);
79: PetscOptionsReal("-tao_lcl_told", "Tolerance for second adjoint solve", "", lclP->tau[3], &lclP->tau[3], NULL);
80: PetscOptionsHeadEnd();
81: TaoLineSearchSetFromOptions(tao->linesearch);
82: MatSetFromOptions(lclP->R);
83: return 0;
84: }
86: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
87: {
88: return 0;
89: }
91: static PetscErrorCode TaoSetup_LCL(Tao tao)
92: {
93: TAO_LCL *lclP = (TAO_LCL *)tao->data;
94: PetscInt lo, hi, nlocalstate, nlocaldesign;
95: IS is_state, is_design;
98: VecDuplicate(tao->solution, &tao->gradient);
99: VecDuplicate(tao->solution, &tao->stepdirection);
100: VecDuplicate(tao->solution, &lclP->W);
101: VecDuplicate(tao->solution, &lclP->X0);
102: VecDuplicate(tao->solution, &lclP->G0);
103: VecDuplicate(tao->solution, &lclP->GL);
104: VecDuplicate(tao->solution, &lclP->GAugL);
106: VecDuplicate(tao->constraints, &lclP->lambda);
107: VecDuplicate(tao->constraints, &lclP->WL);
108: VecDuplicate(tao->constraints, &lclP->lambda0);
109: VecDuplicate(tao->constraints, &lclP->con1);
111: VecSet(lclP->lambda, 0.0);
113: VecGetSize(tao->solution, &lclP->n);
114: VecGetSize(tao->constraints, &lclP->m);
116: VecCreate(((PetscObject)tao)->comm, &lclP->U);
117: VecCreate(((PetscObject)tao)->comm, &lclP->V);
118: ISGetLocalSize(tao->state_is, &nlocalstate);
119: ISGetLocalSize(tao->design_is, &nlocaldesign);
120: VecSetSizes(lclP->U, nlocalstate, lclP->m);
121: VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m);
122: VecSetType(lclP->U, ((PetscObject)(tao->solution))->type_name);
123: VecSetType(lclP->V, ((PetscObject)(tao->solution))->type_name);
124: VecSetFromOptions(lclP->U);
125: VecSetFromOptions(lclP->V);
126: VecDuplicate(lclP->U, &lclP->DU);
127: VecDuplicate(lclP->U, &lclP->U0);
128: VecDuplicate(lclP->U, &lclP->GU);
129: VecDuplicate(lclP->U, &lclP->GU0);
130: VecDuplicate(lclP->U, &lclP->GAugL_U);
131: VecDuplicate(lclP->U, &lclP->GL_U);
132: VecDuplicate(lclP->U, &lclP->GAugL_U0);
133: VecDuplicate(lclP->U, &lclP->GL_U0);
134: VecDuplicate(lclP->U, &lclP->WU);
135: VecDuplicate(lclP->U, &lclP->r);
136: VecDuplicate(lclP->V, &lclP->V0);
137: VecDuplicate(lclP->V, &lclP->V1);
138: VecDuplicate(lclP->V, &lclP->DV);
139: VecDuplicate(lclP->V, &lclP->s);
140: VecDuplicate(lclP->V, &lclP->GV);
141: VecDuplicate(lclP->V, &lclP->GV0);
142: VecDuplicate(lclP->V, &lclP->dbar);
143: VecDuplicate(lclP->V, &lclP->GAugL_V);
144: VecDuplicate(lclP->V, &lclP->GL_V);
145: VecDuplicate(lclP->V, &lclP->GAugL_V0);
146: VecDuplicate(lclP->V, &lclP->GL_V0);
147: VecDuplicate(lclP->V, &lclP->WV);
148: VecDuplicate(lclP->V, &lclP->g1);
149: VecDuplicate(lclP->V, &lclP->g2);
151: /* create scatters for state, design subvecs */
152: VecGetOwnershipRange(lclP->U, &lo, &hi);
153: ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state);
154: VecGetOwnershipRange(lclP->V, &lo, &hi);
155: if (0) {
156: PetscInt sizeU, sizeV;
157: VecGetSize(lclP->U, &sizeU);
158: VecGetSize(lclP->V, &sizeV);
159: PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV);
160: }
161: ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design);
162: VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter);
163: VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter);
164: ISDestroy(&is_state);
165: ISDestroy(&is_design);
166: return 0;
167: }
169: static PetscErrorCode TaoSolve_LCL(Tao tao)
170: {
171: TAO_LCL *lclP = (TAO_LCL *)tao->data;
172: PetscInt phase2_iter, nlocal, its;
173: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
174: PetscReal step = 1.0, f, descent, aldescent;
175: PetscReal cnorm, mnorm;
176: PetscReal adec, r2, rGL_U, rWU;
177: PetscBool set, pset, flag, pflag, symmetric;
179: lclP->rho = lclP->rho0;
180: VecGetLocalSize(lclP->U, &nlocal);
181: VecGetLocalSize(lclP->V, &nlocal);
182: MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m);
183: MatLMVMAllocate(lclP->R, lclP->V, lclP->V);
184: lclP->recompute_jacobian_flag = PETSC_TRUE;
186: /* Scatter to U,V */
187: LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
189: /* Evaluate Function, Gradient, Constraints, and Jacobian */
190: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
191: TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
192: TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design);
193: TaoComputeConstraints(tao, tao->solution, tao->constraints);
195: /* Scatter gradient to GU,GV */
196: LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);
198: /* Evaluate Lagrangian function and gradient */
199: /* p0 */
200: VecSet(lclP->lambda, 0.0); /* Initial guess in CG */
201: MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
202: if (tao->jacobian_state_pre) {
203: MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
204: } else {
205: pset = pflag = PETSC_TRUE;
206: }
207: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
208: else symmetric = PETSC_FALSE;
210: lclP->solve_type = LCL_ADJOINT2;
211: if (tao->jacobian_state_inv) {
212: if (symmetric) {
213: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
214: } else {
215: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
216: }
217: } else {
218: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre);
219: if (symmetric) {
220: KSPSolve(tao->ksp, lclP->GU, lclP->lambda);
221: } else {
222: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda);
223: }
224: KSPGetIterationNumber(tao->ksp, &its);
225: tao->ksp_its += its;
226: tao->ksp_tot_its += its;
227: }
228: VecCopy(lclP->lambda, lclP->lambda0);
229: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);
231: LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
232: LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);
234: /* Evaluate constraint norm */
235: VecNorm(tao->constraints, NORM_2, &cnorm);
236: VecNorm(lclP->GAugL, NORM_2, &mnorm);
238: /* Monitor convergence */
239: tao->reason = TAO_CONTINUE_ITERATING;
240: TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
241: TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
242: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
244: while (tao->reason == TAO_CONTINUE_ITERATING) {
245: /* Call general purpose update function */
246: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
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) PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent);
310: PetscInfo(tao, "Using steepest descent direction instead.\n");
311: VecSet(lclP->r, 0.0);
312: VecAXPY(lclP->r, -1.0, lclP->WU);
313: VecDot(lclP->r, lclP->r, &rWU);
314: VecNorm(lclP->r, NORM_2, &r2);
315: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
316: VecDot(lclP->r, lclP->GAugL_U, &descent);
317: adec = lclP->eps1 * r2;
318: }
320: /*
321: Check descent for aug. lagrangian
322: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
323: GL_U = GUk - Ak'*yk
324: WU = Ak'*con
325: adec=eps1||r||^(2+eps2)
327: ==>
328: Check r'GL_U - rho*r'WU <= adec
329: */
331: VecDot(lclP->r, lclP->GL_U, &rGL_U);
332: aldescent = rGL_U - lclP->rho * rWU;
333: if (aldescent > -adec) {
334: if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent);
335: PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent);
336: lclP->rho = (rGL_U - adec) / rWU;
337: if (lclP->rho > lclP->rhomax) {
338: lclP->rho = lclP->rhomax;
339: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)", (double)lclP->rho);
340: }
341: if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, " Increasing penalty parameter to %g\n", (double)lclP->rho);
342: PetscInfo(tao, " Increasing penalty parameter to %g\n", (double)lclP->rho);
343: }
345: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);
346: LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
347: LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);
349: /* We now minimize the augmented Lagrangian along the Newton direction */
350: VecScale(lclP->r, -1.0);
351: LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection);
352: VecScale(lclP->r, -1.0);
353: LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
354: LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0);
356: lclP->recompute_jacobian_flag = PETSC_TRUE;
358: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
359: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
360: TaoLineSearchSetFromOptions(tao->linesearch);
361: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
362: if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step);
364: LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
365: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
366: LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);
368: LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);
370: /* Check convergence */
371: VecNorm(lclP->GAugL, NORM_2, &mnorm);
372: VecNorm(tao->constraints, NORM_2, &cnorm);
373: TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
374: TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
375: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
376: if (tao->reason != TAO_CONTINUE_ITERATING) break;
378: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
379: for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
380: /* We now minimize the objective function starting from the fraction of
381: the Newton point accepted by applying one step of a reduced-space
382: method. The optimization problem is:
384: minimize f(u+du, v+dv)
385: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
387: In particular, we have that
388: du = -inv(A)*(Bdv + alpha g) */
390: TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
391: TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design);
393: /* Store V and constraints */
394: VecCopy(lclP->V, lclP->V1);
395: VecCopy(tao->constraints, lclP->con1);
397: /* Compute multipliers */
398: /* p1 */
399: VecSet(lclP->lambda, 0.0); /* Initial guess in CG */
400: lclP->solve_type = LCL_ADJOINT1;
401: MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
402: if (tao->jacobian_state_pre) {
403: MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
404: } else {
405: pset = pflag = PETSC_TRUE;
406: }
407: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
408: else symmetric = PETSC_FALSE;
410: if (tao->jacobian_state_inv) {
411: if (symmetric) {
412: MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda);
413: } else {
414: MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda);
415: }
416: } else {
417: if (symmetric) {
418: KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda);
419: } else {
420: KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda);
421: }
422: KSPGetIterationNumber(tao->ksp, &its);
423: tao->ksp_its += its;
424: tao->ksp_tot_its += its;
425: }
426: MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1);
427: VecAXPY(lclP->g1, -1.0, lclP->GAugL_V);
429: /* Compute the limited-memory quasi-newton direction */
430: if (tao->niter > 0) {
431: MatSolve(lclP->R, lclP->g1, lclP->s);
432: VecDot(lclP->s, lclP->g1, &descent);
433: if (descent <= 0) {
434: if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent);
435: VecCopy(lclP->g1, lclP->s);
436: }
437: } else {
438: VecCopy(lclP->g1, lclP->s);
439: }
440: VecScale(lclP->g1, -1.0);
442: /* Recover the full space direction */
443: MatMult(tao->jacobian_design, lclP->s, lclP->WU);
444: /* VecSet(lclP->r,0.0); */ /* Initial guess in CG */
445: lclP->solve_type = LCL_FORWARD2;
446: if (tao->jacobian_state_inv) {
447: MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r);
448: } else {
449: KSPSolve(tao->ksp, lclP->WU, lclP->r);
450: KSPGetIterationNumber(tao->ksp, &its);
451: tao->ksp_its += its;
452: tao->ksp_tot_its += its;
453: }
455: /* We now minimize the augmented Lagrangian along the direction -r,s */
456: VecScale(lclP->r, -1.0);
457: LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection);
458: VecScale(lclP->r, -1.0);
459: lclP->recompute_jacobian_flag = PETSC_TRUE;
461: TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
462: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT);
463: TaoLineSearchSetFromOptions(tao->linesearch);
464: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
465: if (lclP->verbose) PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength = %g\n", (double)step);
467: LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
468: LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V);
469: LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V);
470: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
471: LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);
473: /* Compute the reduced gradient at the new point */
475: TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
476: TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design);
478: /* p2 */
479: /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
480: if (phase2_iter == 0) {
481: VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0);
482: } else {
483: VecAXPY(lclP->lambda, -lclP->rho, tao->constraints);
484: }
486: MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
487: if (tao->jacobian_state_pre) {
488: MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
489: } else {
490: pset = pflag = PETSC_TRUE;
491: }
492: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
493: else symmetric = PETSC_FALSE;
495: lclP->solve_type = LCL_ADJOINT2;
496: if (tao->jacobian_state_inv) {
497: if (symmetric) {
498: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
499: } else {
500: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda);
501: }
502: } else {
503: if (symmetric) {
504: KSPSolve(tao->ksp, lclP->GU, lclP->lambda);
505: } else {
506: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda);
507: }
508: KSPGetIterationNumber(tao->ksp, &its);
509: tao->ksp_its += its;
510: tao->ksp_tot_its += its;
511: }
513: MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2);
514: VecAXPY(lclP->g2, -1.0, lclP->GV);
516: VecScale(lclP->g2, -1.0);
518: /* Update the quasi-newton approximation */
519: MatLMVMUpdate(lclP->R, lclP->V, lclP->g2);
520: /* 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 */
521: }
523: VecCopy(lclP->lambda, lclP->lambda0);
525: /* Evaluate Function, Gradient, Constraints, and Jacobian */
526: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
527: LCLScatter(lclP, tao->solution, lclP->U, lclP->V);
528: LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV);
530: TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
531: TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design);
532: TaoComputeConstraints(tao, tao->solution, tao->constraints);
534: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao);
536: VecNorm(lclP->GAugL, NORM_2, &mnorm);
538: /* Evaluate constraint norm */
539: VecNorm(tao->constraints, NORM_2, &cnorm);
541: /* Monitor convergence */
542: tao->niter++;
543: TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its);
544: TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step);
545: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
546: }
547: return 0;
548: }
550: /*MC
551: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
553: + -tao_lcl_eps1 - epsilon 1 tolerance
554: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
555: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
556: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
557: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
558: . -tao_lcl_verbose - Print verbose output if True
559: . -tao_lcl_tola - Tolerance for first forward solve
560: . -tao_lcl_tolb - Tolerance for first adjoint solve
561: . -tao_lcl_tolc - Tolerance for second forward solve
562: - -tao_lcl_told - Tolerance for second adjoint solve
564: Level: beginner
565: M*/
566: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
567: {
568: TAO_LCL *lclP;
569: const char *morethuente_type = TAOLINESEARCHMT;
571: tao->ops->setup = TaoSetup_LCL;
572: tao->ops->solve = TaoSolve_LCL;
573: tao->ops->view = TaoView_LCL;
574: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
575: tao->ops->destroy = TaoDestroy_LCL;
576: PetscNew(&lclP);
577: tao->data = (void *)lclP;
579: /* Override default settings (unless already changed) */
580: if (!tao->max_it_changed) tao->max_it = 200;
581: if (!tao->catol_changed) tao->catol = 1.0e-4;
582: if (!tao->gatol_changed) tao->gttol = 1.0e-4;
583: if (!tao->grtol_changed) tao->gttol = 1.0e-4;
584: if (!tao->gttol_changed) tao->gttol = 1.0e-4;
585: lclP->rho0 = 1.0e-4;
586: lclP->rhomax = 1e5;
587: lclP->eps1 = 1.0e-8;
588: lclP->eps2 = 0.0;
589: lclP->solve_type = 2;
590: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
591: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
592: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
593: TaoLineSearchSetType(tao->linesearch, morethuente_type);
594: TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
596: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao);
597: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
598: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
599: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
600: KSPSetFromOptions(tao->ksp);
602: MatCreate(((PetscObject)tao)->comm, &lclP->R);
603: MatSetType(lclP->R, MATLMVMBFGS);
604: return 0;
605: }
607: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
608: {
609: Tao tao = (Tao)ptr;
610: TAO_LCL *lclP = (TAO_LCL *)tao->data;
611: PetscBool set, pset, flag, pflag, symmetric;
612: PetscReal cdotl;
614: TaoComputeObjectiveAndGradient(tao, X, f, G);
615: LCLScatter(lclP, G, lclP->GU, lclP->GV);
616: if (lclP->recompute_jacobian_flag) {
617: TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv);
618: TaoComputeJacobianDesign(tao, X, tao->jacobian_design);
619: }
620: TaoComputeConstraints(tao, X, tao->constraints);
621: MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
622: if (tao->jacobian_state_pre) {
623: MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
624: } else {
625: pset = pflag = PETSC_TRUE;
626: }
627: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
628: else symmetric = PETSC_FALSE;
630: VecDot(lclP->lambda0, tao->constraints, &cdotl);
631: lclP->lgn = *f - cdotl;
633: /* Gradient of Lagrangian GL = G - J' * lambda */
634: /* WU = A' * WL
635: WV = B' * WL */
636: if (symmetric) {
637: MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U);
638: } else {
639: MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U);
640: }
641: MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V);
642: VecScale(lclP->GL_U, -1.0);
643: VecScale(lclP->GL_V, -1.0);
644: VecAXPY(lclP->GL_U, 1.0, lclP->GU);
645: VecAXPY(lclP->GL_V, 1.0, lclP->GV);
646: LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL);
648: f[0] = lclP->lgn;
649: VecCopy(lclP->GL, G);
650: return 0;
651: }
653: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
654: {
655: Tao tao = (Tao)ptr;
656: TAO_LCL *lclP = (TAO_LCL *)tao->data;
657: PetscReal con2;
658: PetscBool flag, pflag, set, pset, symmetric;
660: LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao);
661: LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V);
662: VecDot(tao->constraints, tao->constraints, &con2);
663: lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;
665: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
666: /* WU = A' * c
667: WV = B' * c */
668: MatIsSymmetricKnown(tao->jacobian_state, &set, &flag);
669: if (tao->jacobian_state_pre) {
670: MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag);
671: } else {
672: pset = pflag = PETSC_TRUE;
673: }
674: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
675: else symmetric = PETSC_FALSE;
677: if (symmetric) {
678: MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U);
679: } else {
680: MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U);
681: }
683: MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V);
684: VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U);
685: VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V);
686: LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL);
688: f[0] = lclP->aug;
689: VecCopy(lclP->GAugL, G);
690: return 0;
691: }
693: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
694: {
695: VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
696: VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
697: VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
698: VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
699: return 0;
700: }
701: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
702: {
703: VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
704: VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
705: VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
706: VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
707: return 0;
708: }