Actual source code: ipm.c
petsc-3.5.4 2015-05-23
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/constrained/impls/ipm/ipm.h> /*I "ipm.h" I*/
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);
35: static PetscErrorCode TaoSolve_IPM(Tao tao)
36: {
37: PetscErrorCode ierr;
38: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
39: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
40: PetscInt iter = 0,its,i;
41: PetscScalar stepsize=1.0;
42: PetscScalar step_s,step_l,alpha,tau,sigma,phi_target;
45: /* Push initial point away from bounds */
46: IPMInitializeBounds(tao);
47: IPMPushInitialPoint(tao);
48: VecCopy(tao->solution,ipmP->rhs_x);
49: IPMEvaluate(tao);
50: IPMComputeKKT(tao);
51: TaoMonitor(tao,iter++,ipmP->kkt_f,ipmP->phi,0.0,1.0,&reason);
53: while (reason == TAO_CONTINUE_ITERATING) {
54: IPMUpdateK(tao);
55: /*
56: rhs.x = -rd
57: rhs.lame = -rpe
58: rhs.lami = -rpi
59: rhs.com = -com
60: */
62: VecCopy(ipmP->rd,ipmP->rhs_x);
63: if (ipmP->me > 0) {
64: VecCopy(ipmP->rpe,ipmP->rhs_lamdae);
65: }
66: if (ipmP->nb > 0) {
67: VecCopy(ipmP->rpi,ipmP->rhs_lamdai);
68: VecCopy(ipmP->complementarity,ipmP->rhs_s);
69: }
70: IPMGatherRHS(tao,ipmP->bigrhs,ipmP->rhs_x,ipmP->rhs_lamdae,ipmP->rhs_lamdai,ipmP->rhs_s);
71: VecScale(ipmP->bigrhs,-1.0);
73: /* solve K * step = rhs */
74: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
75: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
77: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
78: KSPGetIterationNumber(tao->ksp,&its);
79: tao->ksp_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;
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,iter,ipmP->kkt_f,ipmP->phi,0.0,stepsize,&reason);
191: iter++;
192: }
193: return(0);
194: }
198: static PetscErrorCode TaoSetup_IPM(Tao tao)
199: {
200: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
204: ipmP->nb = ipmP->mi = ipmP->me = 0;
205: ipmP->K=0;
206: VecGetSize(tao->solution,&ipmP->n);
207: if (!tao->gradient) {
208: VecDuplicate(tao->solution, &tao->gradient);
209: VecDuplicate(tao->solution, &tao->stepdirection);
210: VecDuplicate(tao->solution, &ipmP->rd);
211: VecDuplicate(tao->solution, &ipmP->rhs_x);
212: VecDuplicate(tao->solution, &ipmP->work);
213: VecDuplicate(tao->solution, &ipmP->save_x);
214: }
215: if (tao->constraints_equality) {
216: VecGetSize(tao->constraints_equality,&ipmP->me);
217: VecDuplicate(tao->constraints_equality,&ipmP->lamdae);
218: VecDuplicate(tao->constraints_equality,&ipmP->dlamdae);
219: VecDuplicate(tao->constraints_equality,&ipmP->rhs_lamdae);
220: VecDuplicate(tao->constraints_equality,&ipmP->save_lamdae);
221: VecDuplicate(tao->constraints_equality,&ipmP->rpe);
222: VecDuplicate(tao->constraints_equality,&tao->DE);
223: }
224: if (tao->constraints_inequality) {
225: VecDuplicate(tao->constraints_inequality,&tao->DI);
226: }
227: return(0);
228: }
232: static PetscErrorCode IPMInitializeBounds(Tao tao)
233: {
234: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
235: Vec xtmp;
236: PetscInt xstart,xend;
237: PetscInt ucstart,ucend; /* user ci */
238: PetscInt ucestart,uceend; /* user ce */
239: PetscInt sstart,send;
240: PetscInt bigsize;
241: PetscInt i,counter,nloc;
242: PetscInt *cind,*xind,*ucind,*uceind,*stepind;
243: VecType vtype;
244: const PetscInt *xli,*xui;
245: PetscInt xl_offset,xu_offset;
246: IS bigxl,bigxu,isuc,isc,isx,sis,is1;
248: MPI_Comm comm;
251: cind=xind=ucind=uceind=stepind=0;
252: ipmP->mi=0;
253: ipmP->nxlb=0;
254: ipmP->nxub=0;
255: ipmP->nb=0;
256: ipmP->nslack=0;
258: VecDuplicate(tao->solution,&xtmp);
259: if (!tao->XL && !tao->XU && tao->ops->computebounds) {
260: TaoComputeVariableBounds(tao);
261: }
262: if (tao->XL) {
263: VecSet(xtmp,PETSC_NINFINITY);
264: VecWhichGreaterThan(tao->XL,xtmp,&ipmP->isxl);
265: ISGetSize(ipmP->isxl,&ipmP->nxlb);
266: } else {
267: ipmP->nxlb=0;
268: }
269: if (tao->XU) {
270: VecSet(xtmp,PETSC_INFINITY);
271: VecWhichLessThan(tao->XU,xtmp,&ipmP->isxu);
272: ISGetSize(ipmP->isxu,&ipmP->nxub);
273: } else {
274: ipmP->nxub=0;
275: }
276: VecDestroy(&xtmp);
277: if (tao->constraints_inequality) {
278: VecGetSize(tao->constraints_inequality,&ipmP->mi);
279: } else {
280: ipmP->mi = 0;
281: }
282: ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;
284: comm = ((PetscObject)(tao->solution))->comm;
286: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
287: PetscMalloc1(bigsize,&stepind);
288: PetscMalloc1(ipmP->n,&xind);
289: PetscMalloc1(ipmP->me,&uceind);
290: VecGetOwnershipRange(tao->solution,&xstart,&xend);
292: if (ipmP->nb > 0) {
293: VecCreate(comm,&ipmP->s);
294: VecSetSizes(ipmP->s,PETSC_DECIDE,ipmP->nb);
295: VecSetFromOptions(ipmP->s);
296: VecDuplicate(ipmP->s,&ipmP->ds);
297: VecDuplicate(ipmP->s,&ipmP->rhs_s);
298: VecDuplicate(ipmP->s,&ipmP->complementarity);
299: VecDuplicate(ipmP->s,&ipmP->ci);
301: VecDuplicate(ipmP->s,&ipmP->lamdai);
302: VecDuplicate(ipmP->s,&ipmP->dlamdai);
303: VecDuplicate(ipmP->s,&ipmP->rhs_lamdai);
304: VecDuplicate(ipmP->s,&ipmP->save_lamdai);
306: VecDuplicate(ipmP->s,&ipmP->save_s);
307: VecDuplicate(ipmP->s,&ipmP->rpi);
308: VecDuplicate(ipmP->s,&ipmP->Zero_nb);
309: VecSet(ipmP->Zero_nb,0.0);
310: VecDuplicate(ipmP->s,&ipmP->One_nb);
311: VecSet(ipmP->One_nb,1.0);
312: VecDuplicate(ipmP->s,&ipmP->Inf_nb);
313: VecSet(ipmP->Inf_nb,PETSC_INFINITY);
315: PetscMalloc1(ipmP->nb,&cind);
316: PetscMalloc1(ipmP->mi,&ucind);
317: VecGetOwnershipRange(ipmP->s,&sstart,&send);
319: if (ipmP->mi > 0) {
320: VecGetOwnershipRange(tao->constraints_inequality,&ucstart,&ucend);
321: counter=0;
322: for (i=ucstart;i<ucend;i++) {
323: cind[counter++] = i;
324: }
325: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isuc);
326: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isc);
327: VecScatterCreate(tao->constraints_inequality,isuc,ipmP->ci,isc,&ipmP->ci_scat);
329: ISDestroy(&isuc);
330: ISDestroy(&isc);
331: }
332: /* need to know how may xbound indices are on each process */
333: /* TODO better way */
334: if (ipmP->nxlb) {
335: ISAllGather(ipmP->isxl,&bigxl);
336: ISGetIndices(bigxl,&xli);
337: /* find offsets for this processor */
338: xl_offset = ipmP->mi;
339: for (i=0;i<ipmP->nxlb;i++) {
340: if (xli[i] < xstart) {
341: xl_offset++;
342: } else break;
343: }
344: ISRestoreIndices(bigxl,&xli);
346: ISGetIndices(ipmP->isxl,&xli);
347: ISGetLocalSize(ipmP->isxl,&nloc);
348: for (i=0;i<nloc;i++) {
349: xind[i] = xli[i];
350: cind[i] = xl_offset+i;
351: }
353: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
354: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
355: VecScatterCreate(tao->XL,isx,ipmP->ci,isc,&ipmP->xl_scat);
356: ISDestroy(&isx);
357: ISDestroy(&isc);
358: ISDestroy(&bigxl);
359: }
361: if (ipmP->nxub) {
362: ISAllGather(ipmP->isxu,&bigxu);
363: ISGetIndices(bigxu,&xui);
364: /* find offsets for this processor */
365: xu_offset = ipmP->mi + ipmP->nxlb;
366: for (i=0;i<ipmP->nxub;i++) {
367: if (xui[i] < xstart) {
368: xu_offset++;
369: } else break;
370: }
371: ISRestoreIndices(bigxu,&xui);
373: ISGetIndices(ipmP->isxu,&xui);
374: ISGetLocalSize(ipmP->isxu,&nloc);
375: for (i=0;i<nloc;i++) {
376: xind[i] = xui[i];
377: cind[i] = xu_offset+i;
378: }
380: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
381: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
382: VecScatterCreate(tao->XU,isx,ipmP->ci,isc,&ipmP->xu_scat);
383: ISDestroy(&isx);
384: ISDestroy(&isc);
385: ISDestroy(&bigxu);
386: }
387: }
388: VecCreate(comm,&ipmP->bigrhs);
389: VecGetType(tao->solution,&vtype);
390: VecSetType(ipmP->bigrhs,vtype);
391: VecSetSizes(ipmP->bigrhs,PETSC_DECIDE,bigsize);
392: VecSetFromOptions(ipmP->bigrhs);
393: VecDuplicate(ipmP->bigrhs,&ipmP->bigstep);
395: /* create scatters for step->x and x->rhs */
396: for (i=xstart;i<xend;i++) {
397: stepind[i-xstart] = i;
398: xind[i-xstart] = i;
399: }
400: ISCreateGeneral(comm,xend-xstart,stepind,PETSC_COPY_VALUES,&sis);
401: ISCreateGeneral(comm,xend-xstart,xind,PETSC_COPY_VALUES,&is1);
402: VecScatterCreate(ipmP->bigstep,sis,tao->solution,is1,&ipmP->step1);
403: VecScatterCreate(tao->solution,is1,ipmP->bigrhs,sis,&ipmP->rhs1);
404: ISDestroy(&sis);
405: ISDestroy(&is1);
407: if (ipmP->nb > 0) {
408: for (i=sstart;i<send;i++) {
409: stepind[i-sstart] = i+ipmP->n;
410: cind[i-sstart] = i;
411: }
412: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
413: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
414: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step2);
415: ISDestroy(&sis);
417: for (i=sstart;i<send;i++) {
418: stepind[i-sstart] = i+ipmP->n+ipmP->me;
419: cind[i-sstart] = i;
420: }
421: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
422: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs3);
423: ISDestroy(&sis);
424: ISDestroy(&is1);
425: }
427: if (ipmP->me > 0) {
428: VecGetOwnershipRange(tao->constraints_equality,&ucestart,&uceend);
429: for (i=ucestart;i<uceend;i++) {
430: stepind[i-ucestart] = i + ipmP->n+ipmP->nb;
431: uceind[i-ucestart] = i;
432: }
434: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
435: ISCreateGeneral(comm,uceend-ucestart,uceind,PETSC_COPY_VALUES,&is1);
436: VecScatterCreate(ipmP->bigstep,sis,tao->constraints_equality,is1,&ipmP->step3);
437: ISDestroy(&sis);
439: for (i=ucestart;i<uceend;i++) {
440: stepind[i-ucestart] = i + ipmP->n;
441: }
443: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
444: VecScatterCreate(tao->constraints_equality,is1,ipmP->bigrhs,sis,&ipmP->rhs2);
445: ISDestroy(&sis);
446: ISDestroy(&is1);
447: }
449: if (ipmP->nb > 0) {
450: for (i=sstart;i<send;i++) {
451: stepind[i-sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
452: cind[i-sstart] = i;
453: }
454: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
455: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
456: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step4);
457: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs4);
458: ISDestroy(&sis);
459: ISDestroy(&is1);
460: }
462: PetscFree(stepind);
463: PetscFree(cind);
464: PetscFree(ucind);
465: PetscFree(uceind);
466: PetscFree(xind);
467: return(0);
468: }
472: static PetscErrorCode TaoDestroy_IPM(Tao tao)
473: {
474: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
478: VecDestroy(&ipmP->rd);
479: VecDestroy(&ipmP->rpe);
480: VecDestroy(&ipmP->rpi);
481: VecDestroy(&ipmP->work);
482: VecDestroy(&ipmP->lamdae);
483: VecDestroy(&ipmP->lamdai);
484: VecDestroy(&ipmP->s);
485: VecDestroy(&ipmP->ds);
486: VecDestroy(&ipmP->ci);
488: VecDestroy(&ipmP->rhs_x);
489: VecDestroy(&ipmP->rhs_lamdae);
490: VecDestroy(&ipmP->rhs_lamdai);
491: VecDestroy(&ipmP->rhs_s);
493: VecDestroy(&ipmP->save_x);
494: VecDestroy(&ipmP->save_lamdae);
495: VecDestroy(&ipmP->save_lamdai);
496: VecDestroy(&ipmP->save_s);
498: VecScatterDestroy(&ipmP->step1);
499: VecScatterDestroy(&ipmP->step2);
500: VecScatterDestroy(&ipmP->step3);
501: VecScatterDestroy(&ipmP->step4);
503: VecScatterDestroy(&ipmP->rhs1);
504: VecScatterDestroy(&ipmP->rhs2);
505: VecScatterDestroy(&ipmP->rhs3);
506: VecScatterDestroy(&ipmP->rhs4);
508: VecScatterDestroy(&ipmP->ci_scat);
509: VecScatterDestroy(&ipmP->xl_scat);
510: VecScatterDestroy(&ipmP->xu_scat);
512: VecDestroy(&ipmP->dlamdai);
513: VecDestroy(&ipmP->dlamdae);
514: VecDestroy(&ipmP->Zero_nb);
515: VecDestroy(&ipmP->One_nb);
516: VecDestroy(&ipmP->Inf_nb);
517: VecDestroy(&ipmP->complementarity);
519: VecDestroy(&ipmP->bigrhs);
520: VecDestroy(&ipmP->bigstep);
521: MatDestroy(&ipmP->Ai);
522: MatDestroy(&ipmP->K);
523: ISDestroy(&ipmP->isxu);
524: ISDestroy(&ipmP->isxl);
525: PetscFree(tao->data);
526: return(0);
527: }
531: static PetscErrorCode TaoSetFromOptions_IPM(Tao tao)
532: {
533: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
535: PetscBool flg;
538: PetscOptionsHead("IPM method for constrained optimization");
539: PetscOptionsBool("-tao_ipm_monitorkkt","monitor kkt status",NULL,ipmP->monitorkkt,&ipmP->monitorkkt,&flg);
540: PetscOptionsReal("-tao_ipm_pushs","parameter to push initial slack variables away from bounds",NULL,ipmP->pushs,&ipmP->pushs,&flg);
541: PetscOptionsReal("-tao_ipm_pushnu","parameter to push initial (inequality) dual variables away from bounds",NULL,ipmP->pushnu,&ipmP->pushnu,&flg);
542: PetscOptionsTail();
543: ierr =KSPSetFromOptions(tao->ksp);
544: return(0);
545: }
549: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
550: {
551: return 0;
552: }
554: /* IPMObjectiveAndGradient()
555: f = d'x + 0.5 * x' * H * x
556: rd = H*x + d + Ae'*lame - Ai'*lami
557: rpe = Ae*x - be
558: rpi = Ai*x - yi - bi
559: mu = yi' * lami/mi;
560: com = yi.*lami
562: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
563: */
564: /*
567: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
568: {
569: Tao tao = (Tao)tptr;
570: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
573: IPMComputeKKT(tao);
574: *f = ipmP->phi;
575: return(0);
576: }
577: */
579: /*
580: f = d'x + 0.5 * x' * H * x
581: rd = H*x + d + Ae'*lame - Ai'*lami
582: Ai = jac_ineq
583: I (w/lb)
584: -I (w/ub)
586: rpe = ce
587: rpi = ci - s;
588: com = s.*lami
589: mu = yi' * lami/mi;
591: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
592: */
595: static PetscErrorCode IPMComputeKKT(Tao tao)
596: {
597: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
598: PetscScalar norm;
602: VecCopy(tao->gradient,ipmP->rd);
604: if (ipmP->me > 0) {
605: /* rd = gradient + Ae'*lamdae */
606: MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
607: VecAXPY(ipmP->rd, 1.0, ipmP->work);
609: /* rpe = ce(x) */
610: VecCopy(tao->constraints_equality,ipmP->rpe);
611: }
612: if (ipmP->nb > 0) {
613: /* rd = rd - Ai'*lamdai */
614: MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
615: VecAXPY(ipmP->rd, -1.0, ipmP->work);
617: /* rpi = cin - s */
618: VecCopy(ipmP->ci,ipmP->rpi);
619: VecAXPY(ipmP->rpi, -1.0, ipmP->s);
621: /* com = s .* lami */
622: VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
623: }
624: /* phi = ||rd; rpe; rpi; com|| */
625: VecDot(ipmP->rd,ipmP->rd,&norm);
626: ipmP->phi = norm;
627: if (ipmP->me > 0 ) {
628: VecDot(ipmP->rpe,ipmP->rpe,&norm);
629: ipmP->phi += norm;
630: }
631: if (ipmP->nb > 0) {
632: VecDot(ipmP->rpi,ipmP->rpi,&norm);
633: ipmP->phi += norm;
634: VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
635: ipmP->phi += norm;
636: /* mu = s'*lami/nb */
637: VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
638: ipmP->mu /= ipmP->nb;
639: } else {
640: ipmP->mu = 1.0;
641: }
643: ipmP->phi = PetscSqrtScalar(ipmP->phi);
644: return(0);
645: }
649: /* evaluate user info at current point */
650: PetscErrorCode IPMEvaluate(Tao tao)
651: {
652: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
656: TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
657: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
658: if (ipmP->me > 0) {
659: TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
660: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
661: }
662: if (ipmP->mi > 0) {
663: TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
664: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
665: }
666: if (ipmP->nb > 0) {
667: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
668: IPMUpdateAi(tao);
669: }
670: return(0);
671: }
675: /* Push initial point away from bounds */
676: PetscErrorCode IPMPushInitialPoint(Tao tao)
677: {
678: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
682: TaoComputeVariableBounds(tao);
683: if (tao->XL && tao->XU) {
684: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
685: }
686: if (ipmP->nb > 0) {
687: VecSet(ipmP->s,ipmP->pushs);
688: VecSet(ipmP->lamdai,ipmP->pushnu);
689: if (ipmP->mi > 0) {
690: VecSet(tao->DI,ipmP->pushnu);
691: }
692: }
693: if (ipmP->me > 0) {
694: VecSet(tao->DE,1.0);
695: VecSet(ipmP->lamdae,1.0);
696: }
697: return(0);
698: }
702: PetscErrorCode IPMUpdateAi(Tao tao)
703: {
704: /* Ai = Ji
705: I (w/lb)
706: -I (w/ub) */
708: /* Ci = user->ci
709: Xi - lb (w/lb)
710: -Xi + ub (w/ub) */
712: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
713: MPI_Comm comm;
714: PetscInt i;
715: PetscScalar newval;
716: PetscInt newrow,newcol,ncols;
717: const PetscScalar *vals;
718: const PetscInt *cols;
719: PetscInt astart,aend,jstart,jend;
720: PetscInt *nonzeros;
721: PetscInt r2,r3,r4;
722: PetscMPIInt size;
723: PetscErrorCode ierr;
726: r2 = ipmP->mi;
727: r3 = r2 + ipmP->nxlb;
728: r4 = r3 + ipmP->nxub;
730: if (!ipmP->nb) return(0);
732: /* Create Ai matrix if it doesn't exist yet */
733: if (!ipmP->Ai) {
734: comm = ((PetscObject)(tao->solution))->comm;
735: PetscMalloc1(ipmP->nb,&nonzeros);
736: MPI_Comm_size(comm,&size);
737: if (size == 1) {
738: for (i=0;i<ipmP->mi;i++) {
739: MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
740: nonzeros[i] = ncols;
741: MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
742: }
743: for (i=r2;i<r4;i++) {
744: nonzeros[i] = 1;
745: }
746: }
747: MatCreate(comm,&ipmP->Ai);
748: MatSetType(ipmP->Ai,MATAIJ);
749: MatSetSizes(ipmP->Ai,PETSC_DECIDE,PETSC_DECIDE,ipmP->nb,ipmP->n);
750: MatSetFromOptions(ipmP->Ai);
751: MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
752: MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
753: if (size ==1) {
754: PetscFree(nonzeros);
755: }
756: }
758: /* Copy values from user jacobian to Ai */
759: MatGetOwnershipRange(ipmP->Ai,&astart,&aend);
761: /* Ai w/lb */
762: if (ipmP->mi) {
763: MatZeroEntries(ipmP->Ai);
764: MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
765: for (i=jstart;i<jend;i++) {
766: MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
767: newrow = i;
768: MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
769: MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
770: }
771: }
773: /* I w/ xlb */
774: if (ipmP->nxlb) {
775: for (i=0;i<ipmP->nxlb;i++) {
776: if (i>=astart && i<aend) {
777: newrow = i+r2;
778: newcol = i;
779: newval = 1.0;
780: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
781: }
782: }
783: }
784: if (ipmP->nxub) {
785: /* I w/ xub */
786: for (i=0;i<ipmP->nxub;i++) {
787: if (i>=astart && i<aend) {
788: newrow = i+r3;
789: newcol = i;
790: newval = -1.0;
791: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
792: }
793: }
794: }
796: MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
797: MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
798: CHKMEMQ;
800: VecSet(ipmP->ci,0.0);
802: /* user ci */
803: if (ipmP->mi > 0) {
804: VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
805: VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
806: }
807: if (!ipmP->work){
808: VecDuplicate(tao->solution,&ipmP->work);
809: }
810: VecCopy(tao->solution,ipmP->work);
811: if (tao->XL) {
812: VecAXPY(ipmP->work,-1.0,tao->XL);
814: /* lower bounds on variables */
815: if (ipmP->nxlb > 0) {
816: VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
817: VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
818: }
819: }
820: if (tao->XU) {
821: /* upper bounds on variables */
822: VecCopy(tao->solution,ipmP->work);
823: VecScale(ipmP->work,-1.0);
824: VecAXPY(ipmP->work,1.0,tao->XU);
825: if (ipmP->nxub > 0) {
826: VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
827: VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
828: }
829: }
830: return(0);
831: }
835: /* create K = [ Hlag , 0 , Ae', -Ai'];
836: [Ae , 0, 0 , 0];
837: [Ai ,-I, 0 , 0];
838: [ 0 , S , 0, Y ]; */
839: PetscErrorCode IPMUpdateK(Tao tao)
840: {
841: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
842: MPI_Comm comm;
843: PetscMPIInt size;
844: PetscErrorCode ierr;
845: PetscInt i,j,row;
846: PetscInt ncols,newcol,newcols[2],newrow;
847: const PetscInt *cols;
848: const PetscReal *vals;
849: const PetscReal *l,*y;
850: PetscReal *newvals;
851: PetscReal newval;
852: PetscInt subsize;
853: const PetscInt *indices;
854: PetscInt *nonzeros,*d_nonzeros,*o_nonzeros;
855: PetscInt bigsize;
856: PetscInt r1,r2,r3;
857: PetscInt c1,c2,c3;
858: PetscInt klocalsize;
859: PetscInt hstart,hend,kstart,kend;
860: PetscInt aistart,aiend,aestart,aeend;
861: PetscInt sstart,send;
864: comm = ((PetscObject)(tao->solution))->comm;
865: MPI_Comm_size(comm,&size);
866: IPMUpdateAi(tao);
868: /* allocate workspace */
869: subsize = PetscMax(ipmP->n,ipmP->nb);
870: subsize = PetscMax(ipmP->me,subsize);
871: subsize = PetscMax(2,subsize);
872: PetscMalloc1(subsize,&indices);
873: PetscMalloc1(subsize,&newvals);
875: r1 = c1 = ipmP->n;
876: r2 = r1 + ipmP->me; c2 = c1 + ipmP->nb;
877: r3 = c3 = r2 + ipmP->nb;
879: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
880: VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
881: MatGetOwnershipRange(tao->hessian,&hstart,&hend);
882: klocalsize = kend-kstart;
883: if (!ipmP->K) {
884: if (size == 1) {
885: PetscMalloc1((kend-kstart),&nonzeros);
886: for (i=0;i<bigsize;i++) {
887: if (i<r1) {
888: MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
889: nonzeros[i] = ncols;
890: MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
891: nonzeros[i] += ipmP->me+ipmP->nb;
892: } else if (i<r2) {
893: nonzeros[i-kstart] = ipmP->n;
894: } else if (i<r3) {
895: nonzeros[i-kstart] = ipmP->n+1;
896: } else if (i<bigsize) {
897: nonzeros[i-kstart] = 2;
898: }
899: }
900: MatCreate(comm,&ipmP->K);
901: MatSetType(ipmP->K,MATSEQAIJ);
902: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
903: MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
904: MatSetFromOptions(ipmP->K);
905: PetscFree(nonzeros);
906: } else {
907: PetscMalloc1((kend-kstart),&d_nonzeros);
908: PetscMalloc1((kend-kstart),&o_nonzeros);
909: for (i=kstart;i<kend;i++) {
910: if (i<r1) {
911: /* TODO fix preallocation for mpi mats */
912: d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
913: o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
914: } else if (i<r2) {
915: d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
916: o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
917: } else if (i<r3) {
918: d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
919: o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
920: } else {
921: d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
922: o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
923: }
924: }
925: MatCreate(comm,&ipmP->K);
926: MatSetType(ipmP->K,MATMPIAIJ);
927: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
928: MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
929: PetscFree(d_nonzeros);
930: PetscFree(o_nonzeros);
931: MatSetFromOptions(ipmP->K);
932: }
933: }
935: MatZeroEntries(ipmP->K);
936: /* Copy H */
937: for (i=hstart;i<hend;i++) {
938: MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
939: if (ncols > 0) {
940: MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
941: }
942: MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
943: }
945: /* Copy Ae and Ae' */
946: if (ipmP->me > 0) {
947: MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
948: for (i=aestart;i<aeend;i++) {
949: MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
950: if (ncols > 0) {
951: /*Ae*/
952: row = i+r1;
953: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
954: /*Ae'*/
955: for (j=0;j<ncols;j++) {
956: newcol = i + c2;
957: newrow = cols[j];
958: newval = vals[j];
959: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
960: }
961: }
962: MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
963: }
964: }
966: if (ipmP->nb > 0) {
967: MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
968: /* Copy Ai,and Ai' */
969: for (i=aistart;i<aiend;i++) {
970: row = i+r2;
971: MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
972: if (ncols > 0) {
973: /*Ai*/
974: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
975: /*-Ai'*/
976: for (j=0;j<ncols;j++) {
977: newcol = i + c3;
978: newrow = cols[j];
979: newval = -vals[j];
980: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
981: }
982: }
983: MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
984: }
986: /* -I */
987: for (i=kstart;i<kend;i++) {
988: if (i>=r2 && i<r3) {
989: newrow = i;
990: newcol = i-r2+c1;
991: newval = -1.0;
992: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
993: }
994: }
996: /* Copy L,Y */
997: VecGetOwnershipRange(ipmP->s,&sstart,&send);
998: VecGetArrayRead(ipmP->lamdai,&l);
999: VecGetArrayRead(ipmP->s,&y);
1001: for (i=sstart;i<send;i++) {
1002: newcols[0] = c1+i;
1003: newcols[1] = c3+i;
1004: newvals[0] = l[i-sstart];
1005: newvals[1] = y[i-sstart];
1006: newrow = r3+i;
1007: MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
1008: }
1010: VecRestoreArrayRead(ipmP->lamdai,&l);
1011: VecRestoreArrayRead(ipmP->s,&y);
1012: }
1014: PetscFree(indices);
1015: PetscFree(newvals);
1016: MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
1017: MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
1018: return(0);
1019: }
1023: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
1024: {
1025: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1029: /* rhs = [x1 (n)
1030: x2 (me)
1031: x3 (nb)
1032: x4 (nb)] */
1033: if (X1) {
1034: VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1035: VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1036: }
1037: if (ipmP->me > 0 && X2) {
1038: VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1039: VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1040: }
1041: if (ipmP->nb > 0) {
1042: if (X3) {
1043: VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1044: VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1045: }
1046: if (X4) {
1047: VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1048: VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1049: }
1050: }
1051: return(0);
1052: }
1056: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1057: {
1058: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1062: CHKMEMQ;
1063: /* [x1 (n)
1064: x2 (nb) may be 0
1065: x3 (me) may be 0
1066: x4 (nb) may be 0 */
1067: if (X1) {
1068: VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1069: VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1070: }
1071: if (X2 && ipmP->nb > 0) {
1072: VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1073: VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1074: }
1075: if (X3 && ipmP->me > 0) {
1076: VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1077: VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1078: }
1079: if (X4 && ipmP->nb > 0) {
1080: VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1081: VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1082: }
1083: CHKMEMQ;
1084: return(0);
1085: }
1087: /*MC
1088: TAOIPM - Interior point algorithm for generally constrained optimization.
1090: Option Database Keys:
1091: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1092: . -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1094: 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.
1095: Level: beginner
1097: M*/
1099: EXTERN_C_BEGIN
1102: PetscErrorCode TaoCreate_IPM(Tao tao)
1103: {
1104: TAO_IPM *ipmP;
1108: tao->ops->setup = TaoSetup_IPM;
1109: tao->ops->solve = TaoSolve_IPM;
1110: tao->ops->view = TaoView_IPM;
1111: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1112: tao->ops->destroy = TaoDestroy_IPM;
1113: /* tao->ops->computedual = TaoComputeDual_IPM; */
1115: PetscNewLog(tao,&ipmP);
1116: tao->data = (void*)ipmP;
1117: tao->max_it = 200;
1118: tao->max_funcs = 500;
1119: tao->fatol = 1e-4;
1120: tao->frtol = 1e-4;
1121: ipmP->dec = 10000; /* line search critera */
1122: ipmP->taumin = 0.995;
1123: ipmP->monitorkkt = PETSC_FALSE;
1124: ipmP->pushs = 100;
1125: ipmP->pushnu = 100;
1126: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1127: return(0);
1128: }
1129: EXTERN_C_END