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