Actual source code: ipm.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/constrained/impls/ipm/ipm.h> /*I "ipm.h" I*/

  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);

 35: static PetscErrorCode TaoSolve_IPM(Tao tao)
 36: {
 37:   PetscErrorCode     ierr;
 38:   TAO_IPM            *ipmP = (TAO_IPM*)tao->data;
 39:   TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
 40:   PetscInt           its,i;
 41:   PetscScalar        stepsize=1.0;
 42:   PetscScalar        step_s,step_l,alpha,tau,sigma,phi_target;

 45:   /* Push initial point away from bounds */
 46:   IPMInitializeBounds(tao);
 47:   IPMPushInitialPoint(tao);
 48:   VecCopy(tao->solution,ipmP->rhs_x);
 49:   IPMEvaluate(tao);
 50:   IPMComputeKKT(tao);
 51:   TaoMonitor(tao,tao->niter++,ipmP->kkt_f,ipmP->phi,0.0,1.0,&reason);

 53:   while (reason == TAO_CONTINUE_ITERATING) {
 54:     tao->ksp_its=0;
 55:     IPMUpdateK(tao);
 56:     /*
 57:        rhs.x    = -rd
 58:        rhs.lame = -rpe
 59:        rhs.lami = -rpi
 60:        rhs.com  = -com
 61:     */

 63:     VecCopy(ipmP->rd,ipmP->rhs_x);
 64:     if (ipmP->me > 0) {
 65:       VecCopy(ipmP->rpe,ipmP->rhs_lamdae);
 66:     }
 67:     if (ipmP->nb > 0) {
 68:       VecCopy(ipmP->rpi,ipmP->rhs_lamdai);
 69:       VecCopy(ipmP->complementarity,ipmP->rhs_s);
 70:     }
 71:     IPMGatherRHS(tao,ipmP->bigrhs,ipmP->rhs_x,ipmP->rhs_lamdae,ipmP->rhs_lamdai,ipmP->rhs_s);
 72:     VecScale(ipmP->bigrhs,-1.0);

 74:     /* solve K * step = rhs */
 75:     KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
 76:     KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);

 78:     IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
 79:     KSPGetIterationNumber(tao->ksp,&its);
 80:     tao->ksp_its += its;
 81:     tao->ksp_tot_its+=its;
 82:      /* Find distance along step direction to closest bound */
 83:     if (ipmP->nb > 0) {
 84:       VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
 85:       VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
 86:       alpha = PetscMin(step_s,step_l);
 87:       alpha = PetscMin(alpha,1.0);
 88:       ipmP->alpha1 = alpha;
 89:     } else {
 90:       ipmP->alpha1 = alpha = 1.0;
 91:     }

 93:     /* x_aff = x + alpha*d */
 94:     VecCopy(tao->solution,ipmP->save_x);
 95:     if (ipmP->me > 0) {
 96:       VecCopy(ipmP->lamdae,ipmP->save_lamdae);
 97:     }
 98:     if (ipmP->nb > 0) {
 99:       VecCopy(ipmP->lamdai,ipmP->save_lamdai);
100:       VecCopy(ipmP->s,ipmP->save_s);
101:     }

103:     VecAXPY(tao->solution,alpha,tao->stepdirection);
104:     if (ipmP->me > 0) {
105:       VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
106:     }
107:     if (ipmP->nb > 0) {
108:       VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
109:       VecAXPY(ipmP->s,alpha,ipmP->ds);
110:     }

112:     /* Recompute kkt to find centering parameter sigma = (new_mu/old_mu)^3 */
113:     if (ipmP->mu == 0.0) {
114:       sigma = 0.0;
115:     } else {
116:       sigma = 1.0/ipmP->mu;
117:     }
118:     IPMComputeKKT(tao);
119:     sigma *= ipmP->mu;
120:     sigma*=sigma*sigma;

122:     /* revert kkt info */
123:     VecCopy(ipmP->save_x,tao->solution);
124:     if (ipmP->me > 0) {
125:       VecCopy(ipmP->save_lamdae,ipmP->lamdae);
126:     }
127:     if (ipmP->nb > 0) {
128:       VecCopy(ipmP->save_lamdai,ipmP->lamdai);
129:       VecCopy(ipmP->save_s,ipmP->s);
130:     }
131:     IPMComputeKKT(tao);

133:     /* update rhs with new complementarity vector */
134:     if (ipmP->nb > 0) {
135:       VecCopy(ipmP->complementarity,ipmP->rhs_s);
136:       VecScale(ipmP->rhs_s,-1.0);
137:       VecShift(ipmP->rhs_s,sigma*ipmP->mu);
138:     }
139:     IPMGatherRHS(tao,ipmP->bigrhs,NULL,NULL,NULL,ipmP->rhs_s);

141:     /* solve K * step = rhs */
142:     KSPSetOperators(tao->ksp,ipmP->K,ipmP->K);
143:     KSPSolve(tao->ksp,ipmP->bigrhs,ipmP->bigstep);

145:     IPMScatterStep(tao,ipmP->bigstep,tao->stepdirection,ipmP->ds,ipmP->dlamdae,ipmP->dlamdai);
146:     KSPGetIterationNumber(tao->ksp,&its);
147:     tao->ksp_its += its;
148:     tao->ksp_tot_its+=its;
149:     if (ipmP->nb > 0) {
150:       /* Get max step size and apply frac-to-boundary */
151:       tau = PetscMax(ipmP->taumin,1.0-ipmP->mu);
152:       tau = PetscMin(tau,1.0);
153:       if (tau != 1.0) {
154:         VecScale(ipmP->s,tau);
155:         VecScale(ipmP->lamdai,tau);
156:       }
157:       VecStepBoundInfo(ipmP->s,ipmP->ds,ipmP->Zero_nb,ipmP->Inf_nb,&step_s,NULL,NULL);
158:       VecStepBoundInfo(ipmP->lamdai,ipmP->dlamdai,ipmP->Zero_nb,ipmP->Inf_nb,&step_l,NULL,NULL);
159:       if (tau != 1.0) {
160:         VecCopy(ipmP->save_s,ipmP->s);
161:         VecCopy(ipmP->save_lamdai,ipmP->lamdai);
162:       }
163:       alpha = PetscMin(step_s,step_l);
164:       alpha = PetscMin(alpha,1.0);
165:     } else {
166:       alpha = 1.0;
167:     }
168:     ipmP->alpha2 = alpha;
169:     /* TODO make phi_target meaningful */
170:     phi_target = ipmP->dec * ipmP->phi;
171:     for (i=0; i<11;i++) {
172:       VecAXPY(tao->solution,alpha,tao->stepdirection);
173:       if (ipmP->nb > 0) {
174:         VecAXPY(ipmP->s,alpha,ipmP->ds);
175:         VecAXPY(ipmP->lamdai,alpha,ipmP->dlamdai);
176:       }
177:       if (ipmP->me > 0) {
178:         VecAXPY(ipmP->lamdae,alpha,ipmP->dlamdae);
179:       }

181:       /* update dual variables */
182:       if (ipmP->me > 0) {
183:         VecCopy(ipmP->lamdae,tao->DE);
184:       }

186:       IPMEvaluate(tao);
187:       IPMComputeKKT(tao);
188:       if (ipmP->phi <= phi_target) break;
189:       alpha /= 2.0;
190:     }

192:     TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,stepsize,&reason);
193:     tao->niter++;
194:   }
195:   return(0);
196: }

