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(VecSet(lclP->lambda, 0.0));
115: PetscCall(VecGetSize(tao->solution, &lclP->n));
116: PetscCall(VecGetSize(tao->constraints, &lclP->m));
118: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->U));
119: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->V));
120: PetscCall(ISGetLocalSize(tao->state_is, &nlocalstate));
121: PetscCall(ISGetLocalSize(tao->design_is, &nlocaldesign));
122: PetscCall(VecSetSizes(lclP->U, nlocalstate, lclP->m));
123: PetscCall(VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m));
124: PetscCall(VecSetType(lclP->U, ((PetscObject)tao->solution)->type_name));
125: PetscCall(VecSetType(lclP->V, ((PetscObject)tao->solution)->type_name));
126: PetscCall(VecSetFromOptions(lclP->U));
127: PetscCall(VecSetFromOptions(lclP->V));
128: PetscCall(VecDuplicate(lclP->U, &lclP->DU));
129: PetscCall(VecDuplicate(lclP->U, &lclP->U0));
130: PetscCall(VecDuplicate(lclP->U, &lclP->GU));
131: PetscCall(VecDuplicate(lclP->U, &lclP->GU0));
132: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U));
133: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U));
134: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U0));
135: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U0));
136: PetscCall(VecDuplicate(lclP->U, &lclP->WU));
137: PetscCall(VecDuplicate(lclP->U, &lclP->r));
138: PetscCall(VecDuplicate(lclP->V, &lclP->V0));
139: PetscCall(VecDuplicate(lclP->V, &lclP->V1));
140: PetscCall(VecDuplicate(lclP->V, &lclP->DV));
141: PetscCall(VecDuplicate(lclP->V, &lclP->s));
142: PetscCall(VecDuplicate(lclP->V, &lclP->GV));
143: PetscCall(VecDuplicate(lclP->V, &lclP->GV0));
144: PetscCall(VecDuplicate(lclP->V, &lclP->dbar));
145: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V));
146: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V));
147: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V0));
148: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V0));
149: PetscCall(VecDuplicate(lclP->V, &lclP->WV));
150: PetscCall(VecDuplicate(lclP->V, &lclP->g1));
151: PetscCall(VecDuplicate(lclP->V, &lclP->g2));
153: /* create scatters for state, design subvecs */
154: PetscCall(VecGetOwnershipRange(lclP->U, &lo, &hi));
155: PetscCall(ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state));
156: PetscCall(VecGetOwnershipRange(lclP->V, &lo, &hi));
157: if (0) {
158: PetscInt sizeU, sizeV;
159: PetscCall(VecGetSize(lclP->U, &sizeU));
160: PetscCall(VecGetSize(lclP->V, &sizeV));
161: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV));
162: }
163: PetscCall(ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design));
164: PetscCall(VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter));
165: PetscCall(VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter));
166: PetscCall(ISDestroy(&is_state));
167: PetscCall(ISDestroy(&is_design));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode TaoSolve_LCL(Tao tao)
172: {
173: TAO_LCL *lclP = (TAO_LCL *)tao->data;
174: PetscInt phase2_iter, nlocal, its;
175: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
176: PetscReal step = 1.0, f, descent, aldescent;
177: PetscReal cnorm, mnorm;
178: PetscReal adec, r2, rGL_U, rWU;
179: PetscBool set, pset, flag, pflag, symmetric;
181: PetscFunctionBegin;
182: lclP->rho = lclP->rho0;
183: PetscCall(VecGetLocalSize(lclP->U, &nlocal));
184: PetscCall(VecGetLocalSize(lclP->V, &nlocal));
185: PetscCall(MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m));
186: PetscCall(MatLMVMAllocate(lclP->R, lclP->V, lclP->V));
187: lclP->recompute_jacobian_flag = PETSC_TRUE;
189: /* Scatter to U,V */
190: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
192: /* Evaluate Function, Gradient, Constraints, and Jacobian */
193: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
194: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
195: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
196: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
198: /* Scatter gradient to GU,GV */
199: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
201: /* Evaluate Lagrangian function and gradient */
202: /* p0 */
203: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
204: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
205: if (tao->jacobian_state_pre) {
206: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
207: } else {
208: pset = pflag = PETSC_TRUE;
209: }
210: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
211: else symmetric = PETSC_FALSE;
213: lclP->solve_type = LCL_ADJOINT2;
214: if (tao->jacobian_state_inv) {
215: if (symmetric) {
216: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
217: } else {
218: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
219: }
220: } else {
221: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
222: if (symmetric) {
223: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
224: } else {
225: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
226: }
227: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
228: tao->ksp_its += its;
229: tao->ksp_tot_its += its;
230: }
231: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
232: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
234: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
235: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
237: /* Evaluate constraint norm */
238: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
239: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
241: /* Monitor convergence */
242: tao->reason = TAO_CONTINUE_ITERATING;
243: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
244: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
245: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
247: while (tao->reason == TAO_CONTINUE_ITERATING) {
248: /* Call general purpose update function */
249: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
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: PetscCall(VecCopy(lclP->U, lclP->U0));
257: PetscCall(VecCopy(lclP->V, lclP->V0));
258: PetscCall(VecCopy(lclP->GU, lclP->GU0));
259: PetscCall(VecCopy(lclP->GV, lclP->GV0));
260: PetscCall(VecCopy(lclP->GAugL_U, lclP->GAugL_U0));
261: PetscCall(VecCopy(lclP->GAugL_V, lclP->GAugL_V0));
262: PetscCall(VecCopy(lclP->GL_U, lclP->GL_U0));
263: PetscCall(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: PetscCall(VecSet(lclP->r, 0.0)); /* Initial guess in CG */
277: if (tao->jacobian_state_inv) {
278: PetscCall(MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r));
279: } else {
280: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
281: PetscCall(KSPSolve(tao->ksp, tao->constraints, lclP->r));
282: PetscCall(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: PetscCall(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: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->WU));
298: } else {
299: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->WU));
300: }
301: /* Compute r * Ak' * con */
302: PetscCall(VecDot(lclP->r, lclP->WU, &rWU));
304: /* compute ||r||^(2+eps2) */
305: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
306: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
307: adec = lclP->eps1 * r2;
309: if (rWU < adec) {
310: PetscCall(PetscInfo(tao, "Newton direction not descent for constraint, feasibility phase required\n"));
311: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent));
313: PetscCall(PetscInfo(tao, "Using steepest descent direction instead.\n"));
314: PetscCall(VecSet(lclP->r, 0.0));
315: PetscCall(VecAXPY(lclP->r, -1.0, lclP->WU));
316: PetscCall(VecDot(lclP->r, lclP->r, &rWU));
317: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
318: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
319: PetscCall(VecDot(lclP->r, lclP->GAugL_U, &descent));
320: adec = lclP->eps1 * r2;
321: }
323: /*
324: Check descent for aug. lagrangian
325: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
326: GL_U = GUk - Ak'*yk
327: WU = Ak'*con
328: adec=eps1||r||^(2+eps2)
330: ==>
331: Check r'GL_U - rho*r'WU <= adec
332: */
334: PetscCall(VecDot(lclP->r, lclP->GL_U, &rGL_U));
335: aldescent = rGL_U - lclP->rho * rWU;
336: if (aldescent > -adec) {
337: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent));
338: PetscCall(PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent));
339: lclP->rho = (rGL_U - adec) / rWU;
340: 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);
341: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Increasing penalty parameter to %g\n", (double)lclP->rho));
342: PetscCall(PetscInfo(tao, " Increasing penalty parameter to %g\n", (double)lclP->rho));
343: }
345: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
346: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
347: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
349: /* We now minimize the augmented Lagrangian along the Newton direction */
350: PetscCall(VecScale(lclP->r, -1.0));
351: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
352: PetscCall(VecScale(lclP->r, -1.0));
353: PetscCall(LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL));
354: PetscCall(LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0));
356: lclP->recompute_jacobian_flag = PETSC_TRUE;
358: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
359: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
360: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
361: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
362: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step));
364: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
365: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
366: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
368: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
370: /* Check convergence */
371: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
372: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
373: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
374: PetscCall(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: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
391: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
393: /* Store V and constraints */
394: PetscCall(VecCopy(lclP->V, lclP->V1));
395: PetscCall(VecCopy(tao->constraints, lclP->con1));
397: /* Compute multipliers */
398: /* p1 */
399: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
400: lclP->solve_type = LCL_ADJOINT1;
401: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
402: if (tao->jacobian_state_pre) {
403: PetscCall(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: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
413: } else {
414: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
415: }
416: } else {
417: if (symmetric) {
418: PetscCall(KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda));
419: } else {
420: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda));
421: }
422: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
423: tao->ksp_its += its;
424: tao->ksp_tot_its += its;
425: }
426: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1));
427: PetscCall(VecAXPY(lclP->g1, -1.0, lclP->GAugL_V));
429: /* Compute the limited-memory quasi-newton direction */
430: if (tao->niter > 0) {
431: PetscCall(MatSolve(lclP->R, lclP->g1, lclP->s));
432: PetscCall(VecDot(lclP->s, lclP->g1, &descent));
433: if (descent <= 0) {
434: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent));
435: PetscCall(VecCopy(lclP->g1, lclP->s));
436: }
437: } else {
438: PetscCall(VecCopy(lclP->g1, lclP->s));
439: }
440: PetscCall(VecScale(lclP->g1, -1.0));
442: /* Recover the full space direction */
443: PetscCall(MatMult(tao->jacobian_design, lclP->s, lclP->WU));
444: /* PetscCall(VecSet(lclP->r,0.0)); */ /* Initial guess in CG */
445: lclP->solve_type = LCL_FORWARD2;
446: if (tao->jacobian_state_inv) {
447: PetscCall(MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r));
448: } else {
449: PetscCall(KSPSolve(tao->ksp, lclP->WU, lclP->r));
450: PetscCall(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: PetscCall(VecScale(lclP->r, -1.0));
457: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
458: PetscCall(VecScale(lclP->r, -1.0));
459: lclP->recompute_jacobian_flag = PETSC_TRUE;
461: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
462: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
463: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
464: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
465: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength = %g\n", (double)step));
467: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
468: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
469: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
470: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
471: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
473: /* Compute the reduced gradient at the new point */
475: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
476: PetscCall(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: PetscCall(VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0));
482: } else {
483: PetscCall(VecAXPY(lclP->lambda, -lclP->rho, tao->constraints));
484: }
486: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
487: if (tao->jacobian_state_pre) {
488: PetscCall(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: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
499: } else {
500: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
501: }
502: } else {
503: if (symmetric) {
504: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
505: } else {
506: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
507: }
508: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
509: tao->ksp_its += its;
510: tao->ksp_tot_its += its;
511: }
513: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2));
514: PetscCall(VecAXPY(lclP->g2, -1.0, lclP->GV));
516: PetscCall(VecScale(lclP->g2, -1.0));
518: /* Update the quasi-newton approximation */
519: PetscCall(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: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
525: /* Evaluate Function, Gradient, Constraints, and Jacobian */
526: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
527: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
528: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
530: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
531: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
532: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
534: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
536: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
538: /* Evaluate constraint norm */
539: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
541: /* Monitor convergence */
542: tao->niter++;
543: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
544: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
545: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
546: }
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: /*MC
551: TAOLCL - linearly constrained Lagrangian method for PDE-constrained optimization
553: Option Database Keys:
554: + -tao_lcl_eps1 - epsilon 1 tolerance
555: . -tao_lcl_eps2 - epsilon 2 tolerance
556: . -tao_lcl_rho0 - initial value for rho
557: . -tao_lcl_rhomax - maximum allowed value for rho
558: . -tao_lcl_phase2_niter - Number of phase 2 iterations in the LCL algorithm
559: . -tao_lcl_verbose - Print verbose output if True
560: . -tao_lcl_tola - Tolerance for first forward solve
561: . -tao_lcl_tolb - Tolerance for first adjoint solve
562: . -tao_lcl_tolc - Tolerance for second forward solve
563: - -tao_lcl_told - Tolerance for second adjoint solve
565: Level: beginner
567: .seealso: `Tao`, `TaoType`, `TaoSetStateDesignIS()`, `TaoSetJacobianStateRoutine()`, `TaoSetJacobianDesignRoutine()`
568: M*/
569: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
570: {
571: TAO_LCL *lclP;
572: const char *morethuente_type = TAOLINESEARCHMT;
574: PetscFunctionBegin;
575: tao->ops->setup = TaoSetup_LCL;
576: tao->ops->solve = TaoSolve_LCL;
577: tao->ops->view = TaoView_LCL;
578: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
579: tao->ops->destroy = TaoDestroy_LCL;
580: PetscCall(PetscNew(&lclP));
581: tao->data = (void *)lclP;
582: tao->uses_gradient = PETSC_TRUE;
584: /* Override default settings (unless already changed) */
585: PetscCall(TaoParametersInitialize(tao));
586: PetscObjectParameterSetDefault(tao, max_it, 200);
587: PetscObjectParameterSetDefault(tao, catol, 1.0e-4);
588: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
589: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
590: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
592: lclP->rho0 = 1.0e-4;
593: lclP->rhomax = 1e5;
594: lclP->eps1 = 1.0e-8;
595: lclP->eps2 = 0.0;
596: lclP->solve_type = 2;
597: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
598: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
599: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
600: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
601: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
603: PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao));
604: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
605: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
606: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
607: PetscCall(KSPSetFromOptions(tao->ksp));
609: PetscCall(MatCreate(((PetscObject)tao)->comm, &lclP->R));
610: PetscCall(MatSetType(lclP->R, MATLMVMBFGS));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
615: {
616: Tao tao = (Tao)ptr;
617: TAO_LCL *lclP = (TAO_LCL *)tao->data;
618: PetscBool set, pset, flag, pflag, symmetric;
619: PetscReal cdotl;
621: PetscFunctionBegin;
622: PetscCall(TaoComputeObjectiveAndGradient(tao, X, f, G));
623: PetscCall(LCLScatter(lclP, G, lclP->GU, lclP->GV));
624: if (lclP->recompute_jacobian_flag) {
625: PetscCall(TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
626: PetscCall(TaoComputeJacobianDesign(tao, X, tao->jacobian_design));
627: }
628: PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
629: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
630: if (tao->jacobian_state_pre) {
631: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
632: } else {
633: pset = pflag = PETSC_TRUE;
634: }
635: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
636: else symmetric = PETSC_FALSE;
638: PetscCall(VecDot(lclP->lambda0, tao->constraints, &cdotl));
639: lclP->lgn = *f - cdotl;
641: /* Gradient of Lagrangian GL = G - J' * lambda */
642: /* WU = A' * WL
643: WV = B' * WL */
644: if (symmetric) {
645: PetscCall(MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
646: } else {
647: PetscCall(MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
648: }
649: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V));
650: PetscCall(VecScale(lclP->GL_U, -1.0));
651: PetscCall(VecScale(lclP->GL_V, -1.0));
652: PetscCall(VecAXPY(lclP->GL_U, 1.0, lclP->GU));
653: PetscCall(VecAXPY(lclP->GL_V, 1.0, lclP->GV));
654: PetscCall(LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL));
656: f[0] = lclP->lgn;
657: PetscCall(VecCopy(lclP->GL, G));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
662: {
663: Tao tao = (Tao)ptr;
664: TAO_LCL *lclP = (TAO_LCL *)tao->data;
665: PetscReal con2;
666: PetscBool flag, pflag, set, pset, symmetric;
668: PetscFunctionBegin;
669: PetscCall(LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao));
670: PetscCall(LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V));
671: PetscCall(VecDot(tao->constraints, tao->constraints, &con2));
672: lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;
674: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
675: /* WU = A' * c
676: WV = B' * c */
677: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
678: if (tao->jacobian_state_pre) {
679: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
680: } else {
681: pset = pflag = PETSC_TRUE;
682: }
683: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
684: else symmetric = PETSC_FALSE;
686: if (symmetric) {
687: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
688: } else {
689: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
690: }
692: PetscCall(MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V));
693: PetscCall(VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U));
694: PetscCall(VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V));
695: PetscCall(LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL));
697: f[0] = lclP->aug;
698: PetscCall(VecCopy(lclP->GAugL, G));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: static PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
703: {
704: PetscFunctionBegin;
705: PetscCall(VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
706: PetscCall(VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
707: PetscCall(VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
708: PetscCall(VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
711: static PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
712: {
713: PetscFunctionBegin;
714: PetscCall(VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
715: PetscCall(VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
716: PetscCall(VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
717: PetscCall(VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }