Actual source code: ipm.c

petsc-3.9.4 2018-09-11
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:     tao->ksp_its=0;
 56:     IPMUpdateK(tao);
 57:     /*
 58:        rhs.x    = -rd
 59:        rhs.lame = -rpe
 60:        rhs.lami = -rpi
 61:        rhs.com  = -com
 62:     */

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

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

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

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

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

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

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

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

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

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

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

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

193:     TaoLogConvergenceHistory(tao,ipmP->kkt_f,ipmP->phi,0.0,tao->ksp_its);
194:     TaoMonitor(tao,tao->niter,ipmP->kkt_f,ipmP->phi,0.0,stepsize);
195:     (*tao->ops->convergencetest)(tao,tao->cnvP);
196:     tao->niter++;
197:   }
198:   return(0);
199: }

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

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

233: static PetscErrorCode IPMInitializeBounds(Tao tao)
234: {
235:   TAO_IPM        *ipmP = (TAO_IPM*)tao->data;
236:   Vec            xtmp;
237:   PetscInt       xstart,xend;
238:   PetscInt       ucstart,ucend; /* user ci */
239:   PetscInt       ucestart,uceend; /* user ce */
240:   PetscInt       sstart,send;
241:   PetscInt       bigsize;
242:   PetscInt       i,counter,nloc;
243:   PetscInt       *cind,*xind,*ucind,*uceind,*stepind;
244:   VecType        vtype;
245:   const PetscInt *xli,*xui;
246:   PetscInt       xl_offset,xu_offset;
247:   IS             bigxl,bigxu,isuc,isc,isx,sis,is1;
249:   MPI_Comm       comm;

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

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

285:   PetscObjectGetComm((PetscObject)tao->solution,&comm);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

471: static PetscErrorCode TaoDestroy_IPM(Tao tao)
472: {
473:   TAO_IPM        *ipmP = (TAO_IPM*)tao->data;

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

487:   VecDestroy(&ipmP->rhs_x);
488:   VecDestroy(&ipmP->rhs_lamdae);
489:   VecDestroy(&ipmP->rhs_lamdai);
490:   VecDestroy(&ipmP->rhs_s);

492:   VecDestroy(&ipmP->save_x);
493:   VecDestroy(&ipmP->save_lamdae);
494:   VecDestroy(&ipmP->save_lamdai);
495:   VecDestroy(&ipmP->save_s);

497:   VecScatterDestroy(&ipmP->step1);
498:   VecScatterDestroy(&ipmP->step2);
499:   VecScatterDestroy(&ipmP->step3);
500:   VecScatterDestroy(&ipmP->step4);

502:   VecScatterDestroy(&ipmP->rhs1);
503:   VecScatterDestroy(&ipmP->rhs2);
504:   VecScatterDestroy(&ipmP->rhs3);
505:   VecScatterDestroy(&ipmP->rhs4);

507:   VecScatterDestroy(&ipmP->ci_scat);
508:   VecScatterDestroy(&ipmP->xl_scat);
509:   VecScatterDestroy(&ipmP->xu_scat);

511:   VecDestroy(&ipmP->dlamdai);
512:   VecDestroy(&ipmP->dlamdae);
513:   VecDestroy(&ipmP->Zero_nb);
514:   VecDestroy(&ipmP->One_nb);
515:   VecDestroy(&ipmP->Inf_nb);
516:   VecDestroy(&ipmP->complementarity);

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

528: static PetscErrorCode TaoSetFromOptions_IPM(PetscOptionItems *PetscOptionsObject,Tao tao)
529: {
530:   TAO_IPM        *ipmP = (TAO_IPM*)tao->data;

534:   PetscOptionsHead(PetscOptionsObject,"IPM method for constrained optimization");
535:   PetscOptionsBool("-tao_ipm_monitorkkt","monitor kkt status",NULL,ipmP->monitorkkt,&ipmP->monitorkkt,NULL);
536:   PetscOptionsReal("-tao_ipm_pushs","parameter to push initial slack variables away from bounds",NULL,ipmP->pushs,&ipmP->pushs,NULL);
537:   PetscOptionsReal("-tao_ipm_pushnu","parameter to push initial (inequality) dual variables away from bounds",NULL,ipmP->pushnu,&ipmP->pushnu,NULL);
538:   PetscOptionsTail();
539:   KSPSetFromOptions(tao->ksp);
540:   return(0);
541: }

543: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
544: {
545:   return 0;
546: }

548: /* IPMObjectiveAndGradient()
549:    f = d'x + 0.5 * x' * H * x
550:    rd = H*x + d + Ae'*lame - Ai'*lami
551:    rpe = Ae*x - be
552:    rpi = Ai*x - yi - bi
553:    mu = yi' * lami/mi;
554:    com = yi.*lami

556:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
557: */
558: /*
559: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
560: {
561:   Tao tao = (Tao)tptr;
562:   TAO_IPM *ipmP = (TAO_IPM*)tao->data;
565:   IPMComputeKKT(tao);
566:   *f = ipmP->phi;
567:   return(0);
568: }
569: */

571: /*
572:    f = d'x + 0.5 * x' * H * x
573:    rd = H*x + d + Ae'*lame - Ai'*lami
574:        Ai =   jac_ineq
575:                I (w/lb)
576:               -I (w/ub)

578:    rpe = ce
579:    rpi = ci - s;
580:    com = s.*lami
581:    mu = yi' * lami/mi;

583:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
584: */
585: static PetscErrorCode IPMComputeKKT(Tao tao)
586: {
587:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;
588:   PetscScalar    norm;

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

594:   if (ipmP->me > 0) {
595:     /* rd = gradient + Ae'*lamdae */
596:     MatMultTranspose(tao->jacobian_equality,ipmP->lamdae,ipmP->work);
597:     VecAXPY(ipmP->rd, 1.0, ipmP->work);

599:     /* rpe = ce(x) */
600:     VecCopy(tao->constraints_equality,ipmP->rpe);
601:   }
602:   if (ipmP->nb > 0) {
603:     /* rd = rd - Ai'*lamdai */
604:     MatMultTranspose(ipmP->Ai,ipmP->lamdai,ipmP->work);
605:     VecAXPY(ipmP->rd, -1.0, ipmP->work);

607:     /* rpi = cin - s */
608:     VecCopy(ipmP->ci,ipmP->rpi);
609:     VecAXPY(ipmP->rpi, -1.0, ipmP->s);

611:     /* com = s .* lami */
612:     VecPointwiseMult(ipmP->complementarity, ipmP->s,ipmP->lamdai);
613:   }
614:   /* phi = ||rd; rpe; rpi; com|| */
615:   VecDot(ipmP->rd,ipmP->rd,&norm);
616:   ipmP->phi = norm;
617:   if (ipmP->me > 0 ) {
618:     VecDot(ipmP->rpe,ipmP->rpe,&norm);
619:     ipmP->phi += norm;
620:   }
621:   if (ipmP->nb > 0) {
622:     VecDot(ipmP->rpi,ipmP->rpi,&norm);
623:     ipmP->phi += norm;
624:     VecDot(ipmP->complementarity,ipmP->complementarity,&norm);
625:     ipmP->phi += norm;
626:     /* mu = s'*lami/nb */
627:     VecDot(ipmP->s,ipmP->lamdai,&ipmP->mu);
628:     ipmP->mu /= ipmP->nb;
629:   } else {
630:     ipmP->mu = 1.0;
631:   }

633:   ipmP->phi = PetscSqrtScalar(ipmP->phi);
634:   return(0);
635: }

637: /* evaluate user info at current point */
638: PetscErrorCode IPMEvaluate(Tao tao)
639: {
640:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

644:   TaoComputeObjectiveAndGradient(tao,tao->solution,&ipmP->kkt_f,tao->gradient);
645:   TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
646:   if (ipmP->me > 0) {
647:     TaoComputeEqualityConstraints(tao,tao->solution,tao->constraints_equality);
648:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
649:   }
650:   if (ipmP->mi > 0) {
651:     TaoComputeInequalityConstraints(tao,tao->solution,tao->constraints_inequality);
652:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
653:   }
654:   if (ipmP->nb > 0) {
655:     /* Ai' =   jac_ineq | I (w/lb) | -I (w/ub)  */
656:     IPMUpdateAi(tao);
657:   }
658:   return(0);
659: }

661: /* Push initial point away from bounds */
662: PetscErrorCode IPMPushInitialPoint(Tao tao)
663: {
664:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

668:   TaoComputeVariableBounds(tao);
669:   if (tao->XL && tao->XU) {
670:     VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
671:   }
672:   if (ipmP->nb > 0) {
673:     VecSet(ipmP->s,ipmP->pushs);
674:     VecSet(ipmP->lamdai,ipmP->pushnu);
675:     if (ipmP->mi > 0) {
676:       VecSet(tao->DI,ipmP->pushnu);
677:     }
678:   }
679:   if (ipmP->me > 0) {
680:     VecSet(tao->DE,1.0);
681:     VecSet(ipmP->lamdae,1.0);
682:   }
683:   return(0);
684: }

686: PetscErrorCode IPMUpdateAi(Tao tao)
687: {
688:   /* Ai =     Ji
689:               I (w/lb)
690:              -I (w/ub) */

692:   /* Ci =    user->ci
693:              Xi - lb (w/lb)
694:              -Xi + ub (w/ub)  */

696:   TAO_IPM           *ipmP = (TAO_IPM *)tao->data;
697:   MPI_Comm          comm;
698:   PetscInt          i;
699:   PetscScalar       newval;
700:   PetscInt          newrow,newcol,ncols;
701:   const PetscScalar *vals;
702:   const PetscInt    *cols;
703:   PetscInt          astart,aend,jstart,jend;
704:   PetscInt          *nonzeros;
705:   PetscInt          r2,r3,r4;
706:   PetscMPIInt       size;
707:   PetscErrorCode    ierr;

710:   r2 = ipmP->mi;
711:   r3 = r2 + ipmP->nxlb;
712:   r4 = r3 + ipmP->nxub;

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

716:   /* Create Ai matrix if it doesn't exist yet */
717:   if (!ipmP->Ai) {
718:     comm = ((PetscObject)(tao->solution))->comm;
719:     PetscMalloc1(ipmP->nb,&nonzeros);
720:     MPI_Comm_size(comm,&size);
721:     if (size == 1) {
722:       for (i=0;i<ipmP->mi;i++) {
723:         MatGetRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
724:         nonzeros[i] = ncols;
725:         MatRestoreRow(tao->jacobian_inequality,i,&ncols,NULL,NULL);
726:       }
727:       for (i=r2;i<r4;i++) {
728:         nonzeros[i] = 1;
729:       }
730:     }
731:     MatCreate(comm,&ipmP->Ai);
732:     MatSetType(ipmP->Ai,MATAIJ);
733:     MatSetSizes(ipmP->Ai,PETSC_DECIDE,PETSC_DECIDE,ipmP->nb,ipmP->n);
734:     MatSetFromOptions(ipmP->Ai);
735:     MatMPIAIJSetPreallocation(ipmP->Ai,ipmP->nb,NULL,ipmP->nb,NULL);
736:     MatSeqAIJSetPreallocation(ipmP->Ai,PETSC_DEFAULT,nonzeros);
737:     if (size ==1) {
738:       PetscFree(nonzeros);
739:     }
740:   }

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

745:   /* Ai w/lb */
746:   if (ipmP->mi) {
747:     MatZeroEntries(ipmP->Ai);
748:     MatGetOwnershipRange(tao->jacobian_inequality,&jstart,&jend);
749:     for (i=jstart;i<jend;i++) {
750:       MatGetRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
751:       newrow = i;
752:       MatSetValues(ipmP->Ai,1,&newrow,ncols,cols,vals,INSERT_VALUES);
753:       MatRestoreRow(tao->jacobian_inequality,i,&ncols,&cols,&vals);
754:     }
755:   }

757:   /* I w/ xlb */
758:   if (ipmP->nxlb) {
759:     for (i=0;i<ipmP->nxlb;i++) {
760:       if (i>=astart && i<aend) {
761:         newrow = i+r2;
762:         newcol = i;
763:         newval = 1.0;
764:         MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
765:       }
766:     }
767:   }
768:   if (ipmP->nxub) {
769:     /* I w/ xub */
770:     for (i=0;i<ipmP->nxub;i++) {
771:       if (i>=astart && i<aend) {
772:       newrow = i+r3;
773:       newcol = i;
774:       newval = -1.0;
775:       MatSetValues(ipmP->Ai,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
776:       }
777:     }
778:   }

780:   MatAssemblyBegin(ipmP->Ai,MAT_FINAL_ASSEMBLY);
781:   MatAssemblyEnd(ipmP->Ai,MAT_FINAL_ASSEMBLY);
782:   CHKMEMQ;

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

786:   /* user ci */
787:   if (ipmP->mi > 0) {
788:     VecScatterBegin(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
789:     VecScatterEnd(ipmP->ci_scat,tao->constraints_inequality,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
790:   }
791:   if (!ipmP->work){
792:     VecDuplicate(tao->solution,&ipmP->work);
793:   }
794:   VecCopy(tao->solution,ipmP->work);
795:   if (tao->XL) {
796:     VecAXPY(ipmP->work,-1.0,tao->XL);

798:     /* lower bounds on variables */
799:     if (ipmP->nxlb > 0) {
800:       VecScatterBegin(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
801:       VecScatterEnd(ipmP->xl_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
802:     }
803:   }
804:   if (tao->XU) {
805:     /* upper bounds on variables */
806:     VecCopy(tao->solution,ipmP->work);
807:     VecScale(ipmP->work,-1.0);
808:     VecAXPY(ipmP->work,1.0,tao->XU);
809:     if (ipmP->nxub > 0) {
810:       VecScatterBegin(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
811:       VecScatterEnd(ipmP->xu_scat,ipmP->work,ipmP->ci,INSERT_VALUES,SCATTER_FORWARD);
812:     }
813:   }
814:   return(0);
815: }

817: /* create K = [ Hlag , 0 , Ae', -Ai'];
818:               [Ae , 0,   0  , 0];
819:               [Ai ,-I,   0 ,  0];
820:               [ 0 , S ,  0,   Y ];  */
821: PetscErrorCode IPMUpdateK(Tao tao)
822: {
823:   TAO_IPM         *ipmP = (TAO_IPM *)tao->data;
824:   MPI_Comm        comm;
825:   PetscMPIInt     size;
826:   PetscErrorCode  ierr;
827:   PetscInt        i,j,row;
828:   PetscInt        ncols,newcol,newcols[2],newrow;
829:   const PetscInt  *cols;
830:   const PetscReal *vals;
831:   const PetscReal *l,*y;
832:   PetscReal       *newvals;
833:   PetscReal       newval;
834:   PetscInt        subsize;
835:   const PetscInt  *indices;
836:   PetscInt        *nonzeros,*d_nonzeros,*o_nonzeros;
837:   PetscInt        bigsize;
838:   PetscInt        r1,r2,r3;
839:   PetscInt        c1,c2,c3;
840:   PetscInt        klocalsize;
841:   PetscInt        hstart,hend,kstart,kend;
842:   PetscInt        aistart,aiend,aestart,aeend;
843:   PetscInt        sstart,send;

846:   comm = ((PetscObject)(tao->solution))->comm;
847:   MPI_Comm_size(comm,&size);
848:   IPMUpdateAi(tao);

850:   /* allocate workspace */
851:   subsize = PetscMax(ipmP->n,ipmP->nb);
852:   subsize = PetscMax(ipmP->me,subsize);
853:   subsize = PetscMax(2,subsize);
854:   PetscMalloc1(subsize,(PetscInt**)&indices);
855:   PetscMalloc1(subsize,&newvals);

857:   r1 = c1 = ipmP->n;
858:   r2 = r1 + ipmP->me;  c2 = c1 + ipmP->nb;
859:   r3 = c3 = r2 + ipmP->nb;

861:   bigsize = ipmP->n+2*ipmP->nb+ipmP->me;
862:   VecGetOwnershipRange(ipmP->bigrhs,&kstart,&kend);
863:   MatGetOwnershipRange(tao->hessian,&hstart,&hend);
864:   klocalsize = kend-kstart;
865:   if (!ipmP->K) {
866:     if (size == 1) {
867:       PetscMalloc1(kend-kstart,&nonzeros);
868:       for (i=0;i<bigsize;i++) {
869:         if (i<r1) {
870:           MatGetRow(tao->hessian,i,&ncols,NULL,NULL);
871:           nonzeros[i] = ncols;
872:           MatRestoreRow(tao->hessian,i,&ncols,NULL,NULL);
873:           nonzeros[i] += ipmP->me+ipmP->nb;
874:         } else if (i<r2) {
875:           nonzeros[i-kstart] = ipmP->n;
876:         } else if (i<r3) {
877:           nonzeros[i-kstart] = ipmP->n+1;
878:         } else if (i<bigsize) {
879:           nonzeros[i-kstart] = 2;
880:         }
881:       }
882:       MatCreate(comm,&ipmP->K);
883:       MatSetType(ipmP->K,MATSEQAIJ);
884:       MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
885:       MatSeqAIJSetPreallocation(ipmP->K,0,nonzeros);
886:       MatSetFromOptions(ipmP->K);
887:       PetscFree(nonzeros);
888:     } else {
889:       PetscMalloc1(kend-kstart,&d_nonzeros);
890:       PetscMalloc1(kend-kstart,&o_nonzeros);
891:       for (i=kstart;i<kend;i++) {
892:         if (i<r1) {
893:           /* TODO fix preallocation for mpi mats */
894:           d_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,kend-kstart);
895:           o_nonzeros[i-kstart] = PetscMin(ipmP->n+ipmP->me+ipmP->nb,bigsize-(kend-kstart));
896:         } else if (i<r2) {
897:           d_nonzeros[i-kstart] = PetscMin(ipmP->n,kend-kstart);
898:           o_nonzeros[i-kstart] = PetscMin(ipmP->n,bigsize-(kend-kstart));
899:         } else if (i<r3) {
900:           d_nonzeros[i-kstart] = PetscMin(ipmP->n+2,kend-kstart);
901:           o_nonzeros[i-kstart] = PetscMin(ipmP->n+2,bigsize-(kend-kstart));
902:         } else {
903:           d_nonzeros[i-kstart] = PetscMin(2,kend-kstart);
904:           o_nonzeros[i-kstart] = PetscMin(2,bigsize-(kend-kstart));
905:         }
906:       }
907:       MatCreate(comm,&ipmP->K);
908:       MatSetType(ipmP->K,MATMPIAIJ);
909:       MatSetSizes(ipmP->K,klocalsize,klocalsize,PETSC_DETERMINE,PETSC_DETERMINE);
910:       MatMPIAIJSetPreallocation(ipmP->K,0,d_nonzeros,0,o_nonzeros);
911:       PetscFree(d_nonzeros);
912:       PetscFree(o_nonzeros);
913:       MatSetFromOptions(ipmP->K);
914:     }
915:   }

917:   MatZeroEntries(ipmP->K);
918:   /* Copy H */
919:   for (i=hstart;i<hend;i++) {
920:     MatGetRow(tao->hessian,i,&ncols,&cols,&vals);
921:     if (ncols > 0) {
922:       MatSetValues(ipmP->K,1,&i,ncols,cols,vals,INSERT_VALUES);
923:     }
924:     MatRestoreRow(tao->hessian,i,&ncols,&cols,&vals);
925:   }

927:   /* Copy Ae and Ae' */
928:   if (ipmP->me > 0) {
929:     MatGetOwnershipRange(tao->jacobian_equality,&aestart,&aeend);
930:     for (i=aestart;i<aeend;i++) {
931:       MatGetRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
932:       if (ncols > 0) {
933:         /*Ae*/
934:         row = i+r1;
935:         MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
936:         /*Ae'*/
937:         for (j=0;j<ncols;j++) {
938:           newcol = i + c2;
939:           newrow = cols[j];
940:           newval = vals[j];
941:           MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
942:         }
943:       }
944:       MatRestoreRow(tao->jacobian_equality,i,&ncols,&cols,&vals);
945:     }
946:   }

948:   if (ipmP->nb > 0) {
949:     MatGetOwnershipRange(ipmP->Ai,&aistart,&aiend);
950:     /* Copy Ai,and Ai' */
951:     for (i=aistart;i<aiend;i++) {
952:       row = i+r2;
953:       MatGetRow(ipmP->Ai,i,&ncols,&cols,&vals);
954:       if (ncols > 0) {
955:         /*Ai*/
956:         MatSetValues(ipmP->K,1,&row,ncols,cols,vals,INSERT_VALUES);
957:         /*-Ai'*/
958:         for (j=0;j<ncols;j++) {
959:           newcol = i + c3;
960:           newrow = cols[j];
961:           newval = -vals[j];
962:           MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
963:         }
964:       }
965:       MatRestoreRow(ipmP->Ai,i,&ncols,&cols,&vals);
966:     }

968:     /* -I */
969:     for (i=kstart;i<kend;i++) {
970:       if (i>=r2 && i<r3) {
971:         newrow = i;
972:         newcol = i-r2+c1;
973:         newval = -1.0;
974:         MatSetValues(ipmP->K,1,&newrow,1,&newcol,&newval,INSERT_VALUES);
975:       }
976:     }

978:     /* Copy L,Y */
979:     VecGetOwnershipRange(ipmP->s,&sstart,&send);
980:     VecGetArrayRead(ipmP->lamdai,&l);
981:     VecGetArrayRead(ipmP->s,&y);

983:     for (i=sstart;i<send;i++) {
984:       newcols[0] = c1+i;
985:       newcols[1] = c3+i;
986:       newvals[0] = l[i-sstart];
987:       newvals[1] = y[i-sstart];
988:       newrow = r3+i;
989:       MatSetValues(ipmP->K,1,&newrow,2,newcols,newvals,INSERT_VALUES);
990:     }

992:     VecRestoreArrayRead(ipmP->lamdai,&l);
993:     VecRestoreArrayRead(ipmP->s,&y);
994:   }

996:   PetscFree(indices);
997:   PetscFree(newvals);
998:   MatAssemblyBegin(ipmP->K,MAT_FINAL_ASSEMBLY);
999:   MatAssemblyEnd(ipmP->K,MAT_FINAL_ASSEMBLY);
1000:   return(0);
1001: }

1003: PetscErrorCode IPMGatherRHS(Tao tao,Vec RHS,Vec X1,Vec X2,Vec X3,Vec X4)
1004: {
1005:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

1009:   /* rhs = [x1      (n)
1010:             x2     (me)
1011:             x3     (nb)
1012:             x4     (nb)] */
1013:   if (X1) {
1014:     VecScatterBegin(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1015:     VecScatterEnd(ipmP->rhs1,X1,RHS,INSERT_VALUES,SCATTER_FORWARD);
1016:   }
1017:   if (ipmP->me > 0 && X2) {
1018:     VecScatterBegin(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1019:     VecScatterEnd(ipmP->rhs2,X2,RHS,INSERT_VALUES,SCATTER_FORWARD);
1020:   }
1021:   if (ipmP->nb > 0) {
1022:     if (X3) {
1023:       VecScatterBegin(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1024:       VecScatterEnd(ipmP->rhs3,X3,RHS,INSERT_VALUES,SCATTER_FORWARD);
1025:     }
1026:     if (X4) {
1027:       VecScatterBegin(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1028:       VecScatterEnd(ipmP->rhs4,X4,RHS,INSERT_VALUES,SCATTER_FORWARD);
1029:     }
1030:   }
1031:   return(0);
1032: }

1034: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
1035: {
1036:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;

1040:   CHKMEMQ;
1041:   /*        [x1    (n)
1042:              x2    (nb) may be 0
1043:              x3    (me) may be 0
1044:              x4    (nb) may be 0 */
1045:   if (X1) {
1046:     VecScatterBegin(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1047:     VecScatterEnd(ipmP->step1,STEP,X1,INSERT_VALUES,SCATTER_FORWARD);
1048:   }
1049:   if (X2 && ipmP->nb > 0) {
1050:     VecScatterBegin(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1051:     VecScatterEnd(ipmP->step2,STEP,X2,INSERT_VALUES,SCATTER_FORWARD);
1052:   }
1053:   if (X3 && ipmP->me > 0) {
1054:     VecScatterBegin(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1055:     VecScatterEnd(ipmP->step3,STEP,X3,INSERT_VALUES,SCATTER_FORWARD);
1056:   }
1057:   if (X4 && ipmP->nb > 0) {
1058:     VecScatterBegin(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1059:     VecScatterEnd(ipmP->step4,STEP,X4,INSERT_VALUES,SCATTER_FORWARD);
1060:   }
1061:   CHKMEMQ;
1062:   return(0);
1063: }

1065: /*MC
1066:   TAOIPM - Interior point algorithm for generally constrained optimization.

1068:   Option Database Keys:
1069: +   -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1070: .   -tao_ipm_pushs - parameter to push initial slack variables away from bounds

1072:   Notes: 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.
1073:   Level: beginner

1075: M*/

1077: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1078: {
1079:   TAO_IPM        *ipmP;

1083:   tao->ops->setup = TaoSetup_IPM;
1084:   tao->ops->solve = TaoSolve_IPM;
1085:   tao->ops->view = TaoView_IPM;
1086:   tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1087:   tao->ops->destroy = TaoDestroy_IPM;
1088:   /* tao->ops->computedual = TaoComputeDual_IPM; */

1090:   PetscNewLog(tao,&ipmP);
1091:   tao->data = (void*)ipmP;

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

1097:   ipmP->dec = 10000; /* line search critera */
1098:   ipmP->taumin = 0.995;
1099:   ipmP->monitorkkt = PETSC_FALSE;
1100:   ipmP->pushs = 100;
1101:   ipmP->pushnu = 100;
1102:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1103:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
1104:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1105:   return(0);
1106: }