200: static PetscErrorCode TaoSetup_IPM(Tao tao)
201: {
202:   TAO_IPM        *ipmP = (TAO_IPM*)tao->data;

206:   ipmP->nb = ipmP->mi = ipmP->me = 0;
207:   ipmP->K=0;
208:   VecGetSize(tao->solution,&ipmP->n);
209:   if (!tao->gradient) {
210:     VecDuplicate(tao->solution, &tao->gradient);
211:     VecDuplicate(tao->solution, &tao->stepdirection);
212:     VecDuplicate(tao->solution, &ipmP->rd);
213:     VecDuplicate(tao->solution, &ipmP->rhs_x);
214:     VecDuplicate(tao->solution, &ipmP->work);
215:     VecDuplicate(tao->solution, &ipmP->save_x);
216:   }
217:   if (tao->constraints_equality) {
218:     VecGetSize(tao->constraints_equality,&ipmP->me);
219:     VecDuplicate(tao->constraints_equality,&ipmP->lamdae);
220:     VecDuplicate(tao->constraints_equality,&ipmP->dlamdae);
221:     VecDuplicate(tao->constraints_equality,&ipmP->rhs_lamdae);
222:     VecDuplicate(tao->constraints_equality,&ipmP->save_lamdae);
223:     VecDuplicate(tao->constraints_equality,&ipmP->rpe);
224:     VecDuplicate(tao->constraints_equality,&tao->DE);
225:   }
226:   if (tao->constraints_inequality) {
227:     VecDuplicate(tao->constraints_inequality,&tao->DI);
228:   }
229:   return(0);
230: }

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,send;
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;
250:   MPI_Comm       comm;

253:   cind=xind=ucind=uceind=stepind=0;
254:   ipmP->mi=0;
255:   ipmP->nxlb=0;
256:   ipmP->nxub=0;
257:   ipmP->nb=0;
258:   ipmP->nslack=0;

