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->gradient));
102: PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
103: PetscCall(VecDuplicate(tao->solution, &lclP->W));
104: PetscCall(VecDuplicate(tao->solution, &lclP->X0));
105: PetscCall(VecDuplicate(tao->solution, &lclP->G0));
106: PetscCall(VecDuplicate(tao->solution, &lclP->GL));
107: PetscCall(VecDuplicate(tao->solution, &lclP->GAugL));
109: PetscCall(VecDuplicate(tao->constraints, &lclP->lambda));
110: PetscCall(VecDuplicate(tao->constraints, &lclP->WL));
111: PetscCall(VecDuplicate(tao->constraints, &lclP->lambda0));
112: PetscCall(VecDuplicate(tao->constraints, &lclP->con1));
114: PetscCall(VecSet(lclP->lambda, 0.0));
116: PetscCall(VecGetSize(tao->solution, &lclP->n));
117: PetscCall(VecGetSize(tao->constraints, &lclP->m));
119: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->U));
120: PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->V));
121: PetscCall(ISGetLocalSize(tao->state_is, &nlocalstate));
122: PetscCall(ISGetLocalSize(tao->design_is, &nlocaldesign));
123: PetscCall(VecSetSizes(lclP->U, nlocalstate, lclP->m));
124: PetscCall(VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m));
125: PetscCall(VecSetType(lclP->U, ((PetscObject)tao->solution)->type_name));
126: PetscCall(VecSetType(lclP->V, ((PetscObject)tao->solution)->type_name));
127: PetscCall(VecSetFromOptions(lclP->U));
128: PetscCall(VecSetFromOptions(lclP->V));
129: PetscCall(VecDuplicate(lclP->U, &lclP->DU));
130: PetscCall(VecDuplicate(lclP->U, &lclP->U0));
131: PetscCall(VecDuplicate(lclP->U, &lclP->GU));
132: PetscCall(VecDuplicate(lclP->U, &lclP->GU0));
133: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U));
134: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U));
135: PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U0));
136: PetscCall(VecDuplicate(lclP->U, &lclP->GL_U0));
137: PetscCall(VecDuplicate(lclP->U, &lclP->WU));
138: PetscCall(VecDuplicate(lclP->U, &lclP->r));
139: PetscCall(VecDuplicate(lclP->V, &lclP->V0));
140: PetscCall(VecDuplicate(lclP->V, &lclP->V1));
141: PetscCall(VecDuplicate(lclP->V, &lclP->DV));
142: PetscCall(VecDuplicate(lclP->V, &lclP->s));
143: PetscCall(VecDuplicate(lclP->V, &lclP->GV));
144: PetscCall(VecDuplicate(lclP->V, &lclP->GV0));
145: PetscCall(VecDuplicate(lclP->V, &lclP->dbar));
146: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V));
147: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V));
148: PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V0));
149: PetscCall(VecDuplicate(lclP->V, &lclP->GL_V0));
150: PetscCall(VecDuplicate(lclP->V, &lclP->WV));
151: PetscCall(VecDuplicate(lclP->V, &lclP->g1));
152: PetscCall(VecDuplicate(lclP->V, &lclP->g2));
154: /* create scatters for state, design subvecs */
155: PetscCall(VecGetOwnershipRange(lclP->U, &lo, &hi));
156: PetscCall(ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state));
157: PetscCall(VecGetOwnershipRange(lclP->V, &lo, &hi));
158: if (0) {
159: PetscInt sizeU, sizeV;
160: PetscCall(VecGetSize(lclP->U, &sizeU));
161: PetscCall(VecGetSize(lclP->V, &sizeV));
162: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV));
163: }
164: PetscCall(ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design));
165: PetscCall(VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter));
166: PetscCall(VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter));
167: PetscCall(ISDestroy(&is_state));
168: PetscCall(ISDestroy(&is_design));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode TaoSolve_LCL(Tao tao)
173: {
174: TAO_LCL *lclP = (TAO_LCL *)tao->data;
175: PetscInt phase2_iter, nlocal, its;
176: TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
177: PetscReal step = 1.0, f, descent, aldescent;
178: PetscReal cnorm, mnorm;
179: PetscReal adec, r2, rGL_U, rWU;
180: PetscBool set, pset, flag, pflag, symmetric;
182: PetscFunctionBegin;
183: lclP->rho = lclP->rho0;
184: PetscCall(VecGetLocalSize(lclP->U, &nlocal));
185: PetscCall(VecGetLocalSize(lclP->V, &nlocal));
186: PetscCall(MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m));
187: PetscCall(MatLMVMAllocate(lclP->R, lclP->V, lclP->V));
188: lclP->recompute_jacobian_flag = PETSC_TRUE;
190: /* Scatter to U,V */
191: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
193: /* Evaluate Function, Gradient, Constraints, and Jacobian */
194: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
195: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
196: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
197: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
199: /* Scatter gradient to GU,GV */
200: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
202: /* Evaluate Lagrangian function and gradient */
203: /* p0 */
204: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
205: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
206: if (tao->jacobian_state_pre) {
207: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
208: } else {
209: pset = pflag = PETSC_TRUE;
210: }
211: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
212: else symmetric = PETSC_FALSE;
214: lclP->solve_type = LCL_ADJOINT2;
215: if (tao->jacobian_state_inv) {
216: if (symmetric) {
217: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
218: } else {
219: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
220: }
221: } else {
222: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
223: if (symmetric) {
224: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
225: } else {
226: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
227: }
228: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
229: tao->ksp_its += its;
230: tao->ksp_tot_its += its;
231: }
232: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
233: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
235: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
236: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
238: /* Evaluate constraint norm */
239: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
240: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
242: /* Monitor convergence */
243: tao->reason = TAO_CONTINUE_ITERATING;
244: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
245: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
246: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
248: while (tao->reason == TAO_CONTINUE_ITERATING) {
249: /* Call general purpose update function */
250: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
251: tao->ksp_its = 0;
252: /* Compute a descent direction for the linearly constrained subproblem
253: minimize f(u+du, v+dv)
254: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
256: /* Store the points around the linearization */
257: PetscCall(VecCopy(lclP->U, lclP->U0));
258: PetscCall(VecCopy(lclP->V, lclP->V0));
259: PetscCall(VecCopy(lclP->GU, lclP->GU0));
260: PetscCall(VecCopy(lclP->GV, lclP->GV0));
261: PetscCall(VecCopy(lclP->GAugL_U, lclP->GAugL_U0));
262: PetscCall(VecCopy(lclP->GAugL_V, lclP->GAugL_V0));
263: PetscCall(VecCopy(lclP->GL_U, lclP->GL_U0));
264: PetscCall(VecCopy(lclP->GL_V, lclP->GL_V0));
266: lclP->aug0 = lclP->aug;
267: lclP->lgn0 = lclP->lgn;
269: /* Given the design variables, we need to project the current iterate
270: onto the linearized constraint. We choose to fix the design variables
271: and solve the linear system for the state variables. The resulting
272: point is the Newton direction */
274: /* Solve r = A\con */
275: lclP->solve_type = LCL_FORWARD1;
276: PetscCall(VecSet(lclP->r, 0.0)); /* Initial guess in CG */
278: if (tao->jacobian_state_inv) {
279: PetscCall(MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r));
280: } else {
281: PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
282: PetscCall(KSPSolve(tao->ksp, tao->constraints, lclP->r));
283: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
284: tao->ksp_its += its;
285: tao->ksp_tot_its += tao->ksp_its;
286: }
288: /* Set design step direction dv to zero */
289: PetscCall(VecSet(lclP->s, 0.0));
291: /*
292: Check sufficient descent for constraint merit function .5*||con||^2
293: con' Ak r >= eps1 ||r||^(2+eps2)
294: */
296: /* Compute WU= Ak' * con */
297: if (symmetric) {
298: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->WU));
299: } else {
300: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->WU));
301: }
302: /* Compute r * Ak' * con */
303: PetscCall(VecDot(lclP->r, lclP->WU, &rWU));
305: /* compute ||r||^(2+eps2) */
306: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
307: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
308: adec = lclP->eps1 * r2;
310: if (rWU < adec) {
311: PetscCall(PetscInfo(tao, "Newton direction not descent for constraint, feasibility phase required\n"));
312: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent));
314: PetscCall(PetscInfo(tao, "Using steepest descent direction instead.\n"));
315: PetscCall(VecSet(lclP->r, 0.0));
316: PetscCall(VecAXPY(lclP->r, -1.0, lclP->WU));
317: PetscCall(VecDot(lclP->r, lclP->r, &rWU));
318: PetscCall(VecNorm(lclP->r, NORM_2, &r2));
319: r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
320: PetscCall(VecDot(lclP->r, lclP->GAugL_U, &descent));
321: adec = lclP->eps1 * r2;
322: }
324: /*
325: Check descent for aug. lagrangian
326: r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
327: GL_U = GUk - Ak'*yk
328: WU = Ak'*con
329: adec=eps1||r||^(2+eps2)
331: ==>
332: Check r'GL_U - rho*r'WU <= adec
333: */
335: PetscCall(VecDot(lclP->r, lclP->GL_U, &rGL_U));
336: aldescent = rGL_U - lclP->rho * rWU;
337: if (aldescent > -adec) {
338: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent));
339: PetscCall(PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent));
340: lclP->rho = (rGL_U - adec) / rWU;
341: 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);
342: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Increasing penalty parameter to %g\n", (double)lclP->rho));
343: PetscCall(PetscInfo(tao, " Increasing penalty parameter to %g\n", (double)lclP->rho));
344: }
346: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
347: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
348: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
350: /* We now minimize the augmented Lagrangian along the Newton direction */
351: PetscCall(VecScale(lclP->r, -1.0));
352: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
353: PetscCall(VecScale(lclP->r, -1.0));
354: PetscCall(LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL));
355: PetscCall(LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0));
357: lclP->recompute_jacobian_flag = PETSC_TRUE;
359: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
360: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
361: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
362: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
363: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step));
365: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
366: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
367: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
369: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
371: /* Check convergence */
372: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
373: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
374: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
375: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
376: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
377: if (tao->reason != TAO_CONTINUE_ITERATING) break;
379: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
380: for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
381: /* We now minimize the objective function starting from the fraction of
382: the Newton point accepted by applying one step of a reduced-space
383: method. The optimization problem is:
385: minimize f(u+du, v+dv)
386: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
388: In particular, we have that
389: du = -inv(A)*(Bdv + alpha g) */
391: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
392: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
394: /* Store V and constraints */
395: PetscCall(VecCopy(lclP->V, lclP->V1));
396: PetscCall(VecCopy(tao->constraints, lclP->con1));
398: /* Compute multipliers */
399: /* p1 */
400: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
401: lclP->solve_type = LCL_ADJOINT1;
402: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
403: if (tao->jacobian_state_pre) {
404: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
405: } else {
406: pset = pflag = PETSC_TRUE;
407: }
408: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
409: else symmetric = PETSC_FALSE;
411: if (tao->jacobian_state_inv) {
412: if (symmetric) {
413: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
414: } else {
415: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
416: }
417: } else {
418: if (symmetric) {
419: PetscCall(KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda));
420: } else {
421: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda));
422: }
423: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
424: tao->ksp_its += its;
425: tao->ksp_tot_its += its;
426: }
427: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1));
428: PetscCall(VecAXPY(lclP->g1, -1.0, lclP->GAugL_V));
430: /* Compute the limited-memory quasi-newton direction */
431: if (tao->niter > 0) {
432: PetscCall(MatSolve(lclP->R, lclP->g1, lclP->s));
433: PetscCall(VecDot(lclP->s, lclP->g1, &descent));
434: if (descent <= 0) {
435: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent));
436: PetscCall(VecCopy(lclP->g1, lclP->s));
437: }
438: } else {
439: PetscCall(VecCopy(lclP->g1, lclP->s));
440: }
441: PetscCall(VecScale(lclP->g1, -1.0));
443: /* Recover the full space direction */
444: PetscCall(MatMult(tao->jacobian_design, lclP->s, lclP->WU));
445: /* PetscCall(VecSet(lclP->r,0.0)); */ /* Initial guess in CG */
446: lclP->solve_type = LCL_FORWARD2;
447: if (tao->jacobian_state_inv) {
448: PetscCall(MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r));
449: } else {
450: PetscCall(KSPSolve(tao->ksp, lclP->WU, lclP->r));
451: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
452: tao->ksp_its += its;
453: tao->ksp_tot_its += its;
454: }
456: /* We now minimize the augmented Lagrangian along the direction -r,s */
457: PetscCall(VecScale(lclP->r, -1.0));
458: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
459: PetscCall(VecScale(lclP->r, -1.0));
460: lclP->recompute_jacobian_flag = PETSC_TRUE;
462: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
463: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
464: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
465: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
466: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength = %g\n", (double)step));
468: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
469: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
470: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
471: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
472: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
474: /* Compute the reduced gradient at the new point */
476: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
477: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
479: /* p2 */
480: /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
481: if (phase2_iter == 0) {
482: PetscCall(VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0));
483: } else {
484: PetscCall(VecAXPY(lclP->lambda, -lclP->rho, tao->constraints));
485: }
487: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
488: if (tao->jacobian_state_pre) {
489: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
490: } else {
491: pset = pflag = PETSC_TRUE;
492: }
493: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
494: else symmetric = PETSC_FALSE;
496: lclP->solve_type = LCL_ADJOINT2;
497: if (tao->jacobian_state_inv) {
498: if (symmetric) {
499: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
500: } else {
501: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
502: }
503: } else {
504: if (symmetric) {
505: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
506: } else {
507: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
508: }
509: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
510: tao->ksp_its += its;
511: tao->ksp_tot_its += its;
512: }
514: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2));
515: PetscCall(VecAXPY(lclP->g2, -1.0, lclP->GV));
517: PetscCall(VecScale(lclP->g2, -1.0));
519: /* Update the quasi-newton approximation */
520: PetscCall(MatLMVMUpdate(lclP->R, lclP->V, lclP->g2));
521: /* 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 */
522: }
524: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
526: /* Evaluate Function, Gradient, Constraints, and Jacobian */
527: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
528: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
529: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
531: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
532: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
533: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
535: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
537: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
539: /* Evaluate constraint norm */
540: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
542: /* Monitor convergence */
543: tao->niter++;
544: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
545: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
546: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
547: }
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /*MC
552: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
554: + -tao_lcl_eps1 - epsilon 1 tolerance
555: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
556: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
557: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
558: . -tao_lcl_phase2_niter - Number of phase 2 iterations in 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
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;
581: /* Override default settings (unless already changed) */
582: PetscCall(TaoParametersInitialize(tao));
583: PetscObjectParameterSetDefault(tao, max_it, 200);
584: PetscObjectParameterSetDefault(tao, catol, 1.0e-4);
585: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
586: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
587: PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
589: lclP->rho0 = 1.0e-4;
590: lclP->rhomax = 1e5;
591: lclP->eps1 = 1.0e-8;
592: lclP->eps2 = 0.0;
593: lclP->solve_type = 2;
594: lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
595: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
596: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
597: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
598: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
600: PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao));
601: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
602: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
603: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
604: PetscCall(KSPSetFromOptions(tao->ksp));
606: PetscCall(MatCreate(((PetscObject)tao)->comm, &lclP->R));
607: PetscCall(MatSetType(lclP->R, MATLMVMBFGS));
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
612: {
613: Tao tao = (Tao)ptr;
614: TAO_LCL *lclP = (TAO_LCL *)tao->data;
615: PetscBool set, pset, flag, pflag, symmetric;
616: PetscReal cdotl;
618: PetscFunctionBegin;
619: PetscCall(TaoComputeObjectiveAndGradient(tao, X, f, G));
620: PetscCall(LCLScatter(lclP, G, lclP->GU, lclP->GV));
621: if (lclP->recompute_jacobian_flag) {
622: PetscCall(TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
623: PetscCall(TaoComputeJacobianDesign(tao, X, tao->jacobian_design));
624: }
625: PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
626: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
627: if (tao->jacobian_state_pre) {
628: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
629: } else {
630: pset = pflag = PETSC_TRUE;
631: }
632: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
633: else symmetric = PETSC_FALSE;
635: PetscCall(VecDot(lclP->lambda0, tao->constraints, &cdotl));
636: lclP->lgn = *f - cdotl;
638: /* Gradient of Lagrangian GL = G - J' * lambda */
639: /* WU = A' * WL
640: WV = B' * WL */
641: if (symmetric) {
642: PetscCall(MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
643: } else {
644: PetscCall(MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
645: }
646: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V));
647: PetscCall(VecScale(lclP->GL_U, -1.0));
648: PetscCall(VecScale(lclP->GL_V, -1.0));
649: PetscCall(VecAXPY(lclP->GL_U, 1.0, lclP->GU));
650: PetscCall(VecAXPY(lclP->GL_V, 1.0, lclP->GV));
651: PetscCall(LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL));
653: f[0] = lclP->lgn;
654: PetscCall(VecCopy(lclP->GL, G));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
659: {
660: Tao tao = (Tao)ptr;
661: TAO_LCL *lclP = (TAO_LCL *)tao->data;
662: PetscReal con2;
663: PetscBool flag, pflag, set, pset, symmetric;
665: PetscFunctionBegin;
666: PetscCall(LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao));
667: PetscCall(LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V));
668: PetscCall(VecDot(tao->constraints, tao->constraints, &con2));
669: lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;
671: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
672: /* WU = A' * c
673: WV = B' * c */
674: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
675: if (tao->jacobian_state_pre) {
676: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
677: } else {
678: pset = pflag = PETSC_TRUE;
679: }
680: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
681: else symmetric = PETSC_FALSE;
683: if (symmetric) {
684: PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
685: } else {
686: PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
687: }
689: PetscCall(MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V));
690: PetscCall(VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U));
691: PetscCall(VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V));
692: PetscCall(LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL));
694: f[0] = lclP->aug;
695: PetscCall(VecCopy(lclP->GAugL, G));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: static PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
700: {
701: PetscFunctionBegin;
702: PetscCall(VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
703: PetscCall(VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
704: PetscCall(VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
705: PetscCall(VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
708: static PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
709: {
710: PetscFunctionBegin;
711: PetscCall(VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
712: PetscCall(VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
713: PetscCall(VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
714: PetscCall(VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }