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: if (lclP->rho > lclP->rhomax) {
342: lclP->rho = lclP->rhomax;
343: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "rho=%g > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)", (double)lclP->rho);
344: }
345: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Increasing penalty parameter to %g\n", (double)lclP->rho));
346: PetscCall(PetscInfo(tao, " Increasing penalty parameter to %g\n", (double)lclP->rho));
347: }
349: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
350: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
351: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
353: /* We now minimize the augmented Lagrangian along the Newton direction */
354: PetscCall(VecScale(lclP->r, -1.0));
355: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
356: PetscCall(VecScale(lclP->r, -1.0));
357: PetscCall(LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL));
358: PetscCall(LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0));
360: lclP->recompute_jacobian_flag = PETSC_TRUE;
362: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
363: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
364: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
365: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
366: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step));
368: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
369: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
370: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
372: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
374: /* Check convergence */
375: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
376: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
377: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
378: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
379: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
380: if (tao->reason != TAO_CONTINUE_ITERATING) break;
382: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
383: for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
384: /* We now minimize the objective function starting from the fraction of
385: the Newton point accepted by applying one step of a reduced-space
386: method. The optimization problem is:
388: minimize f(u+du, v+dv)
389: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
391: In particular, we have that
392: du = -inv(A)*(Bdv + alpha g) */
394: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
395: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
397: /* Store V and constraints */
398: PetscCall(VecCopy(lclP->V, lclP->V1));
399: PetscCall(VecCopy(tao->constraints, lclP->con1));
401: /* Compute multipliers */
402: /* p1 */
403: PetscCall(VecSet(lclP->lambda, 0.0)); /* Initial guess in CG */
404: lclP->solve_type = LCL_ADJOINT1;
405: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
406: if (tao->jacobian_state_pre) {
407: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
408: } else {
409: pset = pflag = PETSC_TRUE;
410: }
411: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
412: else symmetric = PETSC_FALSE;
414: if (tao->jacobian_state_inv) {
415: if (symmetric) {
416: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
417: } else {
418: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
419: }
420: } else {
421: if (symmetric) {
422: PetscCall(KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda));
423: } else {
424: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda));
425: }
426: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
427: tao->ksp_its += its;
428: tao->ksp_tot_its += its;
429: }
430: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1));
431: PetscCall(VecAXPY(lclP->g1, -1.0, lclP->GAugL_V));
433: /* Compute the limited-memory quasi-newton direction */
434: if (tao->niter > 0) {
435: PetscCall(MatSolve(lclP->R, lclP->g1, lclP->s));
436: PetscCall(VecDot(lclP->s, lclP->g1, &descent));
437: if (descent <= 0) {
438: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent));
439: PetscCall(VecCopy(lclP->g1, lclP->s));
440: }
441: } else {
442: PetscCall(VecCopy(lclP->g1, lclP->s));
443: }
444: PetscCall(VecScale(lclP->g1, -1.0));
446: /* Recover the full space direction */
447: PetscCall(MatMult(tao->jacobian_design, lclP->s, lclP->WU));
448: /* PetscCall(VecSet(lclP->r,0.0)); */ /* Initial guess in CG */
449: lclP->solve_type = LCL_FORWARD2;
450: if (tao->jacobian_state_inv) {
451: PetscCall(MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r));
452: } else {
453: PetscCall(KSPSolve(tao->ksp, lclP->WU, lclP->r));
454: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
455: tao->ksp_its += its;
456: tao->ksp_tot_its += its;
457: }
459: /* We now minimize the augmented Lagrangian along the direction -r,s */
460: PetscCall(VecScale(lclP->r, -1.0));
461: PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
462: PetscCall(VecScale(lclP->r, -1.0));
463: lclP->recompute_jacobian_flag = PETSC_TRUE;
465: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
466: PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
467: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
468: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
469: if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength = %g\n", (double)step));
471: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
472: PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
473: PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
474: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
475: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
477: /* Compute the reduced gradient at the new point */
479: PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
480: PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));
482: /* p2 */
483: /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
484: if (phase2_iter == 0) {
485: PetscCall(VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0));
486: } else {
487: PetscCall(VecAXPY(lclP->lambda, -lclP->rho, tao->constraints));
488: }
490: PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
491: if (tao->jacobian_state_pre) {
492: PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
493: } else {
494: pset = pflag = PETSC_TRUE;
495: }
496: if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
497: else symmetric = PETSC_FALSE;
499: lclP->solve_type = LCL_ADJOINT2;
500: if (tao->jacobian_state_inv) {
501: if (symmetric) {
502: PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
503: } else {
504: PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
505: }
506: } else {
507: if (symmetric) {
508: PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
509: } else {
510: PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
511: }
512: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
513: tao->ksp_its += its;
514: tao->ksp_tot_its += its;
515: }
517: PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2));
518: PetscCall(VecAXPY(lclP->g2, -1.0, lclP->GV));
520: PetscCall(VecScale(lclP->g2, -1.0));
522: /* Update the quasi-newton approximation */
523: PetscCall(MatLMVMUpdate(lclP->R, lclP->V, lclP->g2));
524: /* 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 */
525: }
527: PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
529: /* Evaluate Function, Gradient, Constraints, and Jacobian */
530: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
531: PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
532: PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));
534: PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
535: PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
536: PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));
538: PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
540: PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
542: /* Evaluate constraint norm */
543: PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
545: /* Monitor convergence */
546: tao->niter++;
547: PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
548: PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
549: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
550: }
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: /*MC
555: TAOLCL - linearly constrained lagrangian method for pde-constrained optimization
557: + -tao_lcl_eps1 - epsilon 1 tolerance
558: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
559: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
560: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
561: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
562: . -tao_lcl_verbose - Print verbose output if True
563: . -tao_lcl_tola - Tolerance for first forward solve
564: . -tao_lcl_tolb - Tolerance for first adjoint solve
565: . -tao_lcl_tolc - Tolerance for second forward solve
566: - -tao_lcl_told - Tolerance for second adjoint solve
568: Level: beginner
569: M*/
570: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
571: {
572: TAO_LCL *lclP;
573: const char *morethuente_type = TAOLINESEARCHMT;
575: PetscFunctionBegin;
576: tao->ops->setup = TaoSetup_LCL;
577: tao->ops->solve = TaoSolve_LCL;
578: tao->ops->view = TaoView_LCL;
579: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
580: tao->ops->destroy = TaoDestroy_LCL;
581: PetscCall(PetscNew(&lclP));
582: tao->data = (void *)lclP;
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: }