260:   VecDuplicate(tao->solution,&xtmp);
261:   if (!tao->XL && !tao->XU && tao->ops->computebounds) {
262:     TaoComputeVariableBounds(tao);
263:   }
264:   if (tao->XL) {
265:     VecSet(xtmp,PETSC_NINFINITY);
266:     VecWhichGreaterThan(tao->XL,xtmp,&ipmP->isxl);
267:     ISGetSize(ipmP->isxl,&ipmP->nxlb);
268:   } else {
269:     ipmP->nxlb=0;
270:   }
271:   if (tao->XU) {
272:     VecSet(xtmp,PETSC_INFINITY);
273:     VecWhichLessThan(tao->XU,xtmp,&ipmP->isxu);
274:     ISGetSize(ipmP->isxu,&ipmP->nxub);
275:   } else {
276:     ipmP->nxub=0;
277:   }
278:   VecDestroy(&xtmp);
279:   if (tao->constraints_inequality) {
280:     VecGetSize(tao->constraints_inequality,&ipmP->mi);
281:   } else {
282:     ipmP->mi = 0;
283:   }
284:   ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;

286:   comm = ((PetscObject)(tao->solution))->comm;

288:   bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
289:   PetscMalloc1(bigsize,&stepind);
290:   PetscMalloc1(ipmP->n,&xind);
291:   PetscMalloc1(ipmP->me,&uceind);
292:   VecGetOwnershipRange(tao->solution,&xstart,&xend);

294:   if (ipmP->nb > 0) {
295:     VecCreate(comm,&ipmP->s);
296:     VecSetSizes(ipmP->s,PETSC_DECIDE,ipmP->nb);
297:     VecSetFromOptions(ipmP->s);
298:     VecDuplicate(ipmP->s,&ipmP->ds);
299:     VecDuplicate(ipmP->s,&ipmP->rhs_s);
300:     VecDuplicate(ipmP->s,&ipmP->complementarity);
301:     VecDuplicate(ipmP->s,&ipmP->ci);

303:     VecDuplicate(ipmP->s,&ipmP->lamdai);
304:     VecDuplicate(ipmP->s,&ipmP->dlamdai);
305:     VecDuplicate(ipmP->s,&ipmP->rhs_lamdai);
306:     VecDuplicate(ipmP->s,&ipmP->save_lamdai);

308:     VecDuplicate(ipmP->s,&ipmP->save_s);
309:     VecDuplicate(ipmP->s,&ipmP->rpi);
310:     VecDuplicate(ipmP->s,&ipmP->Zero_nb);
311:     VecSet(ipmP->Zero_nb,0.0);
312:     VecDuplicate(ipmP->s,&ipmP->One_nb);
313:     VecSet(ipmP->One_nb,1.0);
314:     VecDuplicate(ipmP->s,&ipmP->Inf_nb);
315:     VecSet(ipmP->Inf_nb,PETSC_INFINITY);

317:     PetscMalloc1(ipmP->nb,&cind);
318:     PetscMalloc1(ipmP->mi,&ucind);
319:     VecGetOwnershipRange(ipmP->s,&sstart,&send);

321:     if (ipmP->mi > 0) {
322:       VecGetOwnershipRange(tao->constraints_inequality,&ucstart,&ucend);
323:       counter=0;
324:       for (i=ucstart;i<ucend;i++) {
325:         cind[counter++] = i;
326:       }
327:       ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isuc);
328:       ISCreateGeneral(comm,counter,cind,PETSC_COPY_VALUES,&isc);
329:       VecScatterCreate(tao->constraints_inequality,isuc,ipmP->ci,isc,&ipmP->ci_scat);

331:       ISDestroy(&isuc);
332:       ISDestroy(&isc);
333:     }
334:     /* need to know how may xbound indices are on each process */
335:     /* TODO better way */
336:     if (ipmP->nxlb) {
337:       ISAllGather(ipmP->isxl,&bigxl);
338:       ISGetIndices(bigxl,&xli);
339:       /* find offsets for this processor */
340:       xl_offset = ipmP->mi;
341:       for (i=0;i<ipmP->nxlb;i++) {
342:         if (xli[i] < xstart) {
343:           xl_offset++;
344:         } else break;
345:       }
346:       ISRestoreIndices(bigxl,&xli);

348:       ISGetIndices(ipmP->isxl,&xli);
349:       ISGetLocalSize(ipmP->isxl,&nloc);
350:       for (i=0;i<nloc;i++) {
351:         xind[i] = xli[i];
352:         cind[i] = xl_offset+i;
353:       }

355:       ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
356:       ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
357:       VecScatterCreate(tao->XL,isx,ipmP->ci,isc,&ipmP->xl_scat);
358:       ISDestroy(&isx);
359:       ISDestroy(&isc);
360:       ISDestroy(&bigxl);
361:     }

