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: PetscFunctionBegin;
12: if (tao->setupcalled) {
13: PetscCall(MatDestroy(&lclP->R));
14: PetscCall(VecDestroy(&lclP->lambda));
15: PetscCall(VecDestroy(&lclP->lambda0));
16: PetscCall(VecDestroy(&lclP->WL));
17: PetscCall(VecDestroy(&lclP->W));
18: PetscCall(VecDestroy(&lclP->X0));
19: PetscCall(VecDestroy(&lclP->G0));
20: PetscCall(VecDestroy(&lclP->GL));
21: PetscCall(VecDestroy(&lclP->GAugL));
22: PetscCall(VecDestroy(&lclP->dbar));
23: PetscCall(VecDestroy(&lclP->U));
24: PetscCall(VecDestroy(&lclP->U0));
25: PetscCall(VecDestroy(&lclP->V));
26: PetscCall(VecDestroy(&lclP->V0));
27: PetscCall(VecDestroy(&lclP->V1));
28: PetscCall(VecDestroy(&lclP->GU));
29: PetscCall(VecDestroy(&lclP->GV));
30: PetscCall(VecDestroy(&lclP->GU0));
31: PetscCall(VecDestroy(&lclP->GV0));
32: PetscCall(VecDestroy(&lclP->GL_U));
33: PetscCall(VecDestroy(&lclP->GL_V));
34: PetscCall(VecDestroy(&lclP->GAugL_U));
35: PetscCall(VecDestroy(&lclP->GAugL_V));
36: PetscCall(VecDestroy(&lclP->GL_U0));
37: PetscCall(VecDestroy(&lclP->GL_V0));
38: PetscCall(VecDestroy(&lclP->GAugL_U0));
39: PetscCall(VecDestroy(&lclP->GAugL_V0));
40: PetscCall(VecDestroy(&lclP->DU));
41: PetscCall(VecDestroy(&lclP->DV));
42: PetscCall(VecDestroy(&lclP->WU));
43: PetscCall(VecDestroy(&lclP->WV));
44: PetscCall(VecDestroy(&lclP->g1));
45: PetscCall(VecDestroy(&lclP->g2));
46: PetscCall(VecDestroy(&lclP->con1));
48: PetscCall(VecDestroy(&lclP->r));
49: PetscCall(VecDestroy(&lclP->s));
51: PetscCall(ISDestroy(&tao->state_is));
52: PetscCall(ISDestroy(&tao->design_is));
54: PetscCall(VecScatterDestroy(&lclP->state_scatter));
55: PetscCall(VecScatterDestroy(&lclP->design_scatter));
56: }
57: PetscCall(MatDestroy(&lclP->R));
58: PetscCall(KSPDestroy(&tao->ksp));
59: PetscCall(PetscFree(tao->data));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode TaoSetFromOptions_LCL(Tao tao, PetscOptionItems PetscOptionsObject)
64: {
65: TAO_LCL *lclP = (TAO_LCL *)tao->data;
67: PetscFunctionBegin;
68: PetscOptionsHeadBegin(PetscOptionsObject, "Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
69: PetscCall(PetscOptionsReal("-tao_lcl_eps1", "epsilon 1 tolerance", "", lclP->eps1, &lclP->eps1, NULL));
70: PetscCall(PetscOptionsReal("-tao_lcl_eps2", "epsilon 2 tolerance", "", lclP->eps2, &lclP->eps2, NULL));
71: PetscCall(PetscOptionsReal("-tao_lcl_rho0", "init value for rho", "", lclP->rho0, &lclP->rho0, NULL));
72: PetscCall(PetscOptionsReal("-tao_lcl_rhomax", "max value for rho", "", lclP->rhomax, &lclP->rhomax, NULL));
73: lclP->phase2_niter = 1;
74: PetscCall(PetscOptionsInt("-tao_lcl_phase2_niter", "Number of phase 2 iterations in LCL algorithm", "", lclP->phase2_niter, &lclP->phase2_niter, NULL));
75: lclP->verbose = PETSC_FALSE;
76: PetscCall(PetscOptionsBool("-tao_lcl_verbose", "Print verbose output", "", lclP->verbose, &lclP->verbose, NULL));
77: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
78: PetscCall(PetscOptionsReal("-tao_lcl_tola", "Tolerance for first forward solve", "", lclP->tau[0], &lclP->tau[0], NULL));
79: PetscCall(PetscOptionsReal("-tao_lcl_tolb", "Tolerance for first adjoint solve", "", lclP->tau[1], &lclP->tau[1], NULL));
80: PetscCall(PetscOptionsReal("-tao_lcl_tolc", "Tolerance for second forward solve", "", lclP->tau[2], &lclP->tau[2], NULL));
81: PetscCall(PetscOptionsReal("-tao_lcl_told", "Tolerance for second adjoint solve", "", lclP->tau[3], &lclP->tau[3], NULL));
82: PetscOptionsHeadEnd();
83: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
84: PetscCall(MatSetFromOptions(lclP->R));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
89: {
90: return PETSC_SUCCESS;
91: }
93: static PetscErrorCode TaoSetup_LCL(Tao tao)
94: {
95: TAO_LCL *lclP = (TAO_LCL *)tao->data;
96: PetscInt lo, hi, nlocalstate, nlocaldesign;
97: IS is_state, is_design;
99: PetscFunctionBegin;
100: PetscCheck(tao->state_is, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "LCL Solver requires an initial state index set -- use TaoSetStateIS()");
101: PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
102: PetscCall(VecDuplicate(tao->solution, &lclP->W));
103: PetscCall(VecDuplicate(tao->solution, &lclP->X0));
104: PetscCall(VecDuplicate(tao->solution, &lclP->G0));
105: PetscCall(VecDuplicate(tao->solution, &lclP->GL));
106: PetscCall(VecDuplicate(tao->solution, &lclP->GAugL));
108: PetscCall(VecDuplicate(tao->constraints, &lclP->lambda));
109: PetscCall(VecDuplicate(tao->constraints, &lclP->WL));
110: PetscCall(VecDuplicate(tao->constraints, &lclP->lambda0));
111: PetscCall(VecDuplicate(tao->constraints, &lclP->con1));
113: PetscCall(VecGetSize(tao->solution, &lclP->n));
114: PetscCall(VecGetSize(tao->constraints, &lclP->m));
116: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->U));
117: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->V));
118: PetscCall(ISGetLocalSize(tao->state_is, &nlocalstate));
119: PetscCall(ISGetLocalSize(tao->design_is, &nlocaldesign));
120: PetscCall(VecSetSizes(lclP->U, nlocalstate, lclP->m));
121: PetscCall(VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m));
122: PetscCall(VecSetType(lclP->U, ((PetscObject)tao->solution)->type_name));
123: PetscCall(VecSetType(lclP->V, ((PetscObject)tao->solution)->type_name));
124: PetscCall(VecSetFromOptions(lclP->U));
125: PetscCall(VecSetFromOptions(lclP->V));
126: PetscCall(VecDuplicate(lclP->U, &lclP->DU));
127: PetscCall(VecDuplicate(lclP->U, &lclP->U0));
128: PetscCall(VecDuplicate(lclP->U, &lclP->GU));
129: PetscCall(VecDuplicate(lclP->U, &lclP->GU0));
130: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U));
131: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U));
132: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U0));
133: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U0));
134: PetscCall(VecDuplicate(lclP->U, &lclP->WU));
135: PetscCall(VecDuplicate(lclP->U, &lclP->r));
136: PetscCall(VecDuplicate(lclP->V, &lclP->V0));
137: PetscCall(VecDuplicate(lclP->V, &lclP->V1));
138: PetscCall(VecDuplicate(lclP->V, &lclP->DV));
139: PetscCall(VecDuplicate(lclP->V, &lclP->s));
140: PetscCall(VecDuplicate(lclP->V, &lclP->GV));
141: PetscCall(VecDuplicate(lclP->V, &lclP->GV0));
142: PetscCall(VecDuplicate(lclP->V, &lclP->dbar));
143: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V));
144: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V));
145: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V0));
146: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V0));
147: PetscCall(VecDuplicate(lclP->V, &lclP->WV));
148: PetscCall(VecDuplicate(lclP->V, &lclP->g1));
149: PetscCall(VecDuplicate(lclP->V, &lclP->g2));
151: /* create scatters for state, design subvecs */
152: PetscCall(VecGetOwnershipRange(lclP->U, &lo, &hi));
153: PetscCall(ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state));
154: PetscCall(VecGetOwnershipRange(lclP->V, &lo, &hi));
155: if (0) {
156: PetscInt sizeU, sizeV;
157: PetscCall(VecGetSize(lclP->U, &sizeU));
158: PetscCall(VecGetSize(lclP->V, &sizeV));
159: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV));
160: }
161: PetscCall(ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design));
162: PetscCall(VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter));
163: PetscCall(VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter));
164: PetscCall(ISDestroy(&is_state));
165: PetscCall(ISDestroy(&is_design));
166: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
180: lclP->rho = lclP->rho0;
181: PetscCall(VecGetLocalSize(lclP->U, &nlocal));
182: PetscCall(VecGetLocalSize(lclP->V, &nlocal));
183: PetscCall(MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m));
184: PetscCall(MatLMVMAllocate(lclP->R, lclP->V, lclP->V));
185: lclP->recompute_jacobian_flag = PETSC_TRUE;
187: /* Scatter to U,V */
188: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
190: /* Evaluate Function, Gradient, Constraints, and Jacobian */
191: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
192: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
193: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
194: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
196: /* Scatter gradient to GU,GV */
197: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
199: /* Evaluate Lagrangian function and gradient */
200: /* p0 */
201: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
202: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
203: if (tao->jacobian_state_pre) {
204: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
205: } else {
206: pset = pflag = PETSC_TRUE;
207: }
208: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
209: else symmetric = PETSC_FALSE;
211: lclP->solve_type = LCL_ADJOINT2;
212: if (tao->jacobian_state_inv) {
213: if (symmetric) {
214: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
215: } else {
216: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
217: }
218: } else {
219: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
220: if (symmetric) {
221: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
222: } else {
223: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
224: }
225: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
226: tao->ksp_its += its;
227: tao->ksp_tot_its += its;
228: }
229: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
230: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
232: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
233: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
235: /* Evaluate constraint norm */
236: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
237: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
239: /* Monitor convergence */
240: tao->reason = TAO_CONTINUE_ITERATING;
241: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
242: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
243: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
245: while (tao->reason == TAO_CONTINUE_ITERATING) {
246: /* Call general purpose update function */
247: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
248: tao->ksp_its = 0;
249: /* Compute a descent direction for the linearly constrained subproblem
250: minimize f(u+du, v+dv)
251: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
253: /* Store the points around the linearization */
254: PetscCall(VecCopy(lclP->U, lclP->U0));
255: PetscCall(VecCopy(lclP->V, lclP->V0));
256: PetscCall(VecCopy(lclP->GU, lclP->GU0));
257: PetscCall(VecCopy(lclP->GV, lclP->GV0));
258: PetscCall(VecCopy(lclP->GAugL_U, lclP->GAugL_U0));
259: PetscCall(VecCopy(lclP->GAugL_V, lclP->GAugL_V0));
260: PetscCall(VecCopy(lclP->GL_U, lclP->GL_U0));
261: PetscCall(VecCopy(lclP->GL_V, lclP->GL_V0));
263: lclP->aug0 = lclP->aug;
264: lclP->lgn0 = lclP->lgn;
266: /* Given the design variables, we need to project the current iterate
267: onto the linearized constraint. We choose to fix the design variables
268: and solve the linear system for the state variables. The resulting
269: point is the Newton direction */
271: /* Solve r = A\con */
272: lclP->solve_type = LCL_FORWARD1;
273: PetscCall(VecSet(lclP->r, 0.0)); /* Initial guess in CG */
275: if (tao->jacobian_state_inv) {
276: PetscCall(MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r));
277: } else {
278: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
279: PetscCall(KSPSolve(tao->ksp, tao->constraints, lclP->r));
280: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
281: tao->ksp_its += its;
282: tao->ksp_tot_its += tao->ksp_its;
283: }
285: /* Set design step direction dv to zero */
286: PetscCall(VecSet(lclP->s, 0.0));
288: /*
289: Check sufficient descent for constraint merit function .5*||con||^2
290: con' Ak r >= eps1 ||r||^(2+eps2)
291: */
293: /* Compute WU= Ak' * con */
294: if (symmetric) {
295: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->WU));
296: } else {
297: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->WU));
298: }
299: /* Compute r * Ak' * con */
300: PetscCall(VecDot(lclP->r, lclP->WU, &rWU));
302: /* compute ||r||^(2+eps2) */
303: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
304: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
305: adec = lclP->eps1 * r2;
307: if (rWU < adec) {
308: PetscCall(PetscInfo(tao, "Newton direction not descent for constraint, feasibility phase required\n"));
309: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent));
311: PetscCall(PetscInfo(tao, "Using steepest descent direction instead.\n"));
312: PetscCall(VecSet(lclP->r, 0.0));
313: PetscCall(VecAXPY(lclP->r, -1.0, lclP->WU));
314: PetscCall(VecDot(lclP->r, lclP->r, &rWU));
315: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
316: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
317: PetscCall(VecDot(lclP->r, lclP->GAugL_U, &descent));
318: adec = lclP->eps1 * r2;
319: }
321: /*
322: Check descent for aug. lagrangian
323: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
324: GL_U = GUk - Ak'*yk
325: WU = Ak'*con
326: adec=eps1||r||^(2+eps2)
328: ==>
329: Check r'GL_U - rho*r'WU <= adec
330: */
332: PetscCall(VecDot(lclP->r, lclP->GL_U, &rGL_U));
333: aldescent = rGL_U - lclP->rho * rWU;
334: if (aldescent > -adec) {
335: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent));
336: PetscCall(PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent));
337: lclP->rho = (rGL_U - adec) / rWU;
338: PetscCheck(lclP->rho <= lclP->rhomax, PETSC_COMM_WORLD, PETSC_ERR_SUP, "rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)", (double)lclP->rho);
339: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Increasing penalty parameter to %g\n", (double)lclP->rho));
340: PetscCall(PetscInfo(tao, " Increasing penalty parameter to %g\n", (double)lclP->rho));
341: }
343: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
344: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
345: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
347: /* We now minimize the augmented Lagrangian along the Newton direction */
348: PetscCall(VecScale(lclP->r, -1.0));
349: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
350: PetscCall(VecScale(lclP->r, -1.0));
351: PetscCall(LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL));
352: PetscCall(LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0));
354: lclP->recompute_jacobian_flag = PETSC_TRUE;
356: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
357: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
358: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
359: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
360: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step));
362: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
363: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
364: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
366: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
368: /* Check convergence */
369: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
370: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
371: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
372: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
373: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
374: if (tao->reason != TAO_CONTINUE_ITERATING) break;
376: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
377: for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
378: /* We now minimize the objective function starting from the fraction of
379: the Newton point accepted by applying one step of a reduced-space
380: method. The optimization problem is:
382: minimize f(u+du, v+dv)
383: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
385: In particular, we have that
386: du = -inv(A)*(Bdv + alpha g) */
388: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
389: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
391: /* Store V and constraints */
392: PetscCall(VecCopy(lclP->V, lclP->V1));
393: PetscCall(VecCopy(tao->constraints, lclP->con1));
395: /* Compute multipliers */
396: /* p1 */
397: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
398: lclP->solve_type = LCL_ADJOINT1;
399: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
400: if (tao->jacobian_state_pre) {
401: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
402: } else {
403: pset = pflag = PETSC_TRUE;
404: }
405: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
406: else symmetric = PETSC_FALSE;
408: if (tao->jacobian_state_inv) {
409: if (symmetric) {
410: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
411: } else {
412: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
413: }
414: } else {
415: if (symmetric) {
416: PetscCall(KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda));
417: } else {
418: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda));
419: }
420: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
421: tao->ksp_its += its;
422: tao->ksp_tot_its += its;
423: }
424: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1));
425: PetscCall(VecAXPY(lclP->g1, -1.0, lclP->GAugL_V));
427: /* Compute the limited-memory quasi-newton direction */
428: if (tao->niter > 0) {
429: PetscCall(MatSolve(lclP->R, lclP->g1, lclP->s));
430: PetscCall(VecDot(lclP->s, lclP->g1, &descent));
431: if (descent <= 0) {
432: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent));
433: PetscCall(VecCopy(lclP->g1, lclP->s));
434: }
435: } else {
436: PetscCall(VecCopy(lclP->g1, lclP->s));
437: }
438: PetscCall(VecScale(lclP->g1, -1.0));
440: /* Recover the full space direction */
441: PetscCall(MatMult(tao->jacobian_design, lclP->s, lclP->WU));
442: /* PetscCall(VecSet(lclP->r,0.0)); */ /* Initial guess in CG */
443: lclP->solve_type = LCL_FORWARD2;
444: if (tao->jacobian_state_inv) {
445: PetscCall(MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r));
446: } else {
447: PetscCall(KSPSolve(tao->ksp, lclP->WU, lclP->r));
448: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
449: tao->ksp_its += its;
450: tao->ksp_tot_its += its;
451: }
453: /* We now minimize the augmented Lagrangian along the direction -r,s */
454: PetscCall(VecScale(lclP->r, -1.0));
455: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
456: PetscCall(VecScale(lclP->r, -1.0));
457: lclP->recompute_jacobian_flag = PETSC_TRUE;
459: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
460: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
461: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
462: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
463: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength = %g\n", (double)step));
465: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
466: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
467: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
468: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
469: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
471: /* Compute the reduced gradient at the new point */
473: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
474: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
476: /* p2 */
477: /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
478: if (phase2_iter == 0) {
479: PetscCall(VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0));
480: } else {
481: PetscCall(VecAXPY(lclP->lambda, -lclP->rho, tao->constraints));
482: }
484: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
485: if (tao->jacobian_state_pre) {
486: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
487: } else {
488: pset = pflag = PETSC_TRUE;
489: }
490: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
491: else symmetric = PETSC_FALSE;
493: lclP->solve_type = LCL_ADJOINT2;
494: if (tao->jacobian_state_inv) {
495: if (symmetric) {
496: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
497: } else {
498: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
499: }
500: } else {
501: if (symmetric) {
502: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
503: } else {
504: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
505: }
506: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
507: tao->ksp_its += its;
508: tao->ksp_tot_its += its;
509: }
511: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2));
512: PetscCall(VecAXPY(lclP->g2, -1.0, lclP->GV));
514: PetscCall(VecScale(lclP->g2, -1.0));
516: /* Update the quasi-newton approximation */
517: PetscCall(MatLMVMUpdate(lclP->R, lclP->V, lclP->g2));
518: /* 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 */
519: }
521: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
523: /* Evaluate Function, Gradient, Constraints, and Jacobian */
524: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
525: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
526: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
528: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
529: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
530: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
532: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
534: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
536: /* Evaluate constraint norm */
537: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
539: /* Monitor convergence */
540: tao->niter++;
541: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
542: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
543: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
544: }
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: /*MC
549: TAOLCL - linearly constrained Lagrangian method for PDE-constrained optimization
551: Option Database Keys:
552: + -tao_lcl_eps1 - epsilon 1 tolerance
553: . -tao_lcl_eps2 - epsilon 2 tolerance
554: . -tao_lcl_rho0 - initial value for rho
555: . -tao_lcl_rhomax - maximum allowed value for rho
556: . -tao_lcl_phase2_niter - Number of phase 2 iterations in the LCL algorithm
557: . -tao_lcl_verbose - Print verbose output if True
558: . -tao_lcl_tola - Tolerance for first forward solve
559: . -tao_lcl_tolb - Tolerance for first adjoint solve
560: . -tao_lcl_tolc - Tolerance for second forward solve
561: - -tao_lcl_told - Tolerance for second adjoint solve
563: Level: beginner
565: .seealso: `Tao`, `TaoType`, `TaoSetStateDesignIS()`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
566: M*/
567: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
568: {
569: TAO_LCL *lclP;
570: const char *morethuente_type = TAOLINESEARCHMT;
572: PetscFunctionBegin;
573: tao->ops->setup = TaoSetup_LCL;
574: tao->ops->solve = TaoSolve_LCL;
575: tao->ops->view = TaoView_LCL;
576: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
577: tao->ops->destroy = TaoDestroy_LCL;
578: PetscCall(PetscNew(&lclP));
579: tao->data = (void *)lclP;
580: tao->uses_gradient = PETSC_TRUE;
582: /* Override default settings (unless already changed) */
583: PetscCall(TaoParametersInitialize(tao));
584: PetscObjectParameterSetDefault(tao, max_it, 200);
585: PetscObjectParameterSetDefault(tao, catol, 1.0e-4);
586: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
587: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
588: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
590: lclP->rho0 = 1.0e-4;
591: lclP->rhomax = 1e5;
592: lclP->eps1 = 1.0e-8;
593: lclP->eps2 = 0.0;
594: lclP->solve_type = 2;
595: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
596: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
597: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
598: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
599: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
601: PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao));
602: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
603: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
604: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
605: PetscCall(KSPSetFromOptions(tao->ksp));
607: PetscCall(MatCreate(((PetscObject)tao)->comm, &lclP->R));
608: PetscCall(MatSetType(lclP->R, MATLMVMBFGS));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
613: {
614: Tao tao = (Tao)ptr;
615: TAO_LCL *lclP = (TAO_LCL *)tao->data;
616: PetscBool set, pset, flag, pflag, symmetric;
617: PetscReal cdotl;
619: PetscFunctionBegin;
620: PetscCall(TaoComputeObjectiveAndGradient(tao, X, f, G));
621: PetscCall(LCLScatter(lclP, G, lclP->GU, lclP->GV));
622: if (lclP->recompute_jacobian_flag) {
623: PetscCall(TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
624: PetscCall(TaoComputeJacobianDesign(tao, X, tao->jacobian_design));
625: }
626: PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
627: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
628: if (tao->jacobian_state_pre) {
629: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
630: } else {
631: pset = pflag = PETSC_TRUE;
632: }
633: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
634: else symmetric = PETSC_FALSE;
636: PetscCall(VecDot(lclP->lambda0, tao->constraints, &cdotl));
637: lclP->lgn = *f - cdotl;
639: /* Gradient of Lagrangian GL = G - J' * lambda */
640: /* WU = A' * WL
641: WV = B' * WL */
642: if (symmetric) {
643: PetscCall(MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
644: } else {
645: PetscCall(MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
646: }
647: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V));
648: PetscCall(VecScale(lclP->GL_U, -1.0));
649: PetscCall(VecScale(lclP->GL_V, -1.0));
650: PetscCall(VecAXPY(lclP->GL_U, 1.0, lclP->GU));
651: PetscCall(VecAXPY(lclP->GL_V, 1.0, lclP->GV));
652: PetscCall(LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL));
654: f[0] = lclP->lgn;
655: PetscCall(VecCopy(lclP->GL, G));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
660: {
661: Tao tao = (Tao)ptr;
662: TAO_LCL *lclP = (TAO_LCL *)tao->data;
663: PetscReal con2;
664: PetscBool flag, pflag, set, pset, symmetric;
666: PetscFunctionBegin;
667: PetscCall(LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao));
668: PetscCall(LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V));
669: PetscCall(VecDot(tao->constraints, tao->constraints, &con2));
670: lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;
672: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
673: /* WU = A' * c
674: WV = B' * c */
675: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
676: if (tao->jacobian_state_pre) {
677: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
678: } else {
679: pset = pflag = PETSC_TRUE;
680: }
681: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
682: else symmetric = PETSC_FALSE;
684: if (symmetric) {
685: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
686: } else {
687: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
688: }
690: PetscCall(MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V));
691: PetscCall(VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U));
692: PetscCall(VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V));
693: PetscCall(LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL));
695: f[0] = lclP->aug;
696: PetscCall(VecCopy(lclP->GAugL, G));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
701: {
702: PetscFunctionBegin;
703: PetscCall(VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
704: PetscCall(VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
705: PetscCall(VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
706: PetscCall(VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
709: static PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
710: {
711: PetscFunctionBegin;
712: PetscCall(VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
713: PetscCall(VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
714: PetscCall(VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
715: PetscCall(VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }