Actual source code: ipm.c
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: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
36: PetscInt its,i;
37: PetscScalar stepsize=1.0;
38: PetscScalar step_s,step_l,alpha,tau,sigma,phi_target;
40: /* Push initial point away from bounds */
41: IPMInitializeBounds(tao);
42: IPMPushInitialPoint(tao);
43: VecCopy(tao->solution,ipmP->rhs_x);
44: IPMEvaluate(tao);
45: IPMComputeKKT(tao);
47: tao->reason = TAO_CONTINUE_ITERATING;
48: TaoLogConvergenceHistory(tao,ipmP->kkt_f,ipmP->phi,0.0,tao->ksp_its);
49: TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,1.0);
50: (*tao->ops->convergencetest)(tao,tao->cnvP);
52: while (tao->reason == TAO_CONTINUE_ITERATING) {
53: /* Call general purpose update function */
54: if (tao->ops->update) {
55: (*tao->ops->update)(tao, tao->niter, tao->user_update);
56: }
58: tao->ksp_its=0;
59: IPMUpdateK(tao);
60: /*
61: rhs.x = -rd
62: rhs.lame = -rpe
63: rhs.lami = -rpi
64: rhs.com = -com
65: */
67: VecCopy(ipmP->rd,ipmP->rhs_x);
68: if (ipmP->me > 0) {
69: VecCopy(ipmP->rpe,ipmP->rhs_lamdae);
70: }
71: if (ipmP->nb > 0) {
72: VecCopy(ipmP->rpi,ipmP->rhs_lamdai);
73: VecCopy(ipmP->complementarity,ipmP->rhs_s);
74: }
75: IPMGatherRHS(tao,ipmP->bigrhs,ipmP->rhs_x,ipmP->rhs_lamdae,ipmP->rhs_lamdai,ipmP->rhs_s);
76: VecScale(ipmP->bigrhs,-1.0);
78: /* solve K * step = rhs */
79: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
80: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
82: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
83: KSPGetIterationNumber(tao->ksp,&its);
84: tao->ksp_its += its;
85: tao->ksp_tot_its+=its;
86: /* Find distance along step direction to closest bound */
87: if (ipmP->nb > 0) {
88: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
89: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
90: alpha = PetscMin(step_s,step_l);
91: alpha = PetscMin(alpha,1.0);
92: ipmP->alpha1 = alpha;
93: } else {
94: ipmP->alpha1 = alpha = 1.0;
95: }
97: /* x_aff = x + alpha*d */
98: VecCopy(tao->solution,ipmP->save_x);
99: if (ipmP->me > 0) {
100: VecCopy(ipmP->lamdae,ipmP->save_lamdae);
101: }
102: if (ipmP->nb > 0) {
103: VecCopy(ipmP->lamdai,ipmP->save_lamdai);
104: VecCopy(ipmP->s,ipmP->save_s);
105: }
107: VecAXPY(tao->solution,alpha,tao->stepdirection);
108: if (ipmP->me > 0) {
109: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
110: }
111: if (ipmP->nb > 0) {
112: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
113: VecAXPY(ipmP->s,alpha,ipmP->ds);
114: }
116: /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */
117: if (ipmP->mu == 0.0) {
118: sigma = 0.0;
119: } else {
120: sigma = 1.0/ipmP->mu;
121: }
122: IPMComputeKKT(tao);
123: sigma *= ipmP->mu;
124: sigma*=sigma*sigma;
126: /* revert kkt info */
127: VecCopy(ipmP->save_x,tao->solution);
128: if (ipmP->me > 0) {
129: VecCopy(ipmP->save_lamdae,ipmP->lamdae);
130: }
131: if (ipmP->nb > 0) {
132: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
133: VecCopy(ipmP->save_s,ipmP->s);
134: }
135: IPMComputeKKT(tao);
137: /* update rhs with new complementarity vector */
138: if (ipmP->nb > 0) {
139: VecCopy(ipmP->complementarity,ipmP->rhs_s);
140: VecScale(ipmP->rhs_s,-1.0);
141: VecShift(ipmP->rhs_s,sigma*ipmP->mu);
142: }
143: IPMGatherRHS(tao,ipmP->bigrhs,NULL,NULL,NULL,ipmP->rhs_s);
145: /* solve K * step = rhs */
146: KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
147: KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);
149: IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
150: KSPGetIterationNumber(tao->ksp,&its);
151: tao->ksp_its += its;
152: tao->ksp_tot_its+=its;
153: if (ipmP->nb > 0) {
154: /* Get max step size and apply frac-to-boundary */
155: tau = PetscMax(ipmP->taumin,1.0-ipmP->mu);
156: tau = PetscMin(tau,1.0);
157: if (tau != 1.0) {
158: VecScale(ipmP->s,tau);
159: VecScale(ipmP->lamdai,tau);
160: }
161: VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
162: VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
163: if (tau != 1.0) {
164: VecCopy(ipmP->save_s,ipmP->s);
165: VecCopy(ipmP->save_lamdai,ipmP->lamdai);
166: }
167: alpha = PetscMin(step_s,step_l);
168: alpha = PetscMin(alpha,1.0);
169: } else {
170: alpha = 1.0;
171: }
172: ipmP->alpha2 = alpha;
173: /* TODO make phi_target meaningful */
174: phi_target = ipmP->dec * ipmP->phi;
175: for (i=0; i<11;i++) {
176: VecAXPY(tao->solution,alpha,tao->stepdirection);
177: if (ipmP->nb > 0) {
178: VecAXPY(ipmP->s,alpha,ipmP->ds);
179: VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
180: }
181: if (ipmP->me > 0) {
182: VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
183: }
185: /* update dual variables */
186: if (ipmP->me > 0) {
187: VecCopy(ipmP->lamdae,tao->DE);
188: }
190: IPMEvaluate(tao);
191: IPMComputeKKT(tao);
192: if (ipmP->phi <= phi_target) break;
193: alpha /= 2.0;
194: }
196: TaoLogConvergenceHistory(tao,ipmP->kkt_f,ipmP->phi,0.0,tao->ksp_its);
197: TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,stepsize);
198: (*tao->ops->convergencetest)(tao,tao->cnvP);
199: tao->niter++;
200: }
201: return 0;
202: }
204: static PetscErrorCode TaoSetup_IPM(Tao tao)
205: {
206: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
208: ipmP->nb = ipmP->mi = ipmP->me = 0;
209: ipmP->K = NULL;
210: VecGetSize(tao->solution,&ipmP->n);
211: if (!tao->gradient) {
212: VecDuplicate(tao->solution, &tao->gradient);
213: VecDuplicate(tao->solution, &tao->stepdirection);
214: VecDuplicate(tao->solution, &ipmP->rd);
215: VecDuplicate(tao->solution, &ipmP->rhs_x);
216: VecDuplicate(tao->solution, &ipmP->work);
217: VecDuplicate(tao->solution, &ipmP->save_x);
218: }
219: if (tao->constraints_equality) {
220: VecGetSize(tao->constraints_equality,&ipmP->me);
221: VecDuplicate(tao->constraints_equality,&ipmP->lamdae);
222: VecDuplicate(tao->constraints_equality,&ipmP->dlamdae);
223: VecDuplicate(tao->constraints_equality,&ipmP->rhs_lamdae);
224: VecDuplicate(tao->constraints_equality,&ipmP->save_lamdae);
225: VecDuplicate(tao->constraints_equality,&ipmP->rpe);
226: VecDuplicate(tao->constraints_equality,&tao->DE);
227: }
228: if (tao->constraints_inequality) {
229: VecDuplicate(tao->constraints_inequality,&tao->DI);
230: }
231: return 0;
232: }
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 = 0 ,send = 0;
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;
249: MPI_Comm comm;
251: cind=xind=ucind=uceind=stepind=NULL;
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: PetscObjectGetComm((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: }
470: static PetscErrorCode TaoDestroy_IPM(Tao tao)
471: {
472: TAO_IPM *ipmP = (TAO_IPM*)tao->data;
474: VecDestroy(&ipmP->rd);
475: VecDestroy(&ipmP->rpe);
476: VecDestroy(&ipmP->rpi);
477: VecDestroy(&ipmP->work);
478: VecDestroy(&ipmP->lamdae);
479: VecDestroy(&ipmP->lamdai);
480: VecDestroy(&ipmP->s);
481: VecDestroy(&ipmP->ds);
482: VecDestroy(&ipmP->ci);
484: VecDestroy(&ipmP->rhs_x);
485: VecDestroy(&ipmP->rhs_lamdae);
486: VecDestroy(&ipmP->rhs_lamdai);
487: VecDestroy(&ipmP->rhs_s);
489: VecDestroy(&ipmP->save_x);
490: VecDestroy(&ipmP->save_lamdae);
491: VecDestroy(&ipmP->save_lamdai);
492: VecDestroy(&ipmP->save_s);
494: VecScatterDestroy(&ipmP->step1);
495: VecScatterDestroy(&ipmP->step2);
496: VecScatterDestroy(&ipmP->step3);
497: VecScatterDestroy(&ipmP->step4);
499: VecScatterDestroy(&ipmP->rhs1);
500: VecScatterDestroy(&ipmP->rhs2);
501: VecScatterDestroy(&ipmP->rhs3);
502: VecScatterDestroy(&ipmP->rhs4);
504: VecScatterDestroy(&ipmP->ci_scat);
505: VecScatterDestroy(&ipmP->xl_scat);
506: VecScatterDestroy(&ipmP->xu_scat);
508: VecDestroy(&ipmP->dlamdai);
509: VecDestroy(&ipmP->dlamdae);
510: VecDestroy(&ipmP->Zero_nb);
511: VecDestroy(&ipmP->One_nb);
512: VecDestroy(&ipmP->Inf_nb);
513: VecDestroy(&ipmP->complementarity);
515: VecDestroy(&ipmP->bigrhs);
516: VecDestroy(&ipmP->bigstep);
517: MatDestroy(&ipmP->Ai);
518: MatDestroy(&ipmP->K);
519: ISDestroy(&ipmP->isxu);
520: ISDestroy(&ipmP->isxl);
521: PetscFree(tao->data);
522: return 0;
523: }
525: static PetscErrorCode TaoSetFromOptions_IPM(PetscOptionItems *PetscOptionsObject,Tao tao)
526: {
527: 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;
558: IPMComputeKKT(tao);
559: *f = ipmP->phi;
560: return 0;
561: }
562: */
564: /*
565: f = d'x + 0.5 * x' * H * x
566: rd = H*x + d + Ae'*lame - Ai'*lami
567: Ai = jac_ineq
568: I (w/lb)
569: -I (w/ub)
571: rpe = ce
572: rpi = ci - s;
573: com = s.*lami
574: mu = yi' * lami/mi;
576: phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
577: */
578: static PetscErrorCode IPMComputeKKT(Tao tao)
579: {
580: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
581: PetscScalar norm;
583: VecCopy(tao->gradient,ipmP->rd);
585: if (ipmP->me > 0) {
586: /* rd = gradient + Ae'*lamdae */
587: MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
588: VecAXPY(ipmP->rd, 1.0, ipmP->work);
590: /* rpe = ce(x) */
591: VecCopy(tao->constraints_equality,ipmP->rpe);
592: }
593: if (ipmP->nb > 0) {
594: /* rd = rd - Ai'*lamdai */
595: MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
596: VecAXPY(ipmP->rd, -1.0, ipmP->work);
598: /* rpi = cin - s */
599: VecCopy(ipmP->ci,ipmP->rpi);
600: VecAXPY(ipmP->rpi, -1.0, ipmP->s);
602: /* com = s .* lami */
603: VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
604: }
605: /* phi = ||rd; rpe; rpi; com|| */
606: VecDot(ipmP->rd,ipmP->rd,&norm);
607: ipmP->phi = norm;
608: if (ipmP->me > 0) {
609: VecDot(ipmP->rpe,ipmP->rpe,&norm);
610: ipmP->phi += norm;
611: }
612: if (ipmP->nb > 0) {
613: VecDot(ipmP->rpi,ipmP->rpi,&norm);
614: ipmP->phi += norm;
615: VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
616: ipmP->phi += norm;
617: /* mu = s'*lami/nb */
618: VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
619: ipmP->mu /= ipmP->nb;
620: } else {
621: ipmP->mu = 1.0;
622: }
624: ipmP->phi = PetscSqrtScalar(ipmP->phi);
625: return 0;
626: }
628: /* evaluate user info at current point */
629: PetscErrorCode IPMEvaluate(Tao tao)
630: {
631: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
633: TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
634: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
635: if (ipmP->me > 0) {
636: TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
637: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
638: }
639: if (ipmP->mi > 0) {
640: TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
641: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
642: }
643: if (ipmP->nb > 0) {
644: /* Ai' = jac_ineq | I (w/lb) | -I (w/ub) */
645: IPMUpdateAi(tao);
646: }
647: return 0;
648: }
650: /* Push initial point away from bounds */
651: PetscErrorCode IPMPushInitialPoint(Tao tao)
652: {
653: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
655: TaoComputeVariableBounds(tao);
656: if (tao->XL && tao->XU) {
657: VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
658: }
659: if (ipmP->nb > 0) {
660: VecSet(ipmP->s,ipmP->pushs);
661: VecSet(ipmP->lamdai,ipmP->pushnu);
662: if (ipmP->mi > 0) {
663: VecSet(tao->DI,ipmP->pushnu);
664: }
665: }
666: if (ipmP->me > 0) {
667: VecSet(tao->DE,1.0);
668: VecSet(ipmP->lamdae,1.0);
669: }
670: return 0;
671: }
673: PetscErrorCode IPMUpdateAi(Tao tao)
674: {
675: /* Ai = Ji
676: I (w/lb)
677: -I (w/ub) */
679: /* Ci = user->ci
680: Xi - lb (w/lb)
681: -Xi + ub (w/ub) */
683: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
684: MPI_Comm comm;
685: PetscInt i;
686: PetscScalar newval;
687: PetscInt newrow,newcol,ncols;
688: const PetscScalar *vals;
689: const PetscInt *cols;
690: PetscInt astart,aend,jstart,jend;
691: PetscInt *nonzeros;
692: PetscInt r2,r3,r4;
693: PetscMPIInt size;
694: Vec solu;
695: PetscInt nloc;
697: r2 = ipmP->mi;
698: r3 = r2 + ipmP->nxlb;
699: r4 = r3 + ipmP->nxub;
701: if (!ipmP->nb) return 0;
703: /* Create Ai matrix if it doesn't exist yet */
704: if (!ipmP->Ai) {
705: comm = ((PetscObject)(tao->solution))->comm;
706: MPI_Comm_size(comm,&size);
707: if (size == 1) {
708: PetscMalloc1(ipmP->nb,&nonzeros);
709: for (i=0;i<ipmP->mi;i++) {
710: MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
711: nonzeros[i] = ncols;
712: MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
713: }
714: for (i=r2;i<r4;i++) {
715: nonzeros[i] = 1;
716: }
717: }
718: MatCreate(comm,&ipmP->Ai);
719: MatSetType(ipmP->Ai,MATAIJ);
721: TaoGetSolution(tao,&solu);
722: VecGetLocalSize(solu,&nloc);
723: MatSetSizes(ipmP->Ai,PETSC_DECIDE,nloc,ipmP->nb,PETSC_DECIDE);
724: MatSetFromOptions(ipmP->Ai);
725: MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
726: MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
727: if (size ==1) {
728: PetscFree(nonzeros);
729: }
730: }
732: /* Copy values from user jacobian to Ai */
733: MatGetOwnershipRange(ipmP->Ai,&astart,&aend);
735: /* Ai w/lb */
736: if (ipmP->mi) {
737: MatZeroEntries(ipmP->Ai);
738: MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
739: for (i=jstart;i<jend;i++) {
740: MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
741: newrow = i;
742: MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
743: MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
744: }
745: }
747: /* I w/ xlb */
748: if (ipmP->nxlb) {
749: for (i=0;i<ipmP->nxlb;i++) {
750: if (i>=astart && i<aend) {
751: newrow = i+r2;
752: newcol = i;
753: newval = 1.0;
754: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
755: }
756: }
757: }
758: if (ipmP->nxub) {
759: /* I w/ xub */
760: for (i=0;i<ipmP->nxub;i++) {
761: if (i>=astart && i<aend) {
762: newrow = i+r3;
763: newcol = i;
764: newval = -1.0;
765: MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
766: }
767: }
768: }
770: MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
771: MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
772: CHKMEMQ;
774: VecSet(ipmP->ci,0.0);
776: /* user ci */
777: if (ipmP->mi > 0) {
778: VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
779: VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
780: }
781: if (!ipmP->work) {
782: VecDuplicate(tao->solution,&ipmP->work);
783: }
784: VecCopy(tao->solution,ipmP->work);
785: if (tao->XL) {
786: VecAXPY(ipmP->work,-1.0,tao->XL);
788: /* lower bounds on variables */
789: if (ipmP->nxlb > 0) {
790: VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
791: VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
792: }
793: }
794: if (tao->XU) {
795: /* upper bounds on variables */
796: VecCopy(tao->solution,ipmP->work);
797: VecScale(ipmP->work,-1.0);
798: VecAXPY(ipmP->work,1.0,tao->XU);
799: if (ipmP->nxub > 0) {
800: VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
801: VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
802: }
803: }
804: return 0;
805: }
807: /* create K = [ Hlag , 0 , Ae', -Ai'];
808: [Ae , 0, 0 , 0];
809: [Ai ,-I, 0 , 0];
810: [ 0 , S , 0, Y ]; */
811: PetscErrorCode IPMUpdateK(Tao tao)
812: {
813: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
814: MPI_Comm comm;
815: PetscMPIInt size;
816: PetscInt i,j,row;
817: PetscInt ncols,newcol,newcols[2],newrow;
818: const PetscInt *cols;
819: const PetscReal *vals;
820: const PetscReal *l,*y;
821: PetscReal *newvals;
822: PetscReal newval;
823: PetscInt subsize;
824: const PetscInt *indices;
825: PetscInt *nonzeros,*d_nonzeros,*o_nonzeros;
826: PetscInt bigsize;
827: PetscInt r1,r2,r3;
828: PetscInt c1,c2,c3;
829: PetscInt klocalsize;
830: PetscInt hstart,hend,kstart,kend;
831: PetscInt aistart,aiend,aestart,aeend;
832: PetscInt sstart,send;
834: comm = ((PetscObject)(tao->solution))->comm;
835: MPI_Comm_size(comm,&size);
836: IPMUpdateAi(tao);
838: /* allocate workspace */
839: subsize = PetscMax(ipmP->n,ipmP->nb);
840: subsize = PetscMax(ipmP->me,subsize);
841: subsize = PetscMax(2,subsize);
842: PetscMalloc1(subsize,(PetscInt**)&indices);
843: PetscMalloc1(subsize,&newvals);
845: r1 = c1 = ipmP->n;
846: r2 = r1 + ipmP->me; c2 = c1 + ipmP->nb;
847: r3 = c3 = r2 + ipmP->nb;
849: bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
850: VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
851: MatGetOwnershipRange(tao->hessian,&hstart,&hend);
852: klocalsize = kend-kstart;
853: if (!ipmP->K) {
854: if (size == 1) {
855: PetscMalloc1(kend-kstart,&nonzeros);
856: for (i=0;i<bigsize;i++) {
857: if (i<r1) {
858: MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
859: nonzeros[i] = ncols;
860: MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
861: nonzeros[i] += ipmP->me+ipmP->nb;
862: } else if (i<r2) {
863: nonzeros[i-kstart] = ipmP->n;
864: } else if (i<r3) {
865: nonzeros[i-kstart] = ipmP->n+1;
866: } else if (i<bigsize) {
867: nonzeros[i-kstart] = 2;
868: }
869: }
870: MatCreate(comm,&ipmP->K);
871: MatSetType(ipmP->K,MATSEQAIJ);
872: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
873: MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
874: MatSetFromOptions(ipmP->K);
875: PetscFree(nonzeros);
876: } else {
877: PetscMalloc1(kend-kstart,&d_nonzeros);
878: PetscMalloc1(kend-kstart,&o_nonzeros);
879: for (i=kstart;i<kend;i++) {
880: if (i<r1) {
881: /* TODO fix preallocation for mpi mats */
882: d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
883: o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
884: } else if (i<r2) {
885: d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
886: o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
887: } else if (i<r3) {
888: d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
889: o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
890: } else {
891: d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
892: o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
893: }
894: }
895: MatCreate(comm,&ipmP->K);
896: MatSetType(ipmP->K,MATMPIAIJ);
897: MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
898: MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
899: PetscFree(d_nonzeros);
900: PetscFree(o_nonzeros);
901: MatSetFromOptions(ipmP->K);
902: }
903: }
905: MatZeroEntries(ipmP->K);
906: /* Copy H */
907: for (i=hstart;i<hend;i++) {
908: MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
909: if (ncols > 0) {
910: MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
911: }
912: MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
913: }
915: /* Copy Ae and Ae' */
916: if (ipmP->me > 0) {
917: MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
918: for (i=aestart;i<aeend;i++) {
919: MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
920: if (ncols > 0) {
921: /*Ae*/
922: row = i+r1;
923: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
924: /*Ae'*/
925: for (j=0;j<ncols;j++) {
926: newcol = i + c2;
927: newrow = cols[j];
928: newval = vals[j];
929: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
930: }
931: }
932: MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
933: }
934: }
936: if (ipmP->nb > 0) {
937: MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
938: /* Copy Ai,and Ai' */
939: for (i=aistart;i<aiend;i++) {
940: row = i+r2;
941: MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
942: if (ncols > 0) {
943: /*Ai*/
944: MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
945: /*-Ai'*/
946: for (j=0;j<ncols;j++) {
947: newcol = i + c3;
948: newrow = cols[j];
949: newval = -vals[j];
950: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
951: }
952: }
953: MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
954: }
956: /* -I */
957: for (i=kstart;i<kend;i++) {
958: if (i>=r2 && i<r3) {
959: newrow = i;
960: newcol = i-r2+c1;
961: newval = -1.0;
962: MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
963: }
964: }
966: /* Copy L,Y */
967: VecGetOwnershipRange(ipmP->s,&sstart,&send);
968: VecGetArrayRead(ipmP->lamdai,&l);
969: VecGetArrayRead(ipmP->s,&y);
971: for (i=sstart;i<send;i++) {
972: newcols[0] = c1+i;
973: newcols[1] = c3+i;
974: newvals[0] = l[i-sstart];
975: newvals[1] = y[i-sstart];
976: newrow = r3+i;
977: MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
978: }
980: VecRestoreArrayRead(ipmP->lamdai,&l);
981: VecRestoreArrayRead(ipmP->s,&y);
982: }
984: PetscFree(indices);
985: PetscFree(newvals);
986: MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
987: MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
988: return 0;
989: }
991: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
992: {
993: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
995: /* rhs = [x1 (n)
996: x2 (me)
997: x3 (nb)
998: x4 (nb)] */
999: if (X1) {
1000: VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1001: VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1002: }
1003: if (ipmP->me > 0 && X2) {
1004: VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1005: VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1006: }
1007: if (ipmP->nb > 0) {
1008: if (X3) {
1009: VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1010: VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1011: }
1012: if (X4) {
1013: VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1014: VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1015: }
1016: }
1017: return 0;
1018: }
1020: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1021: {
1022: TAO_IPM *ipmP = (TAO_IPM *)tao->data;
1024: CHKMEMQ;
1025: /* [x1 (n)
1026: x2 (nb) may be 0
1027: x3 (me) may be 0
1028: x4 (nb) may be 0 */
1029: if (X1) {
1030: VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1031: VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1032: }
1033: if (X2 && ipmP->nb > 0) {
1034: VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1035: VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1036: }
1037: if (X3 && ipmP->me > 0) {
1038: VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1039: VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1040: }
1041: if (X4 && ipmP->nb > 0) {
1042: VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1043: VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1044: }
1045: CHKMEMQ;
1046: return 0;
1047: }
1049: /*MC
1050: TAOIPM - Interior point algorithm for generally constrained optimization.
1052: Option Database Keys:
1053: + -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1054: - -tao_ipm_pushs - parameter to push initial slack variables away from bounds
1056: Notes:
1057: 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.
1058: Level: beginner
1060: M*/
1062: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1063: {
1064: TAO_IPM *ipmP;
1066: tao->ops->setup = TaoSetup_IPM;
1067: tao->ops->solve = TaoSolve_IPM;
1068: tao->ops->view = TaoView_IPM;
1069: tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1070: tao->ops->destroy = TaoDestroy_IPM;
1071: /* tao->ops->computedual = TaoComputeDual_IPM; */
1073: PetscNewLog(tao,&ipmP);
1074: tao->data = (void*)ipmP;
1076: /* Override default settings (unless already changed) */
1077: if (!tao->max_it_changed) tao->max_it = 200;
1078: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1080: ipmP->dec = 10000; /* line search critera */
1081: ipmP->taumin = 0.995;
1082: ipmP->monitorkkt = PETSC_FALSE;
1083: ipmP->pushs = 100;
1084: ipmP->pushnu = 100;
1085: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1086: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
1087: KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1088: return 0;
1089: }