363:     if (ipmP->nxub) {
364:       ISAllGather(ipmP->isxu,&bigxu);
365:       ISGetIndices(bigxu,&xui);
366:       /* find offsets for this processor */
367:       xu_offset = ipmP->mi + ipmP->nxlb;
368:       for (i=0;i<ipmP->nxub;i++) {
369:         if (xui[i] < xstart) {
370:           xu_offset++;
371:         } else break;
372:       }
373:       ISRestoreIndices(bigxu,&xui);

375:       ISGetIndices(ipmP->isxu,&xui);
376:       ISGetLocalSize(ipmP->isxu,&nloc);
377:       for (i=0;i<nloc;i++) {
378:         xind[i] = xui[i];
379:         cind[i] = xu_offset+i;
380:       }

382:       ISCreateGeneral(comm,nloc,xind,PETSC_COPY_VALUES,&isx);
383:       ISCreateGeneral(comm,nloc,cind,PETSC_COPY_VALUES,&isc);
384:       VecScatterCreate(tao->XU,isx,ipmP->ci,isc,&ipmP->xu_scat);
385:       ISDestroy(&isx);
386:       ISDestroy(&isc);
387:       ISDestroy(&bigxu);
388:     }
389:   }
390:   VecCreate(comm,&ipmP->bigrhs);
391:   VecGetType(tao->solution,&vtype);
392:   VecSetType(ipmP->bigrhs,vtype);
393:   VecSetSizes(ipmP->bigrhs,PETSC_DECIDE,bigsize);
394:   VecSetFromOptions(ipmP->bigrhs);
395:   VecDuplicate(ipmP->bigrhs,&ipmP->bigstep);

397:   /* create scatters for step->x and x->rhs */
398:   for (i=xstart;i<xend;i++) {
399:     stepind[i-xstart] = i;
400:     xind[i-xstart] = i;
401:   }
402:   ISCreateGeneral(comm,xend-xstart,stepind,PETSC_COPY_VALUES,&sis);
403:   ISCreateGeneral(comm,xend-xstart,xind,PETSC_COPY_VALUES,&is1);
404:   VecScatterCreate(ipmP->bigstep,sis,tao->solution,is1,&ipmP->step1);
405:   VecScatterCreate(tao->solution,is1,ipmP->bigrhs,sis,&ipmP->rhs1);
406:   ISDestroy(&sis);
407:   ISDestroy(&is1);

409:   if (ipmP->nb > 0) {
410:     for (i=sstart;i<send;i++) {
411:       stepind[i-sstart] = i+ipmP->n;
412:       cind[i-sstart] = i;
413:     }
414:     ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
415:     ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
416:     VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step2);
417:     ISDestroy(&sis);

419:     for (i=sstart;i<send;i++) {
420:       stepind[i-sstart] = i+ipmP->n+ipmP->me;
421:       cind[i-sstart] = i;
422:     }
423:     ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
424:     VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs3);
425:     ISDestroy(&sis);
426:     ISDestroy(&is1);
427:   }

429:   if (ipmP->me > 0) {
430:     VecGetOwnershipRange(tao->constraints_equality,&ucestart,&uceend);
431:     for (i=ucestart;i<uceend;i++) {
432:       stepind[i-ucestart] = i + ipmP->n+ipmP->nb;
433:       uceind[i-ucestart] = i;
434:     }

436:     ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
437:     ISCreateGeneral(comm,uceend-ucestart,uceind,PETSC_COPY_VALUES,&is1);
438:     VecScatterCreate(ipmP->bigstep,sis,tao->constraints_equality,is1,&ipmP->step3);
439:     ISDestroy(&sis);

441:     for (i=ucestart;i<uceend;i++) {
442:       stepind[i-ucestart] = i + ipmP->n;
443:     }

445:     ISCreateGeneral(comm,uceend-ucestart,stepind,PETSC_COPY_VALUES,&sis);
446:     VecScatterCreate(tao->constraints_equality,is1,ipmP->bigrhs,sis,&ipmP->rhs2);
447:     ISDestroy(&sis);
448:     ISDestroy(&is1);
449:   }

451:   if (ipmP->nb > 0) {
452:     for (i=sstart;i<send;i++) {
453:       stepind[i-sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
454:       cind[i-sstart] = i;
455:     }
456:     ISCreateGeneral(comm,send-sstart,cind,PETSC_COPY_VALUES,&is1);
457:     ISCreateGeneral(comm,send-sstart,stepind,PETSC_COPY_VALUES,&sis);
458:     VecScatterCreate(ipmP->bigstep,sis,ipmP->s,is1,&ipmP->step4);
459:     VecScatterCreate(ipmP->s,is1,ipmP->bigrhs,sis,&ipmP->rhs4);
460:     ISDestroy(&sis);
461:     ISDestroy(&is1);
462:   }

464:   PetscFree(stepind);
465:   PetscFree(cind);
466:   PetscFree(ucind);
467:   PetscFree(uceind);
468:   PetscFree(xind);
469:   return(0);
470: }

