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