Actual source code: ipm.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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:     /* Call general purpose update function */
 56:     if (tao->ops->update) {
 57:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 58:     }
 59:     
 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=0;
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=0;
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: }