474: static PetscErrorCode TaoDestroy_IPM(Tao tao)
475: {
476:   TAO_IPM        *ipmP = (TAO_IPM*)tao->data;

480:   VecDestroy(&ipmP->rd);
481:   VecDestroy(&ipmP->rpe);
482:   VecDestroy(&ipmP->rpi);
483:   VecDestroy(&ipmP->work);
484:   VecDestroy(&ipmP->lamdae);
485:   VecDestroy(&ipmP->lamdai);
486:   VecDestroy(&ipmP->s);
487:   VecDestroy(&ipmP->ds);
488:   VecDestroy(&ipmP->ci);

490:   VecDestroy(&ipmP->rhs_x);
491:   VecDestroy(&ipmP->rhs_lamdae);
492:   VecDestroy(&ipmP->rhs_lamdai);
493:   VecDestroy(&ipmP->rhs_s);

495:   VecDestroy(&ipmP->save_x);
496:   VecDestroy(&ipmP->save_lamdae);
497:   VecDestroy(&ipmP->save_lamdai);
498:   VecDestroy(&ipmP->save_s);

500:   VecScatterDestroy(&ipmP->step1);
501:   VecScatterDestroy(&ipmP->step2);
502:   VecScatterDestroy(&ipmP->step3);
503:   VecScatterDestroy(&ipmP->step4);

505:   VecScatterDestroy(&ipmP->rhs1);
506:   VecScatterDestroy(&ipmP->rhs2);
507:   VecScatterDestroy(&ipmP->rhs3);
508:   VecScatterDestroy(&ipmP->rhs4);

510:   VecScatterDestroy(&ipmP->ci_scat);
511:   VecScatterDestroy(&ipmP->xl_scat);
512:   VecScatterDestroy(&ipmP->xu_scat);

514:   VecDestroy(&ipmP->dlamdai);
515:   VecDestroy(&ipmP->dlamdae);
516:   VecDestroy(&ipmP->Zero_nb);
517:   VecDestroy(&ipmP->One_nb);
518:   VecDestroy(&ipmP->Inf_nb);
519:   VecDestroy(&ipmP->complementarity);

521:   VecDestroy(&ipmP->bigrhs);
522:   VecDestroy(&ipmP->bigstep);
523:   MatDestroy(&ipmP->Ai);
524:   MatDestroy(&ipmP->K);
525:   ISDestroy(&ipmP->isxu);
526:   ISDestroy(&ipmP->isxl);
527:   PetscFree(tao->data);
528:   return(0);
529: }

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

550: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
551: {
552:   return 0;
553: }

555: /* IPMObjectiveAndGradient()
556:    f = d'x + 0.5 * x' * H * x
557:    rd = H*x + d + Ae'*lame - Ai'*lami
558:    rpe = Ae*x - be
559:    rpi = Ai*x - yi - bi
560:    mu = yi' * lami/mi;
561:    com = yi.*lami

563:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
564: */
565: /*
568: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
569: {
570:   Tao tao = (Tao)tptr;
571:   TAO_IPM *ipmP = (TAO_IPM*)tao->data;
574:   IPMComputeKKT(tao);
575:   *f = ipmP->phi;
576:   return(0);
577: }
578: */

580: /*
581:    f = d'x + 0.5 * x' * H * x
582:    rd = H*x + d + Ae'*lame - Ai'*lami
583:        Ai =   jac_ineq
584:                I (w/lb)
585:               -I (w/ub)

587:    rpe = ce
588:    rpi = ci - s;
589:    com = s.*lami
590:    mu = yi' * lami/mi;

592:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
593: */
596: static PetscErrorCode IPMComputeKKT(Tao tao)
597: {
598:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;
599:   PetscScalar    norm;

603:   VecCopy(tao->gradient,ipmP->rd);

605:   if (ipmP->me > 0) {
606:     /* rd = gradient + Ae'*lamdae */
607:     MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
608:     VecAXPY(ipmP->rd, 1.0, ipmP->work);

610:     /* rpe = ce(x) */
611:     VecCopy(tao->constraints_equality,ipmP->rpe);
612:   }
613:   if (ipmP->nb > 0) {
614:     /* rd = rd - Ai'*lamdai */
615:     MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
616:     VecAXPY(ipmP->rd, -1.0, ipmP->work);

618:     /* rpi = cin - s */
619:     VecCopy(ipmP->ci,ipmP->rpi);
620:     VecAXPY(ipmP->rpi, -1.0, ipmP->s);

622:     /* com = s .* lami */
623:     VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
624:   }
625:   /* phi = ||rd; rpe; rpi; com|| */
626:   VecDot(ipmP->rd,ipmP->rd,&norm);
627:   ipmP->phi = norm;
628:   if (ipmP->me > 0 ) {
629:     VecDot(ipmP->rpe,ipmP->rpe,&norm);
630:     ipmP->phi += norm;
631:   }
632:   if (ipmP->nb > 0) {
633:     VecDot(ipmP->rpi,ipmP->rpi,&norm);
634:     ipmP->phi += norm;
635:     VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
636:     ipmP->phi += norm;
637:     /* mu = s'*lami/nb */
638:     VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
639:     ipmP->mu /= ipmP->nb;
640:   } else {
641:     ipmP->mu = 1.0;
642:   }

644:   ipmP->phi = PetscSqrtScalar(ipmP->phi);
645:   return(0);
646: }

