Actual source code: ipm.c
petsc-3.7.5 2017-01-01
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 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,tao->niter++,ipmP->kkt_f,ipmP->phi,0.0,1.0,&reason);
53: while (reason == TAO_CONTINUE_ITERATING) {
54: tao->ksp_its=0;
55: IPMUpdateK(tao);
56: /*
57: rhs.x = -rd
58: rhs.lame = -rpe
59: rhs.lami = -rpi
60: rhs.com = -com
61: */
63: VecCopy(ipmP->rd,ipmP->rhs_x);
64: if (ipmP->me > 0) {
65: VecCopy(ipmP->rpe,ipmP->rhs_lamdae);
66: }
67: if (ipmP->nb > 0) {
68: VecCopy(ipmP->rpi,ipmP->rhs_lamdai);
69: VecCopy(ipmP->complementarity,ipmP->rhs_s);
70: }
71: IPMGatherRHS(tao,ipmP->bigrhs,ipmP->rhs_x,ipmP->rhs_lamdae,ipmP->rhs_lamdai,ipmP->rhs_s);
72: VecScale(ipmP->bigrhs,-1.0);
74: /* solve K * step = rhs */
75: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
76: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
78: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
79: KSPGetIterationNumber(tao->ksp,&its);
80: tao->ksp_its += its;
81: tao->ksp_tot_its+=its;
82: /* Find distance along step direction to closest bound */
83: if (ipmP->nb > 0) {
84: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
85: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
86: alpha = PetscMin(step_s,step_l);
87: alpha = PetscMin(alpha,1.0);
88: ipmP->alpha1 = alpha;
89: } else {
90: ipmP->alpha1 = alpha = 1.0;
91: }
93: /* x_aff = x + alpha*d */
94: VecCopy(tao->solution,ipmP->save_x);
95: if (ipmP->me > 0) {
96: VecCopy(ipmP->lamdae,ipmP->save_lamdae);
97: }
98: if (ipmP->nb > 0) {
99: VecCopy(ipmP->lamdai,ipmP->save_lamdai);
100: VecCopy(ipmP->s,ipmP->save_s);
101: }
103: VecAXPY(tao->solution,alpha,tao->stepdirection);
104: if (ipmP->me > 0) {
105: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
106: }
107: if (ipmP->nb > 0) {
108: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
109: VecAXPY(ipmP->s,alpha,ipmP->ds);
110: }
112: /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */
113: if (ipmP->mu == 0.0) {
114: sigma = 0.0;
115: } else {
116: sigma = 1.0/ipmP->mu;
117: }
118: IPMComputeKKT(tao);
119: sigma *= ipmP->mu;
120: sigma*=sigma*sigma;
122: /* revert kkt info */
123: VecCopy(ipmP->save_x,tao->solution);
124: if (ipmP->me > 0) {
125: VecCopy(ipmP->save_lamdae,ipmP->lamdae);
126: }
127: if (ipmP->nb > 0) {
128: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
129: VecCopy(ipmP->save_s,ipmP->s);
130: }
131: IPMComputeKKT(tao);
133: /* update rhs with new complementarity vector */
134: if (ipmP->nb > 0) {
135: VecCopy(ipmP->complementarity,ipmP->rhs_s);
136: VecScale(ipmP->rhs_s,-1.0);
137: VecShift(ipmP->rhs_s,sigma*ipmP->mu);
138: }
139: IPMGatherRHS(tao,ipmP->bigrhs,NULL,NULL,NULL,ipmP->rhs_s);
141: /* solve K * step = rhs */
142: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
143: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
145: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
146: KSPGetIterationNumber(tao->ksp,&its);
147: tao->ksp_its += its;
148: tao->ksp_tot_its+=its;
149: if (ipmP->nb > 0) {
150: /* Get max step size and apply frac-to-boundary */
151: tau = PetscMax(ipmP->taumin,1.0-ipmP->mu);
152: tau = PetscMin(tau,1.0);
153: if (tau != 1.0) {
154: VecScale(ipmP->s,tau);
155: VecScale(ipmP->lamdai,tau);
156: }
157: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
158: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
159: if (tau != 1.0) {
160: VecCopy(ipmP->save_s,ipmP->s);
161: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
162: }
163: alpha = PetscMin(step_s,step_l);
164: alpha = PetscMin(alpha,1.0);
165: } else {
166: alpha = 1.0;
167: }
168: ipmP->alpha2 = alpha;
169: /* TODO make phi_target meaningful */
170: phi_target = ipmP->dec * ipmP->phi;
171: for (i=0; i<11;i++) {
172: VecAXPY(tao->solution,alpha,tao->stepdirection);
173: if (ipmP->nb > 0) {
174: VecAXPY(ipmP->s,alpha,ipmP->ds);
175: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
176: }
177: if (ipmP->me > 0) {
178: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
179: }
181: /* update dual variables */
182: if (ipmP->me > 0) {
183: VecCopy(ipmP->lamdae,tao->DE);
184: }
186: IPMEvaluate(tao);
187: IPMComputeKKT(tao);
188: if (ipmP->phi <= phi_target) break;
189: alpha /= 2.0;
190: }
192: TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,stepsize,&reason);
193: tao->niter++;
194: }
195: return(0);
196: }
200: static PetscErrorCode TaoSetup_IPM(Tao tao)
201: {
202: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
206: ipmP->nb = ipmP->mi = ipmP->me = 0;
207: ipmP->K=0;
208: VecGetSize(tao->solution,&ipmP->n);
209: if (!tao->gradient) {
210: VecDuplicate(tao->solution, &tao->gradient);
211: VecDuplicate(tao->solution, &tao->stepdirection);
212: VecDuplicate(tao->solution, &ipmP->rd);
213: VecDuplicate(tao->solution, &ipmP->rhs_x);
214: VecDuplicate(tao->solution, &ipmP->work);
215: VecDuplicate(tao->solution, &ipmP->save_x);
216: }
217: if (tao->constraints_equality) {
218: VecGetSize(tao->constraints_equality,&ipmP->me);
219: VecDuplicate(tao->constraints_equality,&ipmP->lamdae);
220: VecDuplicate(tao->constraints_equality,&ipmP->dlamdae);
221: VecDuplicate(tao->constraints_equality,&ipmP->rhs_lamdae);
222: VecDuplicate(tao->constraints_equality,&ipmP->save_lamdae);
223: VecDuplicate(tao->constraints_equality,&ipmP->rpe);
224: VecDuplicate(tao->constraints_equality,&tao->DE);
225: }
226: if (tao->constraints_inequality) {
227: VecDuplicate(tao->constraints_inequality,&tao->DI);
228: }
229: return(0);
230: }
234: static PetscErrorCode IPMInitializeBounds(Tao tao)
235: {
236: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
237: Vec xtmp;
238: PetscInt xstart,xend;
239: PetscInt ucstart,ucend; /* user ci */
240: PetscInt ucestart,uceend; /* user ce */
241: PetscInt sstart,send;
242: PetscInt bigsize;
243: PetscInt i,counter,nloc;
244: PetscInt *cind,*xind,*ucind,*uceind,*stepind;
245: VecType vtype;
246: const PetscInt *xli,*xui;
247: PetscInt xl_offset,xu_offset;
248: IS bigxl,bigxu,isuc,isc,isx,sis,is1;
250: MPI_Comm comm;
253: cind=xind=ucind=uceind=stepind=0;
254: ipmP->mi=0;
255: ipmP->nxlb=0;
256: ipmP->nxub=0;
257: ipmP->nb=0;
258: ipmP->nslack=0;
260: VecDuplicate(tao->solution,&xtmp);
261: if (!tao->XL && !tao->XU && tao->ops->computebounds) {
262: TaoComputeVariableBounds(tao);
263: }
264: if (tao->XL) {
265: VecSet(xtmp,PETSC_NINFINITY);
266: VecWhichGreaterThan(tao->XL,xtmp,&ipmP->isxl);
267: ISGetSize(ipmP->isxl,&ipmP->nxlb);
268: } else {
269: ipmP->nxlb=0;
270: }
271: if (tao->XU) {
272: VecSet(xtmp,PETSC_INFINITY);
273: VecWhichLessThan(tao->XU,xtmp,&ipmP->isxu);
274: ISGetSize(ipmP->isxu,&ipmP->nxub);
275: } else {
276: ipmP->nxub=0;
277: }
278: VecDestroy(&xtmp);
279: if (tao->constraints_inequality) {
280: VecGetSize(tao->constraints_inequality,&ipmP->mi);
281: } else {
282: ipmP->mi = 0;
283: }
284: ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;
286: comm = ((PetscObject)(tao->solution))->comm;
288: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
289: PetscMalloc1(bigsize,&stepind);
290: PetscMalloc1(ipmP->n,&xind);
291: PetscMalloc1(ipmP->me,&uceind);
292: VecGetOwnershipRange(tao->solution,&xstart,&xend);
294: if (ipmP->nb > 0) {
295: VecCreate(comm,&ipmP->s);
296: VecSetSizes(ipmP->s,PETSC_DECIDE,ipmP->nb);
297: VecSetFromOptions(ipmP->s);
298: VecDuplicate(ipmP->s,&ipmP->ds);
299: VecDuplicate(ipmP->s,&ipmP->rhs_s);
300: VecDuplicate(ipmP->s,&ipmP->complementarity);
301: VecDuplicate(ipmP->s,&ipmP->ci);
303: VecDuplicate(ipmP->s,&ipmP->lamdai);
304: VecDuplicate(ipmP->s,&ipmP->dlamdai);
305: VecDuplicate(ipmP->s,&ipmP->rhs_lamdai);
306: VecDuplicate(ipmP->s,&ipmP->save_lamdai);
308: VecDuplicate(ipmP->s,&ipmP->save_s);
309: VecDuplicate(ipmP->s,&ipmP->rpi);
310: VecDuplicate(ipmP->s,&ipmP->Zero_nb);
311: VecSet(ipmP->Zero_nb,0.0);
312: VecDuplicate(ipmP->s,&ipmP->One_nb);
313: VecSet(ipmP->One_nb,1.0);
314: VecDuplicate(ipmP->s,&ipmP->Inf_nb);
315: VecSet(ipmP->Inf_nb,PETSC_INFINITY);
317: PetscMalloc1(ipmP->nb,&cind);
318: PetscMalloc1(ipmP->mi,&ucind);
319: VecGetOwnershipRange(ipmP->s,&sstart,&send);
321: if (ipmP->mi > 0) {
322: VecGetOwnershipRange(tao->constraints_inequality,&ucstart,&ucend);
323: counter=0;
324: for (i=ucstart;i<ucend;i++) {
325: cind[counter++] = i;
326: }
327: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isuc);
328: ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isc);
329: VecScatterCreate(tao->constraints_inequality,isuc,ipmP->ci,isc,&ipmP->ci_scat);
331: ISDestroy(&isuc);
332: ISDestroy(&isc);
333: }
334: /* need to know how may xbound indices are on each process */
335: /* TODO better way */
336: if (ipmP->nxlb) {
337: ISAllGather(ipmP->isxl,&bigxl);
338: ISGetIndices(bigxl,&xli);
339: /* find offsets for this processor */
340: xl_offset = ipmP->mi;
341: for (i=0;i<ipmP->nxlb;i++) {
342: if (xli[i] < xstart) {
343: xl_offset++;
344: } else break;
345: }
346: ISRestoreIndices(bigxl,&xli);
348: ISGetIndices(ipmP->isxl,&xli);
349: ISGetLocalSize(ipmP->isxl,&nloc);
350: for (i=0;i<nloc;i++) {
351: xind[i] = xli[i];
352: cind[i] = xl_offset+i;
353: }
355: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
356: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
357: VecScatterCreate(tao->XL,isx,ipmP->ci,isc,&ipmP->xl_scat);
358: ISDestroy(&isx);
359: ISDestroy(&isc);
360: ISDestroy(&bigxl);
361: }
363: if (ipmP->nxub) {
364: ISAllGather(ipmP->isxu,&bigxu);
365: ISGetIndices(bigxu,&xui);
366: /* find offsets for this processor */
367: xu_offset = ipmP->mi + ipmP->nxlb;
368: for (i=0;i<ipmP->nxub;i++) {
369: if (xui[i] < xstart) {
370: xu_offset++;
371: } else break;
372: }
373: ISRestoreIndices(bigxu,&xui);
375: ISGetIndices(ipmP->isxu,&xui);
376: ISGetLocalSize(ipmP->isxu,&nloc);
377: for (i=0;i<nloc;i++) {
378: xind[i] = xui[i];
379: cind[i] = xu_offset+i;
380: }
382: ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
383: ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
384: VecScatterCreate(tao->XU,isx,ipmP->ci,isc,&ipmP->xu_scat);
385: ISDestroy(&isx);
386: ISDestroy(&isc);
387: ISDestroy(&bigxu);
388: }
389: }
390: VecCreate(comm,&ipmP->bigrhs);
391: VecGetType(tao->solution,&vtype);
392: VecSetType(ipmP->bigrhs,vtype);
393: VecSetSizes(ipmP->bigrhs,PETSC_DECIDE,bigsize);
394: VecSetFromOptions(ipmP->bigrhs);
395: VecDuplicate(ipmP->bigrhs,&ipmP->bigstep);
397: /* create scatters for step->x and x->rhs */
398: for (i=xstart;i<xend;i++) {
399: stepind[i-xstart] = i;
400: xind[i-xstart] = i;
401: }
402: ISCreateGeneral(comm,xend-xstart,stepind,PETSC_COPY_VALUES,&sis);
403: ISCreateGeneral(comm,xend-xstart,xind,PETSC_COPY_VALUES,&is1);
404: VecScatterCreate(ipmP->bigstep,sis,tao->solution,is1,&ipmP->step1);
405: VecScatterCreate(tao->solution,is1,ipmP->bigrhs,sis,&ipmP->rhs1);
406: ISDestroy(&sis);
407: ISDestroy(&is1);
409: if (ipmP->nb > 0) {
410: for (i=sstart;i<send;i++) {
411: stepind[i-sstart] = i+ipmP->n;
412: cind[i-sstart] = i;
413: }
414: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
415: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
416: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step2);
417: ISDestroy(&sis);
419: for (i=sstart;i<send;i++) {
420: stepind[i-sstart] = i+ipmP->n+ipmP->me;
421: cind[i-sstart] = i;
422: }
423: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
424: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs3);
425: ISDestroy(&sis);
426: ISDestroy(&is1);
427: }
429: if (ipmP->me > 0) {
430: VecGetOwnershipRange(tao->constraints_equality,&ucestart,&uceend);
431: for (i=ucestart;i<uceend;i++) {
432: stepind[i-ucestart] = i + ipmP->n+ipmP->nb;
433: uceind[i-ucestart] = i;
434: }
436: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
437: ISCreateGeneral(comm,uceend-ucestart,uceind,PETSC_COPY_VALUES,&is1);
438: VecScatterCreate(ipmP->bigstep,sis,tao->constraints_equality,is1,&ipmP->step3);
439: ISDestroy(&sis);
441: for (i=ucestart;i<uceend;i++) {
442: stepind[i-ucestart] = i + ipmP->n;
443: }
445: ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
446: VecScatterCreate(tao->constraints_equality,is1,ipmP->bigrhs,sis,&ipmP->rhs2);
447: ISDestroy(&sis);
448: ISDestroy(&is1);
449: }
451: if (ipmP->nb > 0) {
452: for (i=sstart;i<send;i++) {
453: stepind[i-sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
454: cind[i-sstart] = i;
455: }
456: ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
457: ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
458: VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step4);
459: VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs4);
460: ISDestroy(&sis);
461: ISDestroy(&is1);
462: }
464: PetscFree(stepind);
465: PetscFree(cind);
466: PetscFree(ucind);
467: PetscFree(uceind);
468: PetscFree(xind);
469: return(0);
470: }
474: static PetscErrorCode TaoDestroy_IPM(Tao tao)
475: {
476: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
480: VecDestroy(&ipmP->rd);
481: VecDestroy(&ipmP->rpe);
482: VecDestroy(&ipmP->rpi);
483: VecDestroy(&ipmP->work);
484: VecDestroy(&ipmP->lamdae);
485: VecDestroy(&ipmP->lamdai);
486: VecDestroy(&ipmP->s);
487: VecDestroy(&ipmP->ds);
488: VecDestroy(&ipmP->ci);
490: VecDestroy(&ipmP->rhs_x);
491: VecDestroy(&ipmP->rhs_lamdae);
492: VecDestroy(&ipmP->rhs_lamdai);
493: VecDestroy(&ipmP->rhs_s);
495: VecDestroy(&ipmP->save_x);
496: VecDestroy(&ipmP->save_lamdae);
497: VecDestroy(&ipmP->save_lamdai);
498: VecDestroy(&ipmP->save_s);
500: VecScatterDestroy(&ipmP->step1);
501: VecScatterDestroy(&ipmP->step2);
502: VecScatterDestroy(&ipmP->step3);
503: VecScatterDestroy(&ipmP->step4);
505: VecScatterDestroy(&ipmP->rhs1);
506: VecScatterDestroy(&ipmP->rhs2);
507: VecScatterDestroy(&ipmP->rhs3);
508: VecScatterDestroy(&ipmP->rhs4);
510: VecScatterDestroy(&ipmP->ci_scat);
511: VecScatterDestroy(&ipmP->xl_scat);
512: VecScatterDestroy(&ipmP->xu_scat);
514: VecDestroy(&ipmP->dlamdai);
515: VecDestroy(&ipmP->dlamdae);
516: VecDestroy(&ipmP->Zero_nb);
517: VecDestroy(&ipmP->One_nb);
518: VecDestroy(&ipmP->Inf_nb);
519: VecDestroy(&ipmP->complementarity);
521: VecDestroy(&ipmP->bigrhs);
522: VecDestroy(&ipmP->bigstep);
523: MatDestroy(&ipmP->Ai);
524: MatDestroy(&ipmP->K);
525: ISDestroy(&ipmP->isxu);
526: ISDestroy(&ipmP->isxl);
527: PetscFree(tao->data);
528: return(0);
529: }
533: static PetscErrorCode TaoSetFromOptions_IPM(PetscOptionItems *PetscOptionsObject,Tao tao)
534: {
535: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
539: PetscOptionsHead(PetscOptionsObject,"IPM method for constrained optimization");
540: PetscOptionsBool("-tao_ipm_monitorkkt","monitor kkt status",NULL,ipmP->monitorkkt,&ipmP->monitorkkt,NULL);
541: PetscOptionsReal("-tao_ipm_pushs","parameter to push initial slack variables away from bounds",NULL,ipmP->pushs,&ipmP->pushs,NULL);
542: PetscOptionsReal("-tao_ipm_pushnu","parameter to push initial (inequality) dual variables away from bounds",NULL,ipmP->pushnu,&ipmP->pushnu,NULL);
543: PetscOptionsTail();
544: KSPSetFromOptions(tao->ksp);
545: return(0);
546: }
550: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
551: {
552: return 0;
553: }
555: /* IPMObjectiveAndGradient()
556: f = d'x + 0.5 * x' * H * x
557: rd = H*x + d + Ae'*lame - Ai'*lami
558: rpe = Ae*x - be
559: rpi = Ai*x - yi - bi
560: mu = yi' * lami/mi;
561: com = yi.*lami
563: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
564: */
565: /*
568: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
569: {
570: Tao tao = (Tao)tptr;
571: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
574: IPMComputeKKT(tao);
575: *f = ipmP->phi;
576: return(0);
577: }
578: */
580: /*
581: f = d'x + 0.5 * x' * H * x
582: rd = H*x + d + Ae'*lame - Ai'*lami
583: Ai = jac_ineq
584: I (w/lb)
585: -I (w/ub)
587: rpe = ce
588: rpi = ci - s;
589: com = s.*lami
590: mu = yi' * lami/mi;
592: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
593: */
596: static PetscErrorCode IPMComputeKKT(Tao tao)
597: {
598: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
599: PetscScalar norm;
603: VecCopy(tao->gradient,ipmP->rd);
605: if (ipmP->me > 0) {
606: /* rd = gradient + Ae'*lamdae */
607: MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
608: VecAXPY(ipmP->rd, 1.0, ipmP->work);
610: /* rpe = ce(x) */
611: VecCopy(tao->constraints_equality,ipmP->rpe);
612: }
613: if (ipmP->nb > 0) {
614: /* rd = rd - Ai'*lamdai */
615: MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
616: VecAXPY(ipmP->rd, -1.0, ipmP->work);
618: /* rpi = cin - s */
619: VecCopy(ipmP->ci,ipmP->rpi);
620: VecAXPY(ipmP->rpi, -1.0, ipmP->s);
622: /* com = s .* lami */
623: VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
624: }
625: /* phi = ||rd; rpe; rpi; com|| */
626: VecDot(ipmP->rd,ipmP->rd,&norm);
627: ipmP->phi = norm;
628: if (ipmP->me > 0 ) {
629: VecDot(ipmP->rpe,ipmP->rpe,&norm);
630: ipmP->phi += norm;
631: }
632: if (ipmP->nb > 0) {
633: VecDot(ipmP->rpi,ipmP->rpi,&norm);
634: ipmP->phi += norm;
635: VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
636: ipmP->phi += norm;
637: /* mu = s'*lami/nb */
638: VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
639: ipmP->mu /= ipmP->nb;
640: } else {
641: ipmP->mu = 1.0;
642: }
644: ipmP->phi = PetscSqrtScalar(ipmP->phi);
645: return(0);
646: }
650: /* evaluate user info at current point */
651: PetscErrorCode IPMEvaluate(Tao tao)
652: {
653: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
657: TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
658: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
659: if (ipmP->me > 0) {
660: TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
661: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
662: }
663: if (ipmP->mi > 0) {
664: TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
665: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
666: }
667: if (ipmP->nb > 0) {
668: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
669: IPMUpdateAi(tao);
670: }
671: return(0);
672: }
676: /* Push initial point away from bounds */
677: PetscErrorCode IPMPushInitialPoint(Tao tao)
678: {
679: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
683: TaoComputeVariableBounds(tao);
684: if (tao->XL && tao->XU) {
685: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
686: }
687: if (ipmP->nb > 0) {
688: VecSet(ipmP->s,ipmP->pushs);
689: VecSet(ipmP->lamdai,ipmP->pushnu);
690: if (ipmP->mi > 0) {
691: VecSet(tao->DI,ipmP->pushnu);
692: }
693: }
694: if (ipmP->me > 0) {
695: VecSet(tao->DE,1.0);
696: VecSet(ipmP->lamdae,1.0);
697: }
698: return(0);
699: }
703: PetscErrorCode IPMUpdateAi(Tao tao)
704: {
705: /* Ai = Ji
706: I (w/lb)
707: -I (w/ub) */
709: /* Ci = user->ci
710: Xi - lb (w/lb)
711: -Xi + ub (w/ub) */
713: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
714: MPI_Comm comm;
715: PetscInt i;
716: PetscScalar newval;
717: PetscInt newrow,newcol,ncols;
718: const PetscScalar *vals;
719: const PetscInt *cols;
720: PetscInt astart,aend,jstart,jend;
721: PetscInt *nonzeros;
722: PetscInt r2,r3,r4;
723: PetscMPIInt size;
724: PetscErrorCode ierr;
727: r2 = ipmP->mi;
728: r3 = r2 + ipmP->nxlb;
729: r4 = r3 + ipmP->nxub;
731: if (!ipmP->nb) return(0);
733: /* Create Ai matrix if it doesn't exist yet */
734: if (!ipmP->Ai) {
735: comm = ((PetscObject)(tao->solution))->comm;
736: PetscMalloc1(ipmP->nb,&nonzeros);
737: MPI_Comm_size(comm,&size);
738: if (size == 1) {
739: for (i=0;i<ipmP->mi;i++) {
740: MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
741: nonzeros[i] = ncols;
742: MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
743: }
744: for (i=r2;i<r4;i++) {
745: nonzeros[i] = 1;
746: }
747: }
748: MatCreate(comm,&ipmP->Ai);
749: MatSetType(ipmP->Ai,MATAIJ);
750: MatSetSizes(ipmP->Ai,PETSC_DECIDE,PETSC_DECIDE,ipmP->nb,ipmP->n);
751: MatSetFromOptions(ipmP->Ai);
752: MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
753: MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
754: if (size ==1) {
755: PetscFree(nonzeros);
756: }
757: }
759: /* Copy values from user jacobian to Ai */
760: MatGetOwnershipRange(ipmP->Ai,&astart,&aend);
762: /* Ai w/lb */
763: if (ipmP->mi) {
764: MatZeroEntries(ipmP->Ai);
765: MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
766: for (i=jstart;i<jend;i++) {
767: MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
768: newrow = i;
769: MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
770: MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
771: }
772: }
774: /* I w/ xlb */
775: if (ipmP->nxlb) {
776: for (i=0;i<ipmP->nxlb;i++) {
777: if (i>=astart && i<aend) {
778: newrow = i+r2;
779: newcol = i;
780: newval = 1.0;
781: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
782: }
783: }
784: }
785: if (ipmP->nxub) {
786: /* I w/ xub */
787: for (i=0;i<ipmP->nxub;i++) {
788: if (i>=astart && i<aend) {
789: newrow = i+r3;
790: newcol = i;
791: newval = -1.0;
792: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
793: }
794: }
795: }
797: MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
798: MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
799: CHKMEMQ;
801: VecSet(ipmP->ci,0.0);
803: /* user ci */
804: if (ipmP->mi > 0) {
805: VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
806: VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
807: }
808: if (!ipmP->work){
809: VecDuplicate(tao->solution,&ipmP->work);
810: }
811: VecCopy(tao->solution,ipmP->work);
812: if (tao->XL) {
813: VecAXPY(ipmP->work,-1.0,tao->XL);
815: /* lower bounds on variables */
816: if (ipmP->nxlb > 0) {
817: VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
818: VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
819: }
820: }
821: if (tao->XU) {
822: /* upper bounds on variables */
823: VecCopy(tao->solution,ipmP->work);
824: VecScale(ipmP->work,-1.0);
825: VecAXPY(ipmP->work,1.0,tao->XU);
826: if (ipmP->nxub > 0) {
827: VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
828: VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
829: }
830: }
831: return(0);
832: }
836: /* create K = [ Hlag , 0 , Ae', -Ai'];
837: [Ae , 0, 0 , 0];
838: [Ai ,-I, 0 , 0];
839: [ 0 , S , 0, Y ]; */
840: PetscErrorCode IPMUpdateK(Tao tao)
841: {
842: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
843: MPI_Comm comm;
844: PetscMPIInt size;
845: PetscErrorCode ierr;
846: PetscInt i,j,row;
847: PetscInt ncols,newcol,newcols[2],newrow;
848: const PetscInt *cols;
849: const PetscReal *vals;
850: const PetscReal *l,*y;
851: PetscReal *newvals;
852: PetscReal newval;
853: PetscInt subsize;
854: const PetscInt *indices;
855: PetscInt *nonzeros,*d_nonzeros,*o_nonzeros;
856: PetscInt bigsize;
857: PetscInt r1,r2,r3;
858: PetscInt c1,c2,c3;
859: PetscInt klocalsize;
860: PetscInt hstart,hend,kstart,kend;
861: PetscInt aistart,aiend,aestart,aeend;
862: PetscInt sstart,send;
865: comm = ((PetscObject)(tao->solution))->comm;
866: MPI_Comm_size(comm,&size);
867: IPMUpdateAi(tao);
869: /* allocate workspace */
870: subsize = PetscMax(ipmP->n,ipmP->nb);
871: subsize = PetscMax(ipmP->me,subsize);
872: subsize = PetscMax(2,subsize);
873: PetscMalloc1(subsize,&indices);
874: PetscMalloc1(subsize,&newvals);
876: r1 = c1 = ipmP->n;
877: r2 = r1 + ipmP->me; c2 = c1 + ipmP->nb;
878: r3 = c3 = r2 + ipmP->nb;
880: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
881: VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
882: MatGetOwnershipRange(tao->hessian,&hstart,&hend);
883: klocalsize = kend-kstart;
884: if (!ipmP->K) {
885: if (size == 1) {
886: PetscMalloc1(kend-kstart,&nonzeros);
887: for (i=0;i<bigsize;i++) {
888: if (i<r1) {
889: MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
890: nonzeros[i] = ncols;
891: MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
892: nonzeros[i] += ipmP->me+ipmP->nb;
893: } else if (i<r2) {
894: nonzeros[i-kstart] = ipmP->n;
895: } else if (i<r3) {
896: nonzeros[i-kstart] = ipmP->n+1;
897: } else if (i<bigsize) {
898: nonzeros[i-kstart] = 2;
899: }
900: }
901: MatCreate(comm,&ipmP->K);
902: MatSetType(ipmP->K,MATSEQAIJ);
903: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
904: MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
905: MatSetFromOptions(ipmP->K);
906: PetscFree(nonzeros);
907: } else {
908: PetscMalloc1(kend-kstart,&d_nonzeros);
909: PetscMalloc1(kend-kstart,&o_nonzeros);
910: for (i=kstart;i<kend;i++) {
911: if (i<r1) {
912: /* TODO fix preallocation for mpi mats */
913: d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
914: o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
915: } else if (i<r2) {
916: d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
917: o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
918: } else if (i<r3) {
919: d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
920: o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
921: } else {
922: d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
923: o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
924: }
925: }
926: MatCreate(comm,&ipmP->K);
927: MatSetType(ipmP->K,MATMPIAIJ);
928: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
929: MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
930: PetscFree(d_nonzeros);
931: PetscFree(o_nonzeros);
932: MatSetFromOptions(ipmP->K);
933: }
934: }
936: MatZeroEntries(ipmP->K);
937: /* Copy H */
938: for (i=hstart;i<hend;i++) {
939: MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
940: if (ncols > 0) {
941: MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
942: }
943: MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
944: }
946: /* Copy Ae and Ae' */
947: if (ipmP->me > 0) {
948: MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
949: for (i=aestart;i<aeend;i++) {
950: MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
951: if (ncols > 0) {
952: /*Ae*/
953: row = i+r1;
954: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
955: /*Ae'*/
956: for (j=0;j<ncols;j++) {
957: newcol = i + c2;
958: newrow = cols[j];
959: newval = vals[j];
960: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
961: }
962: }
963: MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
964: }
965: }
967: if (ipmP->nb > 0) {
968: MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
969: /* Copy Ai,and Ai' */
970: for (i=aistart;i<aiend;i++) {
971: row = i+r2;
972: MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
973: if (ncols > 0) {
974: /*Ai*/
975: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
976: /*-Ai'*/
977: for (j=0;j<ncols;j++) {
978: newcol = i + c3;
979: newrow = cols[j];
980: newval = -vals[j];
981: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
982: }
983: }
984: MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
985: }
987: /* -I */
988: for (i=kstart;i<kend;i++) {
989: if (i>=r2 && i<r3) {
990: newrow = i;
991: newcol = i-r2+c1;
992: newval = -1.0;
993: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
994: }
995: }
997: /* Copy L,Y */
998: VecGetOwnershipRange(ipmP->s,&sstart,&send);
999: VecGetArrayRead(ipmP->lamdai,&l);
1000: VecGetArrayRead(ipmP->s,&y);
1002: for (i=sstart;i<send;i++) {
1003: newcols[0] = c1+i;
1004: newcols[1] = c3+i;
1005: newvals[0] = l[i-sstart];
1006: newvals[1] = y[i-sstart];
1007: newrow = r3+i;
1008: MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
1009: }
1011: VecRestoreArrayRead(ipmP->lamdai,&l);
1012: VecRestoreArrayRead(ipmP->s,&y);
1013: }
1015: PetscFree(indices);
1016: PetscFree(newvals);
1017: MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
1018: MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
1019: return(0);
1020: }
1024: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
1025: {
1026: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1030: /* rhs = [x1 (n)
1031: x2 (me)
1032: x3 (nb)
1033: x4 (nb)] */
1034: if (X1) {
1035: VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1036: VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1037: }
1038: if (ipmP->me > 0 && X2) {
1039: VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1040: VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1041: }
1042: if (ipmP->nb > 0) {
1043: if (X3) {
1044: VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1045: VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1046: }
1047: if (X4) {
1048: VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1049: VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1050: }
1051: }
1052: return(0);
1053: }
1057: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1058: {
1059: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1063: CHKMEMQ;
1064: /* [x1 (n)
1065: x2 (nb) may be 0
1066: x3 (me) may be 0
1067: x4 (nb) may be 0 */
1068: if (X1) {
1069: VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1070: VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1071: }
1072: if (X2 && ipmP->nb > 0) {
1073: VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1074: VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1075: }
1076: if (X3 && ipmP->me > 0) {
1077: VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1078: VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1079: }
1080: if (X4 && ipmP->nb > 0) {
1081: VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1082: VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1083: }
1084: CHKMEMQ;
1085: return(0);
1086: }
1088: /*MC
1089: TAOIPM - Interior point algorithm for generally constrained optimization.
1091: Option Database Keys:
1092: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1093: . -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1095: 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.
1096: Level: beginner
1098: M*/
1102: PETSC_EXTERN 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;
1118: /* Override default settings (unless already changed) */
1119: if (!tao->max_it_changed) tao->max_it = 200;
1120: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1122: ipmP->dec = 10000; /* line search critera */
1123: ipmP->taumin = 0.995;
1124: ipmP->monitorkkt = PETSC_FALSE;
1125: ipmP->pushs = 100;
1126: ipmP->pushnu = 100;
1127: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1128: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1129: return(0);
1130: }