Actual source code: lcl.c
1: #include "lcl.h"
2: #include "src/matrix/lmvmmat.h"
3: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
4: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*);
5: static PetscErrorCode LCLScatter(TAO_LCL*,Vec,Vec,Vec);
6: static PetscErrorCode LCLGather(TAO_LCL*,Vec,Vec,Vec);
10: static PetscErrorCode TaoDestroy_LCL(TaoSolver tao)
11: {
12: TAO_LCL *lclP = (TAO_LCL*)tao->data;
15: if (tao->setupcalled) {
16: MatDestroy(&lclP->R);
17:
18: VecDestroy(&lclP->lamda);
19: VecDestroy(&lclP->lamda0);
20: VecDestroy(&lclP->WL);
21: VecDestroy(&lclP->W);
22: VecDestroy(&lclP->X0);
23: VecDestroy(&lclP->G0);
24: VecDestroy(&lclP->GL);
25: VecDestroy(&lclP->GAugL);
26: VecDestroy(&lclP->dbar);
27: VecDestroy(&lclP->U);
28: VecDestroy(&lclP->U0);
29: VecDestroy(&lclP->V);
30: VecDestroy(&lclP->V0);
31: VecDestroy(&lclP->V1);
32: VecDestroy(&lclP->GU);
33: VecDestroy(&lclP->GV);
34: VecDestroy(&lclP->GU0);
35: VecDestroy(&lclP->GV0);
36: VecDestroy(&lclP->GL_U);
37: VecDestroy(&lclP->GL_V);
38: VecDestroy(&lclP->GAugL_U);
39: VecDestroy(&lclP->GAugL_V);
40: VecDestroy(&lclP->GL_U0);
41: VecDestroy(&lclP->GL_V0);
42: VecDestroy(&lclP->GAugL_U0);
43: VecDestroy(&lclP->GAugL_V0);
44: VecDestroy(&lclP->DU);
45: VecDestroy(&lclP->DV);
46: VecDestroy(&lclP->WU);
47: VecDestroy(&lclP->WV);
48: VecDestroy(&lclP->g1);
49: VecDestroy(&lclP->g2);
50: VecDestroy(&lclP->con1);
52: VecDestroy(&lclP->r);
53: VecDestroy(&lclP->s);
55: ISDestroy(&tao->state_is);
56: ISDestroy(&tao->design_is);
58: VecScatterDestroy(&lclP->state_scatter);
59: VecScatterDestroy(&lclP->design_scatter);
60: }
61: PetscFree(tao->data);
62: tao->data = PETSC_NULL;
64: return(0);
65: }
69: static PetscErrorCode TaoSetFromOptions_LCL(TaoSolver tao)
70: {
71: /*TAO_LCL *lclP = (TAO_LCL*)tao->data;*/
74: TaoLineSearchSetFromOptions(tao->linesearch);
75: return(0);
76: }
80: static PetscErrorCode TaoView_LCL(TaoSolver tao, PetscViewer viewer)
81: {
82: /*
83: TAO_LCL *lclP = (TAO_LCL*)tao->data;
86: return(0);
87: */
88: return 0;
89: }
93: static PetscErrorCode TaoSetup_LCL(TaoSolver tao)
94: {
95: TAO_LCL *lclP = (TAO_LCL*)tao->data;
96: PetscInt lo, hi, nlocalstate, nlocaldesign;
98: PetscBool flag;
99: IS is_state, is_design;
101: /* Check for state/design IS */
102: if (!tao->state_is) {
103: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()");
104: }
105: VecDuplicate(tao->solution, &tao->gradient);
106: VecDuplicate(tao->solution, &tao->stepdirection);
107: VecDuplicate(tao->solution, &lclP->W);
108: VecDuplicate(tao->solution, &lclP->X0);
109: VecDuplicate(tao->solution, &lclP->G0);
110: VecDuplicate(tao->solution, &lclP->GL);
111: VecDuplicate(tao->solution, &lclP->GAugL);
112:
113: VecDuplicate(tao->constraints, &lclP->lamda);
114: VecDuplicate(tao->constraints, &lclP->WL);
115: VecDuplicate(tao->constraints, &lclP->lamda0);
116: VecDuplicate(tao->constraints, &lclP->con1);
118: VecSet(lclP->lamda,0.0);
120: VecGetSize(tao->solution, &lclP->n);
121: VecGetSize(tao->constraints, &lclP->m);
124: VecCreate(((PetscObject)tao)->comm,&lclP->U);
125: VecCreate(((PetscObject)tao)->comm,&lclP->V);
126: ISGetLocalSize(tao->state_is,&nlocalstate);
127: ISGetLocalSize(tao->design_is,&nlocaldesign);
128: VecSetSizes(lclP->U,nlocalstate,lclP->m);
129: VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m);
130: VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name);
131: VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name);
132: VecSetFromOptions(lclP->U);
133: VecSetFromOptions(lclP->V);
134: VecDuplicate(lclP->U,&lclP->DU);
135: VecDuplicate(lclP->U,&lclP->U0);
136: VecDuplicate(lclP->U,&lclP->GU);
137: VecDuplicate(lclP->U,&lclP->GU0);
138: VecDuplicate(lclP->U,&lclP->GAugL_U);
139: VecDuplicate(lclP->U,&lclP->GL_U);
140: VecDuplicate(lclP->U,&lclP->GAugL_U0);
141: VecDuplicate(lclP->U,&lclP->GL_U0);
142: VecDuplicate(lclP->U,&lclP->WU);
143: VecDuplicate(lclP->U,&lclP->r);
144: VecDuplicate(lclP->V,&lclP->V0);
145: VecDuplicate(lclP->V,&lclP->V1);
146: VecDuplicate(lclP->V,&lclP->DV);
147: VecDuplicate(lclP->V,&lclP->s);
148: VecDuplicate(lclP->V,&lclP->GV);
149: VecDuplicate(lclP->V,&lclP->GV0);
150: VecDuplicate(lclP->V,&lclP->dbar);
151: VecDuplicate(lclP->V,&lclP->GAugL_V);
152: VecDuplicate(lclP->V,&lclP->GL_V);
153: VecDuplicate(lclP->V,&lclP->GAugL_V0);
154: VecDuplicate(lclP->V,&lclP->GL_V0);
155: VecDuplicate(lclP->V,&lclP->WV);
156: VecDuplicate(lclP->V,&lclP->g1);
157: VecDuplicate(lclP->V,&lclP->g2);
158:
160:
162: /* create scatters for state, design subvecs */
163: VecGetOwnershipRange(lclP->U,&lo,&hi);
164: ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state);
165: VecGetOwnershipRange(lclP->V,&lo,&hi);
166: if (0) {
167: PetscInt sizeU,sizeV;
168: VecGetSize(lclP->U,&sizeU);
169: VecGetSize(lclP->V,&sizeV);
170: PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV);
171: }
172: ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design);
173: VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter);
174: VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter);
175: ISDestroy(&is_state);
176: ISDestroy(&is_design);
179: lclP->phase2_niter = 1;
180: PetscOptionsInt("-tao_lcl_phase2_niter","Number of phase 2 iterations in LCL algorithm","",lclP->phase2_niter,&lclP->phase2_niter,&flag);
181: lclP->verbose = PETSC_FALSE;
182: PetscOptionsBool("-tao_lcl_verbose","Print verbose output","",lclP->verbose,&lclP->verbose,&flag);
183: lclP->tola = 1e-4;
184: PetscOptionsReal("-tao_lcl_tola","Tolerance for first forward solve","",lclP->tola,&lclP->tola,&flag);
185: lclP->tolb = 1e-4;
186: PetscOptionsReal("-tao_lcl_tolb","Tolerance for first adjoint solve","",lclP->tolb,&lclP->tolb,&flag);
187: lclP->tolc = 1e-4;
188: PetscOptionsReal("-tao_lcl_tolc","Tolerance for second forward solve","",lclP->tolc,&lclP->tolc,&flag);
189: lclP->told = 1e-4;
190: PetscOptionsReal("-tao_lcl_told","Tolerance for second adjoint solve","",lclP->told,&lclP->told,&flag);
193: return(0);
194: }
198: static PetscErrorCode TaoSolve_LCL(TaoSolver tao)
199: {
200: TAO_LCL *lclP = (TAO_LCL*)tao->data;
201: PetscInt iter=0,phase2_iter,nlocal,its;
202: TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
203: TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
204: PetscReal step=1.0,f, descent,con2;
205: PetscReal cnorm, mnorm;
206: PetscReal adec,r2,rGL_U,rWU;
207: PetscBool set,pset,flag,pflag,symmetric;
211: VecGetLocalSize(lclP->U,&nlocal);
212: VecGetLocalSize(lclP->V,&nlocal);
213: MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R);
214: MatLMVMAllocateVectors(lclP->R,lclP->V);
215: lclP->rho = 1.0e-4;
216: lclP->recompute_jacobian_flag = PETSC_TRUE;
217:
218: /* Scatter to U,V */
219: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
220:
221: /* Evaluate Function, Gradient, Constraints, and Jacobian */
222: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
223: TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag);
224: TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design);
225: TaoComputeConstraints(tao,tao->solution, tao->constraints);
227:
228: /* Scatter gradient to GU,GV */
229: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
230:
232: /* Evaluate Lagrangian function and gradient */
233: /* p0 */
234: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
235: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
236: if (tao->jacobian_state_pre) {
237: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
238: } else {
239: pset = pflag = PETSC_TRUE;
240: }
241: if (set && pset && flag && pflag)
242: symmetric = PETSC_TRUE;
243: else
244: symmetric = PETSC_FALSE;
246: lclP->solve_type = LCL_ADJOINT2;
247: if (tao->jacobian_state_inv) {
248: if (symmetric) {
249: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); } else {
250: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
251: }
252: } else {
253: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag);
254: if (symmetric) {
255: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
256: } else {
257: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
258: }
259: KSPGetIterationNumber(tao->ksp,&its);
260: tao->ksp_its += its;
261: }
262:
263: VecCopy(lclP->lamda,lclP->lamda0);
265: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
267: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
268: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
269:
271: /* Evaluate constraint norm */
272: VecNorm(tao->constraints, NORM_2, &cnorm);
273: VecNorm(lclP->GAugL, NORM_2, &mnorm);
274:
276: /* Monitor convergence */
277: TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);
279: while (reason == TAO_CONTINUE_ITERATING) {
280: /* Compute a descent direction for the linearly constrained subproblem
281: minimize f(u+du, v+dv)
282: s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */
284: /* Store the points around the linearization */
285: VecCopy(lclP->U, lclP->U0);
286: VecCopy(lclP->V, lclP->V0);
287: VecCopy(lclP->GU,lclP->GU0);
288: VecCopy(lclP->GV,lclP->GV0);
289: VecCopy(lclP->GAugL_U,lclP->GAugL_U0);
290: VecCopy(lclP->GAugL_V,lclP->GAugL_V0);
291: VecCopy(lclP->GL_U,lclP->GL_U0);
292: VecCopy(lclP->GL_V,lclP->GL_V0);
294: lclP->aug0 = lclP->aug;
295: lclP->lgn0 = lclP->lgn;
297: /* Given the design variables, we need to project the current iterate
298: onto the linearized constraint. We choose to fix the design variables
299: and solve the linear system for the state variables. The resulting
300: point is the Newton direction */
301:
302: /* Solve r = A\con */
303: lclP->solve_type = LCL_FORWARD1;
304: VecSet(lclP->r,0.0); /* Initial guess in CG */
305:
306: if (tao->jacobian_state_inv) {
307: MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r);
308: } else {
309: KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag);
310: KSPSolve(tao->ksp, tao->constraints, lclP->r);
311: KSPGetIterationNumber(tao->ksp,&its);
312: tao->ksp_its += its;
313: }
314:
316: VecNorm(lclP->r,NORM_2,&cnorm);
317: VecSet(lclP->s, 0.0);
319: /* Make sure the Newton direction is a descent direction for the merit function */
320: if (symmetric) {
321: MatMult(tao->jacobian_state,tao->constraints,lclP->WU);
322: } else {
323: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU);
324: }
325:
326: VecDot(lclP->r,lclP->WU,&descent);
327: if (descent <= 0) {
328: PetscInfo1(tao,"Newton direction not descent: %G\n",descent);
329: if (lclP->verbose) {
330: PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent: %G\n",descent);
331: }
332: PetscInfo(tao,"Using steepest descent direction instead.\n");
333: VecSet(lclP->r,0.0);
334: VecAXPY(lclP->r,-1.0,lclP->WU);
335: }
337: /* Check descent for aug. lagrangian */
338: VecDot(lclP->r,lclP->r,&r2);
339: VecDot(lclP->r,lclP->GAugL_U,&descent);
340: adec = 1e-8 * r2;
341: if (descent <= adec) {
342: PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %G",descent);
343: VecDot(tao->constraints,tao->constraints,&con2);
344: VecDot(lclP->r,lclP->GL_U,&rGL_U);
345: VecDot(lclP->r,lclP->WU,&rWU);
346: lclP->rho = 2 * (adec - rGL_U) / rWU;
347: PetscInfo1(tao," Increasing penalty parameter to %G",lclP->rho);
348: if (lclP->verbose) {
349: PetscPrintf(PETSC_COMM_WORLD," Increasing penalty parameter to %G\n",lclP->rho);
350: }
351: }
353: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
354: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
355: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
357: /* We now minimize the augmented Lagrangian along the Newton direction */
358: VecScale(lclP->r,-1.0);
359: LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection);
360: VecScale(lclP->r,-1.0);
361: LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL);
362: LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0);
364: lclP->recompute_jacobian_flag = PETSC_TRUE;
366: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
367: TaoLineSearchSetType(tao->linesearch, TAOLINESEARCH_MT);
368: TaoLineSearchSetFromOptions(tao->linesearch);
369: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason);
370: if (lclP->verbose) {
371: PetscPrintf(PETSC_COMM_WORLD,"Steplength = %G\n",step);
372: }
373:
374: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
375: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
376: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
378: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
380: /* Check convergence */
381: VecNorm(lclP->GAugL, NORM_2, &mnorm);
382: VecNorm(tao->constraints, NORM_2, &cnorm);
383: TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);
384: if (reason != TAO_CONTINUE_ITERATING){
385: break;
386: }
389: /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
390: for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++){
391: /* We now minimize the objective function starting from the fraction of
392: the Newton point accepted by applying one step of a reduced-space
393: method. The optimization problem is:
394:
395: minimize f(u+du, v+dv)
396: s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)
398: In particular, we have that
399: du = -inv(A)*(Bdv + alpha g) */
401: TaoComputeJacobianState(tao,lclP->X0,&tao->jacobian_state,&tao->jacobian_state_pre,&tao->jacobian_state_inv,&lclP->statematflag);
402: TaoComputeJacobianDesign(tao,lclP->X0,&tao->jacobian_design);
404: /* Store V and constraints */
405: VecCopy(lclP->V, lclP->V1);
406: VecCopy(tao->constraints,lclP->con1);
408: /* Compute multipliers */
409: /* p1 */
410: VecSet(lclP->lamda,0.0); /* Initial guess in CG */
411: lclP->solve_type = LCL_ADJOINT1;
412: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
413: if (tao->jacobian_state_pre) {
414: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
415: } else {
416: pset = pflag = PETSC_TRUE;
417: }
418: if (set && pset && flag && pflag)
419: symmetric = PETSC_TRUE;
420: else
421: symmetric = PETSC_FALSE;
422:
423: if (tao->jacobian_state_inv) {
424: if (symmetric) {
425: MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
426: } else {
427: MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda);
428: }
429: } else {
430: if (symmetric) {
431: KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lamda);
432: } else {
433: KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lamda);
434: }
435: KSPGetIterationNumber(tao->ksp,&its);
436: tao->ksp_its += its;
437: }
438: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1);
439: VecAXPY(lclP->g1,-1.0,lclP->GAugL_V);
440:
441: /* Compute the limited-memory quasi-newton direction */
442: if (iter > 0) {
443: MatLMVMSolve(lclP->R,lclP->g1,lclP->s);
444: VecDot(lclP->s,lclP->g1,&descent);
445: if (descent < 0e-8) {
446: if (lclP->verbose) {
447: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %G\n",descent);
448: }
449: VecCopy(lclP->g1,lclP->s);
450: }
451: } else {
452: VecCopy(lclP->g1,lclP->s);
453: }
454: VecScale(lclP->g1,-1.0);
457: /* Recover the full space direction */
458: MatMult(tao->jacobian_design,lclP->s,lclP->WU);
459: //VecSet(lclP->r,0.0); /* Initial guess in CG */
460: lclP->solve_type = LCL_FORWARD2;
461: if (tao->jacobian_state_inv) {
462: MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r);
463: } else {
464: KSPSolve(tao->ksp, lclP->WU, lclP->r);
465: KSPGetIterationNumber(tao->ksp,&its);
466: tao->ksp_its += its;
467: }
469: /* We now minimize the augmented Lagrangian along the direction -r,s */
470: VecScale(lclP->r, -1.0);
471: LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection);
472: VecScale(lclP->r, -1.0);
473: lclP->recompute_jacobian_flag = PETSC_TRUE;
475: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
476: TaoLineSearchSetType(tao->linesearch,TAOLINESEARCH_MT);
477: TaoLineSearchSetFromOptions(tao->linesearch);
478: TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason);
479: if (lclP->verbose){
480: PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength = %G\n",step);
481: }
483: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
484: LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V);
485: LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V);
486: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
487: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
489: /* Compute the reduced gradient at the new point */
491: TaoComputeJacobianState(tao,lclP->X0,&tao->jacobian_state,&tao->jacobian_state_pre,&tao->jacobian_state_inv,&lclP->statematflag);
492: TaoComputeJacobianDesign(tao,lclP->X0,&tao->jacobian_design);
494: /* p2 */
495: /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */
496: if (phase2_iter==0){
497: VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0);
498: }
499: else {
500: VecAXPY(lclP->lamda,-lclP->rho,tao->constraints);
501: }
503: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
504: if (tao->jacobian_state_pre) {
505: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
506: } else {
507: pset = pflag = PETSC_TRUE;
508: }
509: if (set && pset && flag && pflag)
510: symmetric = PETSC_TRUE;
511: else
512: symmetric = PETSC_FALSE;
514: lclP->solve_type = LCL_ADJOINT2;
515: if (tao->jacobian_state_inv) {
516: if (symmetric) {
517: MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
518: } else {
519: MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda);
520: }
521: } else {
522: if (symmetric) {
523: KSPSolve(tao->ksp, lclP->GU, lclP->lamda);
524: } else {
525: KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda);
526: }
527: KSPGetIterationNumber(tao->ksp,&its);
528: tao->ksp_its += its;
529: }
530:
531:
532: MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g2);
533: VecAXPY(lclP->g2,-1.0,lclP->GV);
535: VecScale(lclP->g2,-1.0);
537: /* Update the quasi-newton approximation */
538: if (phase2_iter >= 0){
539: MatLMVMSetPrev(lclP->R,lclP->V1,lclP->g1);
540: }
541: MatLMVMUpdate(lclP->R,lclP->V,lclP->g2);
542: /* 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 */
544: }
546: VecCopy(lclP->lamda,lclP->lamda0);
547:
548: /* Evaluate Function, Gradient, Constraints, and Jacobian */
549: TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient);
550: LCLScatter(lclP,tao->solution,lclP->U,lclP->V);
551: LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV);
553: TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag);
554: TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design);
555: TaoComputeConstraints(tao,tao->solution, tao->constraints);
558: LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao);
560: VecNorm(lclP->GAugL, NORM_2, &mnorm);
562: /* Evaluate constraint norm */
563: VecNorm(tao->constraints, NORM_2, &cnorm);
564:
565: /* Monitor convergence */
566: iter++;
567: TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason);
568:
569: }
571: MatDestroy(&lclP->R);
573: return(0);
574: }
579: PetscErrorCode TaoCreate_LCL(TaoSolver tao)
580: {
581: TAO_LCL *lclP;
583: const char *morethuente_type = TAOLINESEARCH_MT;
585:
586: tao->ops->setup = TaoSetup_LCL;
587: tao->ops->solve = TaoSolve_LCL;
588: tao->ops->view = TaoView_LCL;
589: tao->ops->setfromoptions = TaoSetFromOptions_LCL;
590: tao->ops->destroy = TaoDestroy_LCL;
591:
593: PetscNewLog(tao,TAO_LCL,&lclP);
594: tao->data = (void*)lclP;
597: tao->max_it=200;
598: tao->fatol=1e-8;
599: tao->frtol=1e-8;
600: tao->catol=1e-4;
601: tao->crtol=1e-4;
602: tao->gttol=1e-4;
603: tao->gatol=1e-4;
604: tao->grtol=1e-4;
605:
606: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
607: TaoLineSearchSetType(tao->linesearch, morethuente_type);
609: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao);
610: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
611: KSPSetFromOptions(tao->ksp);
614:
615: return(0);
617: }
623: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
624: {
625: TaoSolver tao = (TaoSolver)ptr;
626: TAO_LCL *lclP = (TAO_LCL*)tao->data;
627: PetscBool set,pset,flag,pflag,symmetric;
628: PetscReal cdotl;
632: TaoComputeObjectiveAndGradient(tao,X,f,G);
633: LCLScatter(lclP,G,lclP->GU,lclP->GV);
634: if (lclP->recompute_jacobian_flag) {
635: TaoComputeJacobianState(tao,X, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag);
636: TaoComputeJacobianDesign(tao,X, &tao->jacobian_design);
637: }
638: TaoComputeConstraints(tao,X, tao->constraints);
639: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
640: if (tao->jacobian_state_pre) {
641: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
642: } else {
643: pset = pflag = PETSC_TRUE;
644: }
645: if (set && pset && flag && pflag)
646: symmetric = PETSC_TRUE;
647: else
648: symmetric = PETSC_FALSE;
650: VecDot(lclP->lamda0, tao->constraints, &cdotl);
651: lclP->lgn = *f - cdotl;
653: /* Gradient of Lagrangian GL = G - J' * lamda */
654: /* WU = A' * WL
655: WV = B' * WL */
656: if (symmetric) {
657: MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
658: } else {
659: MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U);
660: }
661: MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V);
662: VecScale(lclP->GL_U,-1.0);
663: VecScale(lclP->GL_V,-1.0);
664: VecAXPY(lclP->GL_U,1.0,lclP->GU);
665: VecAXPY(lclP->GL_V,1.0,lclP->GV);
666: LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL);
668: f[0] = lclP->lgn;
669: VecCopy(lclP->GL,G);
671: return(0);
672: }
676: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
677: {
678: TaoSolver tao = (TaoSolver)ptr;
679: TAO_LCL *lclP = (TAO_LCL*)tao->data;
680: PetscReal con2;
681: PetscBool flag,pflag,set,pset,symmetric;
685: LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao);
686: LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V);
687: VecDot(tao->constraints,tao->constraints,&con2);
688: lclP->aug = lclP->lgn + 0.5*lclP->rho*con2;
689:
690: /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
691: /* WU = A' * c
692: WV = B' * c */
693: MatIsSymmetricKnown(tao->jacobian_state,&set,&flag);
694: if (tao->jacobian_state_pre) {
695: MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag);
696: } else {
697: pset = pflag = PETSC_TRUE;
698: }
699: if (set && pset && flag && pflag)
700: symmetric = PETSC_TRUE;
701: else
702: symmetric = PETSC_FALSE;
704: if (symmetric) {
705: MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
706: } else {
707: MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U);
708: }
710: MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V);
711: VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U);
712: VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V);
713: LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL);
715: f[0] = lclP->aug;
716: VecCopy(lclP->GAugL,G);
718: return(0);
719: }
723: PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
724: {
727: VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
728: VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE);
729: VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
730: VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE);
731: return(0);
732:
733: }
736: PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
737: {
740: VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
741: VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD);
742: VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
743: VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD);
744: return(0);
745:
746: }