650: /* evaluate user info at current point */
651: PetscErrorCode IPMEvaluate(Tao tao)
652: {
653:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

657:   TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
658:   TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
659:   if (ipmP->me > 0) {
660:     TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
661:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
662:   }
663:   if (ipmP->mi > 0) {
664:     TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
665:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
666:   }
667:   if (ipmP->nb > 0) {
668:     /* Ai' =   jac_ineq | I (w/lb) | -I (w/ub)  */
669:     IPMUpdateAi(tao);
670:   }
671:   return(0);
672: }

676: /* Push initial point away from bounds */
677: PetscErrorCode IPMPushInitialPoint(Tao tao)
678: {
679:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

683:   TaoComputeVariableBounds(tao);
684:   if (tao->XL && tao->XU) {
685:     VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
686:   }
687:   if (ipmP->nb > 0) {
688:     VecSet(ipmP->s,ipmP->pushs);
689:     VecSet(ipmP->lamdai,ipmP->pushnu);
690:     if (ipmP->mi > 0) {
691:       VecSet(tao->DI,ipmP->pushnu);
692:     }
693:   }
694:   if (ipmP->me > 0) {
695:     VecSet(tao->DE,1.0);
696:     VecSet(ipmP->lamdae,1.0);
697:   }
698:   return(0);
699: }

703: PetscErrorCode IPMUpdateAi(Tao tao)
704: {
705:   /* Ai =     Ji
706:               I (w/lb)
707:              -I (w/ub) */

709:   /* Ci =    user->ci
710:              Xi - lb (w/lb)
711:              -Xi + ub (w/ub)  */

713:   TAO_IPM           *ipmP = (TAO_IPM *)tao->data;
714:   MPI_Comm          comm;
715:   PetscInt          i;
716:   PetscScalar       newval;
717:   PetscInt          newrow,newcol,ncols;
718:   const PetscScalar *vals;
719:   const PetscInt    *cols;
720:   PetscInt          astart,aend,jstart,jend;
721:   PetscInt          *nonzeros;
722:   PetscInt          r2,r3,r4;
723:   PetscMPIInt       size;
724:   PetscErrorCode    ierr;

727:   r2 = ipmP->mi;
728:   r3 = r2 + ipmP->nxlb;
729:   r4 = r3 + ipmP->nxub;

731:   if (!ipmP->nb) return(0);

733:   /* Create Ai matrix if it doesn't exist yet */
734:   if (!ipmP->Ai) {
735:     comm = ((PetscObject)(tao->solution))->comm;
736:     PetscMalloc1(ipmP->nb,&nonzeros);
737:     MPI_Comm_size(comm,&size);
738:     if (size == 1) {
739:       for (i=0;i<ipmP->mi;i++) {
740:         MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
741:         nonzeros[i] = ncols;
742:         MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
743:       }
744:       for (i=r2;i<r4;i++) {
745:         nonzeros[i] = 1;
746:       }
747:     }
748:     MatCreate(comm,&ipmP->Ai);
749:     MatSetType(ipmP->Ai,MATAIJ);
750:     MatSetSizes(ipmP->Ai,PETSC_DECIDE,PETSC_DECIDE,ipmP->nb,ipmP->n);
751:     MatSetFromOptions(ipmP->Ai);
752:     MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
753:     MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
754:     if (size ==1) {
755:       PetscFree(nonzeros);
756:     }
757:   }

759:   /* Copy values from user jacobian to Ai */
760:   MatGetOwnershipRange(ipmP->Ai,&astart,&aend);

762:   /* Ai w/lb */
763:   if (ipmP->mi) {
764:     MatZeroEntries(ipmP->Ai);
765:     MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
766:     for (i=jstart;i<jend;i++) {
767:       MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
768:       newrow = i;
769:       MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
770:       MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
771:     }
772:   }

774:   /* I w/ xlb */
775:   if (ipmP->nxlb) {
776:     for (i=0;i<ipmP->nxlb;i++) {
777:       if (i>=astart && i<aend) {
778:         newrow = i+r2;
779:         newcol = i;
780:         newval = 1.0;
781:         MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
782:       }
783:     }
784:   }
785:   if (ipmP->nxub) {
786:     /* I w/ xub */
787:     for (i=0;i<ipmP->nxub;i++) {
788:       if (i>=astart && i<aend) {
789:       newrow = i+r3;
790:       newcol = i;
791:       newval = -1.0;
792:       MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
793:       }
794:     }
795:   }

797:   MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
798:   MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
799:   CHKMEMQ;

801:   VecSet(ipmP->ci,0.0);

803:   /* user ci */
804:   if (ipmP->mi > 0) {
805:     VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
806:     VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
807:   }
808:   if (!ipmP->work){
809:     VecDuplicate(tao->solution,&ipmP->work);
810:   }
811:   VecCopy(tao->solution,ipmP->work);
812:   if (tao->XL) {
813:     VecAXPY(ipmP->work,-1.0,tao->XL);

815:     /* lower bounds on variables */
816:     if (ipmP->nxlb > 0) {
817:       VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
818:       VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
819:     }
820:   }
821:   if (tao->XU) {
822:     /* upper bounds on variables */
823:     VecCopy(tao->solution,ipmP->work);
824:     VecScale(ipmP->work,-1.0);
825:     VecAXPY(ipmP->work,1.0,tao->XU);
826:     if (ipmP->nxub > 0) {
827:       VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
828:       VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
829:     }
830:   }
831:   return(0);
832: }

