Actual source code: ipm.c
petsc-3.8.4 2018-03-24
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/constrained/impls/ipm/ipm.h>
4: /*
5: x,d in R^n
6: f in R
7: nb = mi + nlb+nub
8: s in R^nb is slack vector CI(x) / x-XL / -x+XU
9: bin in R^mi (tao->constraints_inequality)
10: beq in R^me (tao->constraints_equality)
11: lamdai in R^nb (ipmP->lamdai)
12: lamdae in R^me (ipmP->lamdae)
13: Jeq in R^(me x n) (tao->jacobian_equality)
14: Jin in R^(mi x n) (tao->jacobian_inequality)
15: Ai in R^(nb x n) (ipmP->Ai)
16: H in R^(n x n) (tao->hessian)
17: min f=(1/2)*x'*H*x + d'*x
18: s.t. CE(x) == 0
19: CI(x) >= 0
20: x >= tao->XL
21: -x >= -tao->XU
22: */
24: static PetscErrorCode IPMComputeKKT(Tao tao);
25: static PetscErrorCode IPMPushInitialPoint(Tao tao);
26: static PetscErrorCode IPMEvaluate(Tao tao);
27: static PetscErrorCode IPMUpdateK(Tao tao);
28: static PetscErrorCode IPMUpdateAi(Tao tao);
29: static PetscErrorCode IPMGatherRHS(Tao tao,Vec,Vec,Vec,Vec,Vec);
30: static PetscErrorCode IPMScatterStep(Tao tao,Vec,Vec,Vec,Vec,Vec);
31: static PetscErrorCode IPMInitializeBounds(Tao tao);
33: static PetscErrorCode TaoSolve_IPM(Tao tao)
34: {
35: PetscErrorCode ierr;
36: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
37: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
38: PetscInt its,i;
39: PetscScalar stepsize=1.0;
40: PetscScalar step_s,step_l,alpha,tau,sigma,phi_target;
43: /* Push initial point away from bounds */
44: IPMInitializeBounds(tao);
45: IPMPushInitialPoint(tao);
46: VecCopy(tao->solution,ipmP->rhs_x);
47: IPMEvaluate(tao);
48: IPMComputeKKT(tao);
49: TaoMonitor(tao,tao->niter++,ipmP->kkt_f,ipmP->phi,0.0,1.0,&reason);
51: while (reason == TAO_CONTINUE_ITERATING) {
52: tao->ksp_its=0;
53: IPMUpdateK(tao);
54: /*
55: rhs.x = -rd
56: rhs.lame = -rpe
57: rhs.lami = -rpi
58: rhs.com = -com
59: */
61: VecCopy(ipmP->rd,ipmP->rhs_x);
62: if (ipmP->me > 0) {
63: VecCopy(ipmP->rpe,ipmP->rhs_lamdae);
64: }
65: if (ipmP->nb > 0) {
66: VecCopy(ipmP->rpi,ipmP->rhs_lamdai);
67: VecCopy(ipmP->complementarity,ipmP->rhs_s);
68: }
69: IPMGatherRHS(tao,ipmP->bigrhs,ipmP->rhs_x,ipmP->rhs_lamdae,ipmP->rhs_lamdai,ipmP->rhs_s);
70: VecScale(ipmP->bigrhs,-1.0);
72: /* solve K * step = rhs */
73: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
74: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
76: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
77: KSPGetIterationNumber(tao->ksp,&its);
78: tao->ksp_its += its;
79: tao->ksp_tot_its+=its;
80: /* Find distance along step direction to closest bound */
81: if (ipmP->nb > 0) {
82: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
83: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
84: alpha = PetscMin(step_s,step_l);
85: alpha = PetscMin(alpha,1.0);
86: ipmP->alpha1 = alpha;
87: } else {
88: ipmP->alpha1 = alpha = 1.0;
89: }
91: /* x_aff = x + alpha*d */
92: VecCopy(tao->solution,ipmP->save_x);
93: if (ipmP->me > 0) {
94: VecCopy(ipmP->lamdae,ipmP->save_lamdae);
95: }
96: if (ipmP->nb > 0) {
97: VecCopy(ipmP->lamdai,ipmP->save_lamdai);
98: VecCopy(ipmP->s,ipmP->save_s);
99: }
101: VecAXPY(tao->solution,alpha,tao->stepdirection);
102: if (ipmP->me > 0) {
103: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
104: }
105: if (ipmP->nb > 0) {
106: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
107: VecAXPY(ipmP->s,alpha,ipmP->ds);
108: }
110: /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */
111: if (ipmP->mu == 0.0) {
112: sigma = 0.0;
113: } else {
114: sigma = 1.0/ipmP->mu;
115: }
116: IPMComputeKKT(tao);
117: sigma *= ipmP->mu;
118: sigma*=sigma*sigma;
120: /* revert kkt info */
121: VecCopy(ipmP->save_x,tao->solution);
122: if (ipmP->me > 0) {
123: VecCopy(ipmP->save_lamdae,ipmP->lamdae);
124: }
125: if (ipmP->nb > 0) {
126: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
127: VecCopy(ipmP->save_s,ipmP->s);
128: }
129: IPMComputeKKT(tao);
131: /* update rhs with new complementarity vector */
132: if (ipmP->nb > 0) {
133: VecCopy(ipmP->complementarity,ipmP->rhs_s);
134: VecScale(ipmP->rhs_s,-1.0);
135: VecShift(ipmP->rhs_s,sigma*ipmP->mu);
136: }
137: IPMGatherRHS(tao,ipmP->bigrhs,NULL,NULL,NULL,ipmP->rhs_s);
139: /* solve K * step = rhs */
140: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
141: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
143: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
144: KSPGetIterationNumber(tao->ksp,&its);
145: tao->ksp_its += its;
146: tao->ksp_tot_its+=its;
147: if (ipmP->nb > 0) {
148: /* Get max step size and apply frac-to-boundary */
149: tau = PetscMax(ipmP->taumin,1.0-ipmP->mu);
150: tau = PetscMin(tau,1.0);
151: if (tau != 1.0) {
152: VecScale(ipmP->s,tau);
153: VecScale(ipmP->lamdai,tau);
154: }
155: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
156: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
157: if (tau != 1.0) {
158: VecCopy(ipmP->save_s,ipmP->s);
159: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
160: }
161: alpha = PetscMin(step_s,step_l);
162: alpha = PetscMin(alpha,1.0);
163: } else {
164: alpha = 1.0;
165: }
166: ipmP->alpha2 = alpha;
167: /* TODO make phi_target meaningful */
168: phi_target = ipmP->dec * ipmP->phi;
169: for (i=0; i<11;i++) {
170: VecAXPY(tao->solution,alpha,tao->stepdirection);
171: if (ipmP->nb > 0) {
172: VecAXPY(ipmP->s,alpha,ipmP->ds);
173: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
174: }
175: if (ipmP->me > 0) {
176: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
177: }
179: /* update dual variables */
180: if (ipmP->me > 0) {
181: VecCopy(ipmP->lamdae,tao->DE);
182: }
184: IPMEvaluate(tao);
185: IPMComputeKKT(tao);
186: if (ipmP->phi <= phi_target) break;
187: alpha /= 2.0;
188: }
190: TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,stepsize,&reason);
191: tao->niter++;
192: }
193: return(0);
194: }
196: static PetscErrorCode TaoSetup_IPM(Tao tao)
197: {
198: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
202: ipmP->nb = ipmP->mi = ipmP->me = 0;
203: ipmP->K=0;
204: VecGetSize(tao->solution,&ipmP->n);
205: if (!tao->gradient) {
206: VecDuplicate(tao->solution, &tao->gradient);
207: VecDuplicate(tao->solution, &tao->stepdirection);
208: VecDuplicate(tao->solution, &ipmP->rd);
209: VecDuplicate(tao->solution, &ipmP->rhs_x);
210: VecDuplicate(tao->solution, &ipmP->work);
211: VecDuplicate(tao->solution, &ipmP->save_x);
212: }
213: if (tao->constraints_equality) {
214: VecGetSize(tao->constraints_equality,&ipmP->me);
215: VecDuplicate(tao->constraints_equality,&ipmP->lamdae);
216: VecDuplicate(tao->constraints_equality,&ipmP->dlamdae);
217: VecDuplicate(tao->constraints_equality,&ipmP->rhs_lamdae);
218: VecDuplicate(tao->constraints_equality,&ipmP->save_lamdae);
219: VecDuplicate(tao->constraints_equality,&ipmP->rpe);
220: VecDuplicate(tao->constraints_equality,&tao->DE);
221: }
222: if (tao->constraints_inequality) {
223: VecDuplicate(tao->constraints_inequality,&tao->DI);
224: }
225: return(0);
226: }
228: static PetscErrorCode IPMInitializeBounds(Tao tao)
229: {
230: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
231: Vec xtmp;
232: PetscInt xstart,xend;
233: PetscInt ucstart,ucend; /* user ci */
234: PetscInt ucestart,uceend; /* user ce */
235: PetscInt sstart,send;
236: PetscInt bigsize;
237: PetscInt i,counter,nloc;
238: PetscInt *cind,*xind,*ucind,*uceind,*stepind;
239: VecType vtype;
240: const PetscInt *xli,*xui;
241: PetscInt xl_offset,xu_offset;
242: IS bigxl,bigxu,isuc,isc,isx,sis,is1;
244: MPI_Comm comm;
247: cind=xind=ucind=uceind=stepind=0;
248: ipmP->mi=0;
249: ipmP->nxlb=0;
250: ipmP->nxub=0;
251: ipmP->nb=0;
252: ipmP->nslack=0;
254: VecDuplicate(tao->solution,&xtmp);
255: if (!tao->XL && !tao->XU && tao->ops->computebounds) {
256: TaoComputeVariableBounds(tao);
257: }
258: if (tao->XL) {
259: VecSet(xtmp,PETSC_NINFINITY);
260: VecWhichGreaterThan(tao->XL,xtmp,&ipmP->isxl);
261: ISGetSize(ipmP->isxl,&ipmP->nxlb);
262: } else {
263: ipmP->nxlb=0;
264: }
265: if (tao->XU) {
266: VecSet(xtmp,PETSC_INFINITY);
267: VecWhichLessThan(tao->XU,xtmp,&ipmP->isxu);
268: ISGetSize(ipmP->isxu,&ipmP->nxub);
269: } else {
270: ipmP->nxub=0;
271: }
272: VecDestroy(&xtmp);
273: if (tao->constraints_inequality) {
274: VecGetSize(tao->constraints_inequality,&ipmP->mi);
275: } else {
276: ipmP->mi = 0;
277: }
278: ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;
280: PetscObjectGetComm((PetscObject)tao->solution,&comm);
282: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
283: PetscMalloc1(bigsize,&stepind);
284: PetscMalloc1(ipmP->n,&xind);
285: PetscMalloc1(ipmP->me,&uceind);
286: VecGetOwnershipRange(tao->solution,&xstart,&xend);
288: if (ipmP->nb > 0) {
289: VecCreate(comm,&ipmP->s);
290: VecSetSizes(ipmP->s,PETSC_DECIDE,ipmP->nb);
291: VecSetFromOptions(ipmP->s);
292: VecDuplicate(ipmP->s,&ipmP->ds);
293: VecDuplicate(ipmP->s,&ipmP->rhs_s);
294: VecDuplicate(ipmP->s,&ipmP->complementarity);
295: VecDuplicate(ipmP->s,&ipmP->ci);
297: VecDuplicate(ipmP->s,&ipmP->lamdai);
298: VecDuplicate(ipmP->s,&ipmP->dlamdai);
299: VecDuplicate(ipmP->s,&ipmP->rhs_lamdai);
300: VecDuplicate(ipmP->s,&ipmP->save_lamdai);
302: VecDuplicate(ipmP->s,&ipmP->save_s);
303: VecDuplicate(ipmP->s,&ipmP->rpi);
304: VecDuplicate(ipmP->s,&ipmP->Zero_nb);
305: VecSet(ipmP->Zero_nb,0.0);
306: VecDuplicate(ipmP->s,&ipmP->One_nb);
307: VecSet(ipmP->One_nb,1.0);
308: VecDuplicate(ipmP->s,&ipmP->Inf_nb);
309: VecSet(ipmP->Inf_nb,PETSC_INFINITY);
311: PetscMalloc1(ipmP->nb,&cind);
312: PetscMalloc1(ipmP->mi,&ucind);
313: VecGetOwnershipRange(ipmP->s,&sstart,&send);
315: if (ipmP->mi > 0) {
316: VecGetOwnershipRange(tao->constraints_inequality,&ucstart,&ucend);
317: counter=0;
318: for (i=ucstart;i<ucend;i++) {
319: cind[counter++] = i;
320: }
321: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isuc);
322: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isc);
323: VecScatterCreate(tao->constraints_inequality,isuc,ipmP->ci,isc,&ipmP->ci_scat);
325: ISDestroy(&isuc);
326: ISDestroy(&isc);
327: }
328: /* need to know how may xbound indices are on each process */
329: /* TODO better way */
330: if (ipmP->nxlb) {
331: ISAllGather(ipmP->isxl,&bigxl);
332: ISGetIndices(bigxl,&xli);
333: /* find offsets for this processor */
334: xl_offset = ipmP->mi;
335: for (i=0;i<ipmP->nxlb;i++) {
336: if (xli[i] < xstart) {
337: xl_offset++;
338: } else break;
339: }
340: ISRestoreIndices(bigxl,&xli);
342: ISGetIndices(ipmP->isxl,&xli);
343: ISGetLocalSize(ipmP->isxl,&nloc);
344: for (i=0;i<nloc;i++) {
345: xind[i] = xli[i];
346: cind[i] = xl_offset+i;
347: }
349: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
350: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
351: VecScatterCreate(tao->XL,isx,ipmP->ci,isc,&ipmP->xl_scat);
352: ISDestroy(&isx);
353: ISDestroy(&isc);
354: ISDestroy(&bigxl);
355: }
357: if (ipmP->nxub) {
358: ISAllGather(ipmP->isxu,&bigxu);
359: ISGetIndices(bigxu,&xui);
360: /* find offsets for this processor */
361: xu_offset = ipmP->mi + ipmP->nxlb;
362: for (i=0;i<ipmP->nxub;i++) {
363: if (xui[i] < xstart) {
364: xu_offset++;
365: } else break;
366: }
367: ISRestoreIndices(bigxu,&xui);
369: ISGetIndices(ipmP->isxu,&xui);
370: ISGetLocalSize(ipmP->isxu,&nloc);
371: for (i=0;i<nloc;i++) {
372: xind[i] = xui[i];
373: cind[i] = xu_offset+i;
374: }
376: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
377: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
378: VecScatterCreate(tao->XU,isx,ipmP->ci,isc,&ipmP->xu_scat);
379: ISDestroy(&isx);
380: ISDestroy(&isc);
381: ISDestroy(&bigxu);
382: }
383: }
384: VecCreate(comm,&ipmP->bigrhs);
385: VecGetType(tao->solution,&vtype);
386: VecSetType(ipmP->bigrhs,vtype);
387: VecSetSizes(ipmP->bigrhs,PETSC_DECIDE,bigsize);
388: VecSetFromOptions(ipmP->bigrhs);
389: VecDuplicate(ipmP->bigrhs,&ipmP->bigstep);
391: /* create scatters for step->x and x->rhs */
392: for (i=xstart;i<xend;i++) {
393: stepind[i-xstart] = i;
394: xind[i-xstart] = i;
395: }
396: ISCreateGeneral(comm,xend-xstart,stepind,PETSC_COPY_VALUES,&sis);
397: ISCreateGeneral(comm,xend-xstart,xind,PETSC_COPY_VALUES,&is1);
398: VecScatterCreate(ipmP->bigstep,sis,tao->solution,is1,&ipmP->step1);
399: VecScatterCreate(tao->solution,is1,ipmP->bigrhs,sis,&ipmP->rhs1);
400: ISDestroy(&sis);
401: ISDestroy(&is1);
403: if (ipmP->nb > 0) {
404: for (i=sstart;i<send;i++) {
405: stepind[i-sstart] = i+ipmP->n;
406: cind[i-sstart] = i;
407: }
408: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
409: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
410: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step2);
411: ISDestroy(&sis);
413: for (i=sstart;i<send;i++) {
414: stepind[i-sstart] = i+ipmP->n+ipmP->me;
415: cind[i-sstart] = i;
416: }
417: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
418: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs3);
419: ISDestroy(&sis);
420: ISDestroy(&is1);
421: }
423: if (ipmP->me > 0) {
424: VecGetOwnershipRange(tao->constraints_equality,&ucestart,&uceend);
425: for (i=ucestart;i<uceend;i++) {
426: stepind[i-ucestart] = i + ipmP->n+ipmP->nb;
427: uceind[i-ucestart] = i;
428: }
430: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
431: ISCreateGeneral(comm,uceend-ucestart,uceind,PETSC_COPY_VALUES,&is1);
432: VecScatterCreate(ipmP->bigstep,sis,tao->constraints_equality,is1,&ipmP->step3);
433: ISDestroy(&sis);
435: for (i=ucestart;i<uceend;i++) {
436: stepind[i-ucestart] = i + ipmP->n;
437: }
439: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
440: VecScatterCreate(tao->constraints_equality,is1,ipmP->bigrhs,sis,&ipmP->rhs2);
441: ISDestroy(&sis);
442: ISDestroy(&is1);
443: }
445: if (ipmP->nb > 0) {
446: for (i=sstart;i<send;i++) {
447: stepind[i-sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
448: cind[i-sstart] = i;
449: }
450: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
451: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
452: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step4);
453: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs4);
454: ISDestroy(&sis);
455: ISDestroy(&is1);
456: }
458: PetscFree(stepind);
459: PetscFree(cind);
460: PetscFree(ucind);
461: PetscFree(uceind);
462: PetscFree(xind);
463: return(0);
464: }
466: static PetscErrorCode TaoDestroy_IPM(Tao tao)
467: {
468: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
472: VecDestroy(&ipmP->rd);
473: VecDestroy(&ipmP->rpe);
474: VecDestroy(&ipmP->rpi);
475: VecDestroy(&ipmP->work);
476: VecDestroy(&ipmP->lamdae);
477: VecDestroy(&ipmP->lamdai);
478: VecDestroy(&ipmP->s);
479: VecDestroy(&ipmP->ds);
480: VecDestroy(&ipmP->ci);
482: VecDestroy(&ipmP->rhs_x);
483: VecDestroy(&ipmP->rhs_lamdae);
484: VecDestroy(&ipmP->rhs_lamdai);
485: VecDestroy(&ipmP->rhs_s);
487: VecDestroy(&ipmP->save_x);
488: VecDestroy(&ipmP->save_lamdae);
489: VecDestroy(&ipmP->save_lamdai);
490: VecDestroy(&ipmP->save_s);
492: VecScatterDestroy(&ipmP->step1);
493: VecScatterDestroy(&ipmP->step2);
494: VecScatterDestroy(&ipmP->step3);
495: VecScatterDestroy(&ipmP->step4);
497: VecScatterDestroy(&ipmP->rhs1);
498: VecScatterDestroy(&ipmP->rhs2);
499: VecScatterDestroy(&ipmP->rhs3);
500: VecScatterDestroy(&ipmP->rhs4);
502: VecScatterDestroy(&ipmP->ci_scat);
503: VecScatterDestroy(&ipmP->xl_scat);
504: VecScatterDestroy(&ipmP->xu_scat);
506: VecDestroy(&ipmP->dlamdai);
507: VecDestroy(&ipmP->dlamdae);
508: VecDestroy(&ipmP->Zero_nb);
509: VecDestroy(&ipmP->One_nb);
510: VecDestroy(&ipmP->Inf_nb);
511: VecDestroy(&ipmP->complementarity);
513: VecDestroy(&ipmP->bigrhs);
514: VecDestroy(&ipmP->bigstep);
515: MatDestroy(&ipmP->Ai);
516: MatDestroy(&ipmP->K);
517: ISDestroy(&ipmP->isxu);
518: ISDestroy(&ipmP->isxl);
519: PetscFree(tao->data);
520: return(0);
521: }
523: static PetscErrorCode TaoSetFromOptions_IPM(PetscOptionItems *PetscOptionsObject,Tao tao)
524: {
525: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
529: PetscOptionsHead(PetscOptionsObject,"IPM method for constrained optimization");
530: PetscOptionsBool("-tao_ipm_monitorkkt","monitor kkt status",NULL,ipmP->monitorkkt,&ipmP->monitorkkt,NULL);
531: PetscOptionsReal("-tao_ipm_pushs","parameter to push initial slack variables away from bounds",NULL,ipmP->pushs,&ipmP->pushs,NULL);
532: PetscOptionsReal("-tao_ipm_pushnu","parameter to push initial (inequality) dual variables away from bounds",NULL,ipmP->pushnu,&ipmP->pushnu,NULL);
533: PetscOptionsTail();
534: KSPSetFromOptions(tao->ksp);
535: return(0);
536: }
538: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
539: {
540: return 0;
541: }
543: /* IPMObjectiveAndGradient()
544: f = d'x + 0.5 * x' * H * x
545: rd = H*x + d + Ae'*lame - Ai'*lami
546: rpe = Ae*x - be
547: rpi = Ai*x - yi - bi
548: mu = yi' * lami/mi;
549: com = yi.*lami
551: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
552: */
553: /*
554: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
555: {
556: Tao tao = (Tao)tptr;
557: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
560: IPMComputeKKT(tao);
561: *f = ipmP->phi;
562: return(0);
563: }
564: */
566: /*
567: f = d'x + 0.5 * x' * H * x
568: rd = H*x + d + Ae'*lame - Ai'*lami
569: Ai = jac_ineq
570: I (w/lb)
571: -I (w/ub)
573: rpe = ce
574: rpi = ci - s;
575: com = s.*lami
576: mu = yi' * lami/mi;
578: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
579: */
580: static PetscErrorCode IPMComputeKKT(Tao tao)
581: {
582: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
583: PetscScalar norm;
587: VecCopy(tao->gradient,ipmP->rd);
589: if (ipmP->me > 0) {
590: /* rd = gradient + Ae'*lamdae */
591: MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
592: VecAXPY(ipmP->rd, 1.0, ipmP->work);
594: /* rpe = ce(x) */
595: VecCopy(tao->constraints_equality,ipmP->rpe);
596: }
597: if (ipmP->nb > 0) {
598: /* rd = rd - Ai'*lamdai */
599: MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
600: VecAXPY(ipmP->rd, -1.0, ipmP->work);
602: /* rpi = cin - s */
603: VecCopy(ipmP->ci,ipmP->rpi);
604: VecAXPY(ipmP->rpi, -1.0, ipmP->s);
606: /* com = s .* lami */
607: VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
608: }
609: /* phi = ||rd; rpe; rpi; com|| */
610: VecDot(ipmP->rd,ipmP->rd,&norm);
611: ipmP->phi = norm;
612: if (ipmP->me > 0 ) {
613: VecDot(ipmP->rpe,ipmP->rpe,&norm);
614: ipmP->phi += norm;
615: }
616: if (ipmP->nb > 0) {
617: VecDot(ipmP->rpi,ipmP->rpi,&norm);
618: ipmP->phi += norm;
619: VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
620: ipmP->phi += norm;
621: /* mu = s'*lami/nb */
622: VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
623: ipmP->mu /= ipmP->nb;
624: } else {
625: ipmP->mu = 1.0;
626: }
628: ipmP->phi = PetscSqrtScalar(ipmP->phi);
629: return(0);
630: }
632: /* evaluate user info at current point */
633: PetscErrorCode IPMEvaluate(Tao tao)
634: {
635: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
639: TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
640: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
641: if (ipmP->me > 0) {
642: TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
643: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
644: }
645: if (ipmP->mi > 0) {
646: TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
647: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
648: }
649: if (ipmP->nb > 0) {
650: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
651: IPMUpdateAi(tao);
652: }
653: return(0);
654: }
656: /* Push initial point away from bounds */
657: PetscErrorCode IPMPushInitialPoint(Tao tao)
658: {
659: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
663: TaoComputeVariableBounds(tao);
664: if (tao->XL && tao->XU) {
665: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
666: }
667: if (ipmP->nb > 0) {
668: VecSet(ipmP->s,ipmP->pushs);
669: VecSet(ipmP->lamdai,ipmP->pushnu);
670: if (ipmP->mi > 0) {
671: VecSet(tao->DI,ipmP->pushnu);
672: }
673: }
674: if (ipmP->me > 0) {
675: VecSet(tao->DE,1.0);
676: VecSet(ipmP->lamdae,1.0);
677: }
678: return(0);
679: }
681: PetscErrorCode IPMUpdateAi(Tao tao)
682: {
683: /* Ai = Ji
684: I (w/lb)
685: -I (w/ub) */
687: /* Ci = user->ci
688: Xi - lb (w/lb)
689: -Xi + ub (w/ub) */
691: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
692: MPI_Comm comm;
693: PetscInt i;
694: PetscScalar newval;
695: PetscInt newrow,newcol,ncols;
696: const PetscScalar *vals;
697: const PetscInt *cols;
698: PetscInt astart,aend,jstart,jend;
699: PetscInt *nonzeros;
700: PetscInt r2,r3,r4;
701: PetscMPIInt size;
702: PetscErrorCode ierr;
705: r2 = ipmP->mi;
706: r3 = r2 + ipmP->nxlb;
707: r4 = r3 + ipmP->nxub;
709: if (!ipmP->nb) return(0);
711: /* Create Ai matrix if it doesn't exist yet */
712: if (!ipmP->Ai) {
713: comm = ((PetscObject)(tao->solution))->comm;
714: PetscMalloc1(ipmP->nb,&nonzeros);
715: MPI_Comm_size(comm,&size);
716: if (size == 1) {
717: for (i=0;i<ipmP->mi;i++) {
718: MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
719: nonzeros[i] = ncols;
720: MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
721: }
722: for (i=r2;i<r4;i++) {
723: nonzeros[i] = 1;
724: }
725: }
726: MatCreate(comm,&ipmP->Ai);
727: MatSetType(ipmP->Ai,MATAIJ);
728: MatSetSizes(ipmP->Ai,PETSC_DECIDE,PETSC_DECIDE,ipmP->nb,ipmP->n);
729: MatSetFromOptions(ipmP->Ai);
730: MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
731: MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
732: if (size ==1) {
733: PetscFree(nonzeros);
734: }
735: }
737: /* Copy values from user jacobian to Ai */
738: MatGetOwnershipRange(ipmP->Ai,&astart,&aend);
740: /* Ai w/lb */
741: if (ipmP->mi) {
742: MatZeroEntries(ipmP->Ai);
743: MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
744: for (i=jstart;i<jend;i++) {
745: MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
746: newrow = i;
747: MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
748: MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
749: }
750: }
752: /* I w/ xlb */
753: if (ipmP->nxlb) {
754: for (i=0;i<ipmP->nxlb;i++) {
755: if (i>=astart && i<aend) {
756: newrow = i+r2;
757: newcol = i;
758: newval = 1.0;
759: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
760: }
761: }
762: }
763: if (ipmP->nxub) {
764: /* I w/ xub */
765: for (i=0;i<ipmP->nxub;i++) {
766: if (i>=astart && i<aend) {
767: newrow = i+r3;
768: newcol = i;
769: newval = -1.0;
770: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
771: }
772: }
773: }
775: MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
776: MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
777: CHKMEMQ;
779: VecSet(ipmP->ci,0.0);
781: /* user ci */
782: if (ipmP->mi > 0) {
783: VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
784: VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
785: }
786: if (!ipmP->work){
787: VecDuplicate(tao->solution,&ipmP->work);
788: }
789: VecCopy(tao->solution,ipmP->work);
790: if (tao->XL) {
791: VecAXPY(ipmP->work,-1.0,tao->XL);
793: /* lower bounds on variables */
794: if (ipmP->nxlb > 0) {
795: VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
796: VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
797: }
798: }
799: if (tao->XU) {
800: /* upper bounds on variables */
801: VecCopy(tao->solution,ipmP->work);
802: VecScale(ipmP->work,-1.0);
803: VecAXPY(ipmP->work,1.0,tao->XU);
804: if (ipmP->nxub > 0) {
805: VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
806: VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
807: }
808: }
809: return(0);
810: }
812: /* create K = [ Hlag , 0 , Ae', -Ai'];
813: [Ae , 0, 0 , 0];
814: [Ai ,-I, 0 , 0];
815: [ 0 , S , 0, Y ]; */
816: PetscErrorCode IPMUpdateK(Tao tao)
817: {
818: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
819: MPI_Comm comm;
820: PetscMPIInt size;
821: PetscErrorCode ierr;
822: PetscInt i,j,row;
823: PetscInt ncols,newcol,newcols[2],newrow;
824: const PetscInt *cols;
825: const PetscReal *vals;
826: const PetscReal *l,*y;
827: PetscReal *newvals;
828: PetscReal newval;
829: PetscInt subsize;
830: const PetscInt *indices;
831: PetscInt *nonzeros,*d_nonzeros,*o_nonzeros;
832: PetscInt bigsize;
833: PetscInt r1,r2,r3;
834: PetscInt c1,c2,c3;
835: PetscInt klocalsize;
836: PetscInt hstart,hend,kstart,kend;
837: PetscInt aistart,aiend,aestart,aeend;
838: PetscInt sstart,send;
841: comm = ((PetscObject)(tao->solution))->comm;
842: MPI_Comm_size(comm,&size);
843: IPMUpdateAi(tao);
845: /* allocate workspace */
846: subsize = PetscMax(ipmP->n,ipmP->nb);
847: subsize = PetscMax(ipmP->me,subsize);
848: subsize = PetscMax(2,subsize);
849: PetscMalloc1(subsize,&indices);
850: PetscMalloc1(subsize,&newvals);
852: r1 = c1 = ipmP->n;
853: r2 = r1 + ipmP->me; c2 = c1 + ipmP->nb;
854: r3 = c3 = r2 + ipmP->nb;
856: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
857: VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
858: MatGetOwnershipRange(tao->hessian,&hstart,&hend);
859: klocalsize = kend-kstart;
860: if (!ipmP->K) {
861: if (size == 1) {
862: PetscMalloc1(kend-kstart,&nonzeros);
863: for (i=0;i<bigsize;i++) {
864: if (i<r1) {
865: MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
866: nonzeros[i] = ncols;
867: MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
868: nonzeros[i] += ipmP->me+ipmP->nb;
869: } else if (i<r2) {
870: nonzeros[i-kstart] = ipmP->n;
871: } else if (i<r3) {
872: nonzeros[i-kstart] = ipmP->n+1;
873: } else if (i<bigsize) {
874: nonzeros[i-kstart] = 2;
875: }
876: }
877: MatCreate(comm,&ipmP->K);
878: MatSetType(ipmP->K,MATSEQAIJ);
879: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
880: MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
881: MatSetFromOptions(ipmP->K);
882: PetscFree(nonzeros);
883: } else {
884: PetscMalloc1(kend-kstart,&d_nonzeros);
885: PetscMalloc1(kend-kstart,&o_nonzeros);
886: for (i=kstart;i<kend;i++) {
887: if (i<r1) {
888: /* TODO fix preallocation for mpi mats */
889: d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
890: o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
891: } else if (i<r2) {
892: d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
893: o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
894: } else if (i<r3) {
895: d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
896: o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
897: } else {
898: d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
899: o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
900: }
901: }
902: MatCreate(comm,&ipmP->K);
903: MatSetType(ipmP->K,MATMPIAIJ);
904: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
905: MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
906: PetscFree(d_nonzeros);
907: PetscFree(o_nonzeros);
908: MatSetFromOptions(ipmP->K);
909: }
910: }
912: MatZeroEntries(ipmP->K);
913: /* Copy H */
914: for (i=hstart;i<hend;i++) {
915: MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
916: if (ncols > 0) {
917: MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
918: }
919: MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
920: }
922: /* Copy Ae and Ae' */
923: if (ipmP->me > 0) {
924: MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
925: for (i=aestart;i<aeend;i++) {
926: MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
927: if (ncols > 0) {
928: /*Ae*/
929: row = i+r1;
930: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
931: /*Ae'*/
932: for (j=0;j<ncols;j++) {
933: newcol = i + c2;
934: newrow = cols[j];
935: newval = vals[j];
936: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
937: }
938: }
939: MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
940: }
941: }
943: if (ipmP->nb > 0) {
944: MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
945: /* Copy Ai,and Ai' */
946: for (i=aistart;i<aiend;i++) {
947: row = i+r2;
948: MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
949: if (ncols > 0) {
950: /*Ai*/
951: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
952: /*-Ai'*/
953: for (j=0;j<ncols;j++) {
954: newcol = i + c3;
955: newrow = cols[j];
956: newval = -vals[j];
957: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
958: }
959: }
960: MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
961: }
963: /* -I */
964: for (i=kstart;i<kend;i++) {
965: if (i>=r2 && i<r3) {
966: newrow = i;
967: newcol = i-r2+c1;
968: newval = -1.0;
969: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
970: }
971: }
973: /* Copy L,Y */
974: VecGetOwnershipRange(ipmP->s,&sstart,&send);
975: VecGetArrayRead(ipmP->lamdai,&l);
976: VecGetArrayRead(ipmP->s,&y);
978: for (i=sstart;i<send;i++) {
979: newcols[0] = c1+i;
980: newcols[1] = c3+i;
981: newvals[0] = l[i-sstart];
982: newvals[1] = y[i-sstart];
983: newrow = r3+i;
984: MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
985: }
987: VecRestoreArrayRead(ipmP->lamdai,&l);
988: VecRestoreArrayRead(ipmP->s,&y);
989: }
991: PetscFree(indices);
992: PetscFree(newvals);
993: MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
994: MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
995: return(0);
996: }
998: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
999: {
1000: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1004: /* rhs = [x1 (n)
1005: x2 (me)
1006: x3 (nb)
1007: x4 (nb)] */
1008: if (X1) {
1009: VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1010: VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1011: }
1012: if (ipmP->me > 0 && X2) {
1013: VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1014: VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1015: }
1016: if (ipmP->nb > 0) {
1017: if (X3) {
1018: VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1019: VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1020: }
1021: if (X4) {
1022: VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1023: VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1024: }
1025: }
1026: return(0);
1027: }
1029: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1030: {
1031: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1035: CHKMEMQ;
1036: /* [x1 (n)
1037: x2 (nb) may be 0
1038: x3 (me) may be 0
1039: x4 (nb) may be 0 */
1040: if (X1) {
1041: VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1042: VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1043: }
1044: if (X2 && ipmP->nb > 0) {
1045: VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1046: VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1047: }
1048: if (X3 && ipmP->me > 0) {
1049: VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1050: VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1051: }
1052: if (X4 && ipmP->nb > 0) {
1053: VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1054: VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1055: }
1056: CHKMEMQ;
1057: return(0);
1058: }
1060: /*MC
1061: TAOIPM - Interior point algorithm for generally constrained optimization.
1063: Option Database Keys:
1064: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1065: . -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1067: Notes: This algorithm is more of a place-holder for future constrained optimization algorithms and should not yet be used for large problems or production code.
1068: Level: beginner
1070: M*/
1072: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1073: {
1074: TAO_IPM *ipmP;
1078: tao->ops->setup = TaoSetup_IPM;
1079: tao->ops->solve = TaoSolve_IPM;
1080: tao->ops->view = TaoView_IPM;
1081: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1082: tao->ops->destroy = TaoDestroy_IPM;
1083: /* tao->ops->computedual = TaoComputeDual_IPM; */
1085: PetscNewLog(tao,&ipmP);
1086: tao->data = (void*)ipmP;
1088: /* Override default settings (unless already changed) */
1089: if (!tao->max_it_changed) tao->max_it = 200;
1090: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1092: ipmP->dec = 10000; /* line search critera */
1093: ipmP->taumin = 0.995;
1094: ipmP->monitorkkt = PETSC_FALSE;
1095: ipmP->pushs = 100;
1096: ipmP->pushnu = 100;
1097: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1098: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1099: return(0);
1100: }