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: }