836: /* create K = [ Hlag , 0 , Ae', -Ai'];
837:               [Ae , 0,   0  , 0];
838:               [Ai ,-I,   0 ,  0];
839:               [ 0 , S ,  0,   Y ];  */
840: PetscErrorCode IPMUpdateK(Tao tao)
841: {
842:   TAO_IPM         *ipmP = (TAO_IPM *)tao->data;
843:   MPI_Comm        comm;
844:   PetscMPIInt     size;
845:   PetscErrorCode  ierr;
846:   PetscInt        i,j,row;
847:   PetscInt        ncols,newcol,newcols[2],newrow;
848:   const PetscInt  *cols;
849:   const PetscReal *vals;
850:   const PetscReal *l,*y;
851:   PetscReal       *newvals;
852:   PetscReal       newval;
853:   PetscInt        subsize;
854:   const PetscInt  *indices;
855:   PetscInt        *nonzeros,*d_nonzeros,*o_nonzeros;
856:   PetscInt        bigsize;
857:   PetscInt        r1,r2,r3;
858:   PetscInt        c1,c2,c3;
859:   PetscInt        klocalsize;
860:   PetscInt        hstart,hend,kstart,kend;
861:   PetscInt        aistart,aiend,aestart,aeend;
862:   PetscInt        sstart,send;

865:   comm = ((PetscObject)(tao->solution))->comm;
866:   MPI_Comm_size(comm,&size);
867:   IPMUpdateAi(tao);

869:   /* allocate workspace */
870:   subsize = PetscMax(ipmP->n,ipmP->nb);
871:   subsize = PetscMax(ipmP->me,subsize);
872:   subsize = PetscMax(2,subsize);
873:   PetscMalloc1(subsize,&indices);
874:   PetscMalloc1(subsize,&newvals);

876:   r1 = c1 = ipmP->n;
877:   r2 = r1 + ipmP->me;  c2 = c1 + ipmP->nb;
878:   r3 = c3 = r2 + ipmP->nb;

880:   bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
881:   VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
882:   MatGetOwnershipRange(tao->hessian,&hstart,&hend);
883:   klocalsize = kend-kstart;
884:   if (!ipmP->K) {
885:     if (size == 1) {
886:       PetscMalloc1(kend-kstart,&nonzeros);
887:       for (i=0;i<bigsize;i++) {
888:         if (i<r1) {
889:           MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
890:           nonzeros[i] = ncols;
891:           MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
892:           nonzeros[i] += ipmP->me+ipmP->nb;
893:         } else if (i<r2) {
894:           nonzeros[i-kstart] = ipmP->n;
895:         } else if (i<r3) {
896:           nonzeros[i-kstart] = ipmP->n+1;
897:         } else if (i<bigsize) {
898:           nonzeros[i-kstart] = 2;
899:         }
900:       }
901:       MatCreate(comm,&ipmP->K);
902:       MatSetType(ipmP->K,MATSEQAIJ);
903:       MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
904:       MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
905:       MatSetFromOptions(ipmP->K);
906:       PetscFree(nonzeros);
907:     } else {
908:       PetscMalloc1(kend-kstart,&d_nonzeros);
909:       PetscMalloc1(kend-kstart,&o_nonzeros);
910:       for (i=kstart;i<kend;i++) {
911:         if (i<r1) {
912:           /* TODO fix preallocation for mpi mats */
913:           d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
914:           o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
915:         } else if (i<r2) {
916:           d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
917:           o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
918:         } else if (i<r3) {
919:           d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
920:           o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
921:         } else {
922:           d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
923:           o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
924:         }
925:       }
926:       MatCreate(comm,&ipmP->K);
927:       MatSetType(ipmP->K,MATMPIAIJ);
928:       MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
929:       MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
930:       PetscFree(d_nonzeros);
931:       PetscFree(o_nonzeros);
932:       MatSetFromOptions(ipmP->K);
933:     }
934:   }

936:   MatZeroEntries(ipmP->K);
937:   /* Copy H */
938:   for (i=hstart;i<hend;i++) {
939:     MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
940:     if (ncols > 0) {
941:       MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
942:     }
943:     MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
944:   }

946:   /* Copy Ae and Ae' */
947:   if (ipmP->me > 0) {
948:     MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
949:     for (i=aestart;i<aeend;i++) {
950:       MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
951:       if (ncols > 0) {
952:         /*Ae*/
953:         row = i+r1;
954:         MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
955:         /*Ae'*/
956:         for (j=0;j<ncols;j++) {
957:           newcol = i + c2;
958:           newrow = cols[j];
959:           newval = vals[j];
960:           MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
961:         }
962:       }
963:       MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
964:     }
965:   }

967:   if (ipmP->nb > 0) {
968:     MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
969:     /* Copy Ai,and Ai' */
970:     for (i=aistart;i<aiend;i++) {
971:       row = i+r2;
972:       MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
973:       if (ncols > 0) {
974:         /*Ai*/
975:         MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
976:         /*-Ai'*/
977:         for (j=0;j<ncols;j++) {
978:           newcol = i + c3;
979:           newrow = cols[j];
980:           newval = -vals[j];
981:           MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
982:         }
983:       }
984:       MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
985:     }

987:     /* -I */
988:     for (i=kstart;i<kend;i++) {
989:       if (i>=r2 && i<r3) {
990:         newrow = i;
991:         newcol = i-r2+c1;
992:         newval = -1.0;
993:         MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
994:       }
995:     }

997:     /* Copy L,Y */
998:     VecGetOwnershipRange(ipmP->s,&sstart,&send);
999:     VecGetArrayRead(ipmP->lamdai,&l);
1000:     VecGetArrayRead(ipmP->s,&y);

1002:     for (i=sstart;i<send;i++) {
1003:       newcols[0] = c1+i;
1004:       newcols[1] = c3+i;
1005:       newvals[0] = l[i-sstart];
1006:       newvals[1] = y[i-sstart];
1007:       newrow = r3+i;
1008:       MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
1009:     }

1011:     VecRestoreArrayRead(ipmP->lamdai,&l);
1012:     VecRestoreArrayRead(ipmP->s,&y);
1013:   }

1015:   PetscFree(indices);
1016:   PetscFree(newvals);
1017:   MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
1018:   MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
1019:   return(0);
1020: }

1024: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
1025: {
1026:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

1030:   /* rhs = [x1      (n)
1031:             x2     (me)
1032:             x3     (nb)
1033:             x4     (nb)] */
1034:   if (X1) {
1035:     VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1036:     VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1037:   }
1038:   if (ipmP->me > 0 && X2) {
1039:     VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1040:     VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1041:   }
1042:   if (ipmP->nb > 0) {
1043:     if (X3) {
1044:       VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1045:       VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1046:     }
1047:     if (X4) {
1048:       VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1049:       VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1050:     }
1051:   }
1052:   return(0);
1053: }

1057: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1058: {
1059:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

1063:   CHKMEMQ;
1064:   /*        [x1    (n)
1065:              x2    (nb) may be 0
1066:              x3    (me) may be 0
1067:              x4    (nb) may be 0 */
1068:   if (X1) {
1069:     VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1070:     VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1071:   }
1072:   if (X2 && ipmP->nb > 0) {
1073:     VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1074:     VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1075:   }
1076:   if (X3 && ipmP->me > 0) {
1077:     VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1078:     VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1079:   }
1080:   if (X4 && ipmP->nb > 0) {
1081:     VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1082:     VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1083:   }
1084:   CHKMEMQ;
1085:   return(0);
1086: }

1088: /*MC
1089:   TAOIPM - Interior point algorithm for generally constrained optimization.

1091:   Option Database Keys:
1092: +   -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1093: .   -tao_ipm_pushs - parameter to push initial slack variables away from bounds

1095:   Notes: 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.
1096:   Level: beginner

1098: M*/

1102: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1103: {
1104:   TAO_IPM        *ipmP;

1108:   tao->ops->setup = TaoSetup_IPM;
1109:   tao->ops->solve = TaoSolve_IPM;
1110:   tao->ops->view = TaoView_IPM;
1111:   tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1112:   tao->ops->destroy = TaoDestroy_IPM;
1113:   /* tao->ops->computedual = TaoComputeDual_IPM; */

1115:   PetscNewLog(tao,&ipmP);
1116:   tao->data = (void*)ipmP;

1118:   /* Override default settings (unless already changed) */
1119:   if (!tao->max_it_changed) tao->max_it = 200;
1120:   if (!tao->max_funcs_changed) tao->max_funcs = 500;

1122:   ipmP->dec = 10000; /* line search critera */
1123:   ipmP->taumin = 0.995;
1124:   ipmP->monitorkkt = PETSC_FALSE;
1125:   ipmP->pushs = 100;
1126:   ipmP->pushnu = 100;
1127:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1128:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1129:   return(0);
1130: }