Actual source code: ipm.c

  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/constrained/impls/ipm/ipm.h>

  4: /*
  5:    x,d in R^n
  6:    f in R
  7:    nb = mi + nlb+nub
  8:    s in R^nb is slack vector CI(x) / x-XL / -x+XU
  9:    bin in R^mi (tao->constraints_inequality)
 10:    beq in R^me (tao->constraints_equality)
 11:    lambdai in R^nb (ipmP->lambdai)
 12:    lambdae in R^me (ipmP->lambdae)
 13:    Jeq in R^(me x n) (tao->jacobian_equality)
 14:    Jin in R^(mi x n) (tao->jacobian_inequality)
 15:    Ai in  R^(nb x n) (ipmP->Ai)
 16:    H in R^(n x n) (tao->hessian)
 17:    min f=(1/2)*x'*H*x + d'*x
 18:    s.t.  CE(x) == 0
 19:          CI(x) >= 0
 20:          x >= tao->XL
 21:          -x >= -tao->XU
 22: */

 24: static PetscErrorCode IPMComputeKKT(Tao tao);
 25: static PetscErrorCode IPMPushInitialPoint(Tao tao);
 26: static PetscErrorCode IPMEvaluate(Tao tao);
 27: static PetscErrorCode IPMUpdateK(Tao tao);
 28: static PetscErrorCode IPMUpdateAi(Tao tao);
 29: static PetscErrorCode IPMGatherRHS(Tao tao, Vec, Vec, Vec, Vec, Vec);
 30: static PetscErrorCode IPMScatterStep(Tao tao, Vec, Vec, Vec, Vec, Vec);
 31: static PetscErrorCode IPMInitializeBounds(Tao tao);

 33: static PetscErrorCode TaoSolve_IPM(Tao tao)
 34: {
 35:   TAO_IPM    *ipmP = (TAO_IPM *)tao->data;
 36:   PetscInt    its, i;
 37:   PetscScalar stepsize = 1.0;
 38:   PetscScalar step_s, step_l, alpha, tau, sigma, phi_target;

 40:   /* Push initial point away from bounds */
 41:   IPMInitializeBounds(tao);
 42:   IPMPushInitialPoint(tao);
 43:   VecCopy(tao->solution, ipmP->rhs_x);
 44:   IPMEvaluate(tao);
 45:   IPMComputeKKT(tao);

 47:   tao->reason = TAO_CONTINUE_ITERATING;
 48:   TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its);
 49:   TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, 1.0);
 50:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

 52:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 53:     /* Call general purpose update function */
 54:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

 56:     tao->ksp_its = 0;
 57:     IPMUpdateK(tao);
 58:     /*
 59:        rhs.x    = -rd
 60:        rhs.lame = -rpe
 61:        rhs.lami = -rpi
 62:        rhs.com  = -com
 63:     */

 65:     VecCopy(ipmP->rd, ipmP->rhs_x);
 66:     if (ipmP->me > 0) VecCopy(ipmP->rpe, ipmP->rhs_lambdae);
 67:     if (ipmP->nb > 0) {
 68:       VecCopy(ipmP->rpi, ipmP->rhs_lambdai);
 69:       VecCopy(ipmP->complementarity, ipmP->rhs_s);
 70:     }
 71:     IPMGatherRHS(tao, ipmP->bigrhs, ipmP->rhs_x, ipmP->rhs_lambdae, ipmP->rhs_lambdai, 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->dlambdae, ipmP->dlambdai);
 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->lambdai, ipmP->dlambdai, 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) VecCopy(ipmP->lambdae, ipmP->save_lambdae);
 96:     if (ipmP->nb > 0) {
 97:       VecCopy(ipmP->lambdai, ipmP->save_lambdai);
 98:       VecCopy(ipmP->s, ipmP->save_s);
 99:     }

101:     VecAXPY(tao->solution, alpha, tao->stepdirection);
102:     if (ipmP->me > 0) VecAXPY(ipmP->lambdae, alpha, ipmP->dlambdae);
103:     if (ipmP->nb > 0) {
104:       VecAXPY(ipmP->lambdai, alpha, ipmP->dlambdai);
105:       VecAXPY(ipmP->s, alpha, ipmP->ds);
106:     }

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

118:     /* revert kkt info */
119:     VecCopy(ipmP->save_x, tao->solution);
120:     if (ipmP->me > 0) VecCopy(ipmP->save_lambdae, ipmP->lambdae);
121:     if (ipmP->nb > 0) {
122:       VecCopy(ipmP->save_lambdai, ipmP->lambdai);
123:       VecCopy(ipmP->save_s, ipmP->s);
124:     }
125:     IPMComputeKKT(tao);

127:     /* update rhs with new complementarity vector */
128:     if (ipmP->nb > 0) {
129:       VecCopy(ipmP->complementarity, ipmP->rhs_s);
130:       VecScale(ipmP->rhs_s, -1.0);
131:       VecShift(ipmP->rhs_s, sigma * ipmP->mu);
132:     }
133:     IPMGatherRHS(tao, ipmP->bigrhs, NULL, NULL, NULL, ipmP->rhs_s);

135:     /* solve K * step = rhs */
136:     KSPSetOperators(tao->ksp, ipmP->K, ipmP->K);
137:     KSPSolve(tao->ksp, ipmP->bigrhs, ipmP->bigstep);

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

173:       /* update dual variables */
174:       if (ipmP->me > 0) VecCopy(ipmP->lambdae, tao->DE);

176:       IPMEvaluate(tao);
177:       IPMComputeKKT(tao);
178:       if (ipmP->phi <= phi_target) break;
179:       alpha /= 2.0;
180:     }

182:     TaoLogConvergenceHistory(tao, ipmP->kkt_f, ipmP->phi, 0.0, tao->ksp_its);
183:     TaoMonitor(tao, tao->niter, ipmP->kkt_f, ipmP->phi, 0.0, stepsize);
184:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
185:     tao->niter++;
186:   }
187:   return 0;
188: }

190: static PetscErrorCode TaoSetup_IPM(Tao tao)
191: {
192:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

194:   ipmP->nb = ipmP->mi = ipmP->me = 0;
195:   ipmP->K                        = NULL;
196:   VecGetSize(tao->solution, &ipmP->n);
197:   if (!tao->gradient) {
198:     VecDuplicate(tao->solution, &tao->gradient);
199:     VecDuplicate(tao->solution, &tao->stepdirection);
200:     VecDuplicate(tao->solution, &ipmP->rd);
201:     VecDuplicate(tao->solution, &ipmP->rhs_x);
202:     VecDuplicate(tao->solution, &ipmP->work);
203:     VecDuplicate(tao->solution, &ipmP->save_x);
204:   }
205:   if (tao->constraints_equality) {
206:     VecGetSize(tao->constraints_equality, &ipmP->me);
207:     VecDuplicate(tao->constraints_equality, &ipmP->lambdae);
208:     VecDuplicate(tao->constraints_equality, &ipmP->dlambdae);
209:     VecDuplicate(tao->constraints_equality, &ipmP->rhs_lambdae);
210:     VecDuplicate(tao->constraints_equality, &ipmP->save_lambdae);
211:     VecDuplicate(tao->constraints_equality, &ipmP->rpe);
212:     VecDuplicate(tao->constraints_equality, &tao->DE);
213:   }
214:   if (tao->constraints_inequality) VecDuplicate(tao->constraints_inequality, &tao->DI);
215:   return 0;
216: }

218: static PetscErrorCode IPMInitializeBounds(Tao tao)
219: {
220:   TAO_IPM        *ipmP = (TAO_IPM *)tao->data;
221:   Vec             xtmp;
222:   PetscInt        xstart, xend;
223:   PetscInt        ucstart, ucend;   /* user ci */
224:   PetscInt        ucestart, uceend; /* user ce */
225:   PetscInt        sstart = 0, send = 0;
226:   PetscInt        bigsize;
227:   PetscInt        i, counter, nloc;
228:   PetscInt       *cind, *xind, *ucind, *uceind, *stepind;
229:   VecType         vtype;
230:   const PetscInt *xli, *xui;
231:   PetscInt        xl_offset, xu_offset;
232:   IS              bigxl, bigxu, isuc, isc, isx, sis, is1;
233:   MPI_Comm        comm;

235:   cind = xind = ucind = uceind = stepind = NULL;
236:   ipmP->mi                               = 0;
237:   ipmP->nxlb                             = 0;
238:   ipmP->nxub                             = 0;
239:   ipmP->nb                               = 0;
240:   ipmP->nslack                           = 0;

242:   VecDuplicate(tao->solution, &xtmp);
243:   TaoComputeVariableBounds(tao);
244:   if (tao->XL) {
245:     VecSet(xtmp, PETSC_NINFINITY);
246:     VecWhichGreaterThan(tao->XL, xtmp, &ipmP->isxl);
247:     ISGetSize(ipmP->isxl, &ipmP->nxlb);
248:   } else {
249:     ipmP->nxlb = 0;
250:   }
251:   if (tao->XU) {
252:     VecSet(xtmp, PETSC_INFINITY);
253:     VecWhichLessThan(tao->XU, xtmp, &ipmP->isxu);
254:     ISGetSize(ipmP->isxu, &ipmP->nxub);
255:   } else {
256:     ipmP->nxub = 0;
257:   }
258:   VecDestroy(&xtmp);
259:   if (tao->constraints_inequality) {
260:     VecGetSize(tao->constraints_inequality, &ipmP->mi);
261:   } else {
262:     ipmP->mi = 0;
263:   }
264:   ipmP->nb = ipmP->nxlb + ipmP->nxub + ipmP->mi;

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

268:   bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
269:   PetscMalloc1(bigsize, &stepind);
270:   PetscMalloc1(ipmP->n, &xind);
271:   PetscMalloc1(ipmP->me, &uceind);
272:   VecGetOwnershipRange(tao->solution, &xstart, &xend);

274:   if (ipmP->nb > 0) {
275:     VecCreate(comm, &ipmP->s);
276:     VecSetSizes(ipmP->s, PETSC_DECIDE, ipmP->nb);
277:     VecSetFromOptions(ipmP->s);
278:     VecDuplicate(ipmP->s, &ipmP->ds);
279:     VecDuplicate(ipmP->s, &ipmP->rhs_s);
280:     VecDuplicate(ipmP->s, &ipmP->complementarity);
281:     VecDuplicate(ipmP->s, &ipmP->ci);

283:     VecDuplicate(ipmP->s, &ipmP->lambdai);
284:     VecDuplicate(ipmP->s, &ipmP->dlambdai);
285:     VecDuplicate(ipmP->s, &ipmP->rhs_lambdai);
286:     VecDuplicate(ipmP->s, &ipmP->save_lambdai);

288:     VecDuplicate(ipmP->s, &ipmP->save_s);
289:     VecDuplicate(ipmP->s, &ipmP->rpi);
290:     VecDuplicate(ipmP->s, &ipmP->Zero_nb);
291:     VecSet(ipmP->Zero_nb, 0.0);
292:     VecDuplicate(ipmP->s, &ipmP->One_nb);
293:     VecSet(ipmP->One_nb, 1.0);
294:     VecDuplicate(ipmP->s, &ipmP->Inf_nb);
295:     VecSet(ipmP->Inf_nb, PETSC_INFINITY);

297:     PetscMalloc1(ipmP->nb, &cind);
298:     PetscMalloc1(ipmP->mi, &ucind);
299:     VecGetOwnershipRange(ipmP->s, &sstart, &send);

301:     if (ipmP->mi > 0) {
302:       VecGetOwnershipRange(tao->constraints_inequality, &ucstart, &ucend);
303:       counter = 0;
304:       for (i = ucstart; i < ucend; i++) cind[counter++] = i;
305:       ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isuc);
306:       ISCreateGeneral(comm, counter, cind, PETSC_COPY_VALUES, &isc);
307:       VecScatterCreate(tao->constraints_inequality, isuc, ipmP->ci, isc, &ipmP->ci_scat);

309:       ISDestroy(&isuc);
310:       ISDestroy(&isc);
311:     }
312:     /* need to know how may xbound indices are on each process */
313:     /* TODO better way */
314:     if (ipmP->nxlb) {
315:       ISAllGather(ipmP->isxl, &bigxl);
316:       ISGetIndices(bigxl, &xli);
317:       /* find offsets for this processor */
318:       xl_offset = ipmP->mi;
319:       for (i = 0; i < ipmP->nxlb; i++) {
320:         if (xli[i] < xstart) {
321:           xl_offset++;
322:         } else break;
323:       }
324:       ISRestoreIndices(bigxl, &xli);

326:       ISGetIndices(ipmP->isxl, &xli);
327:       ISGetLocalSize(ipmP->isxl, &nloc);
328:       for (i = 0; i < nloc; i++) {
329:         xind[i] = xli[i];
330:         cind[i] = xl_offset + i;
331:       }

333:       ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx);
334:       ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc);
335:       VecScatterCreate(tao->XL, isx, ipmP->ci, isc, &ipmP->xl_scat);
336:       ISDestroy(&isx);
337:       ISDestroy(&isc);
338:       ISDestroy(&bigxl);
339:     }

341:     if (ipmP->nxub) {
342:       ISAllGather(ipmP->isxu, &bigxu);
343:       ISGetIndices(bigxu, &xui);
344:       /* find offsets for this processor */
345:       xu_offset = ipmP->mi + ipmP->nxlb;
346:       for (i = 0; i < ipmP->nxub; i++) {
347:         if (xui[i] < xstart) {
348:           xu_offset++;
349:         } else break;
350:       }
351:       ISRestoreIndices(bigxu, &xui);

353:       ISGetIndices(ipmP->isxu, &xui);
354:       ISGetLocalSize(ipmP->isxu, &nloc);
355:       for (i = 0; i < nloc; i++) {
356:         xind[i] = xui[i];
357:         cind[i] = xu_offset + i;
358:       }

360:       ISCreateGeneral(comm, nloc, xind, PETSC_COPY_VALUES, &isx);
361:       ISCreateGeneral(comm, nloc, cind, PETSC_COPY_VALUES, &isc);
362:       VecScatterCreate(tao->XU, isx, ipmP->ci, isc, &ipmP->xu_scat);
363:       ISDestroy(&isx);
364:       ISDestroy(&isc);
365:       ISDestroy(&bigxu);
366:     }
367:   }
368:   VecCreate(comm, &ipmP->bigrhs);
369:   VecGetType(tao->solution, &vtype);
370:   VecSetType(ipmP->bigrhs, vtype);
371:   VecSetSizes(ipmP->bigrhs, PETSC_DECIDE, bigsize);
372:   VecSetFromOptions(ipmP->bigrhs);
373:   VecDuplicate(ipmP->bigrhs, &ipmP->bigstep);

375:   /* create scatters for step->x and x->rhs */
376:   for (i = xstart; i < xend; i++) {
377:     stepind[i - xstart] = i;
378:     xind[i - xstart]    = i;
379:   }
380:   ISCreateGeneral(comm, xend - xstart, stepind, PETSC_COPY_VALUES, &sis);
381:   ISCreateGeneral(comm, xend - xstart, xind, PETSC_COPY_VALUES, &is1);
382:   VecScatterCreate(ipmP->bigstep, sis, tao->solution, is1, &ipmP->step1);
383:   VecScatterCreate(tao->solution, is1, ipmP->bigrhs, sis, &ipmP->rhs1);
384:   ISDestroy(&sis);
385:   ISDestroy(&is1);

387:   if (ipmP->nb > 0) {
388:     for (i = sstart; i < send; i++) {
389:       stepind[i - sstart] = i + ipmP->n;
390:       cind[i - sstart]    = i;
391:     }
392:     ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis);
393:     ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1);
394:     VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step2);
395:     ISDestroy(&sis);

397:     for (i = sstart; i < send; i++) {
398:       stepind[i - sstart] = i + ipmP->n + ipmP->me;
399:       cind[i - sstart]    = i;
400:     }
401:     ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis);
402:     VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs3);
403:     ISDestroy(&sis);
404:     ISDestroy(&is1);
405:   }

407:   if (ipmP->me > 0) {
408:     VecGetOwnershipRange(tao->constraints_equality, &ucestart, &uceend);
409:     for (i = ucestart; i < uceend; i++) {
410:       stepind[i - ucestart] = i + ipmP->n + ipmP->nb;
411:       uceind[i - ucestart]  = i;
412:     }

414:     ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis);
415:     ISCreateGeneral(comm, uceend - ucestart, uceind, PETSC_COPY_VALUES, &is1);
416:     VecScatterCreate(ipmP->bigstep, sis, tao->constraints_equality, is1, &ipmP->step3);
417:     ISDestroy(&sis);

419:     for (i = ucestart; i < uceend; i++) stepind[i - ucestart] = i + ipmP->n;

421:     ISCreateGeneral(comm, uceend - ucestart, stepind, PETSC_COPY_VALUES, &sis);
422:     VecScatterCreate(tao->constraints_equality, is1, ipmP->bigrhs, sis, &ipmP->rhs2);
423:     ISDestroy(&sis);
424:     ISDestroy(&is1);
425:   }

427:   if (ipmP->nb > 0) {
428:     for (i = sstart; i < send; i++) {
429:       stepind[i - sstart] = i + ipmP->n + ipmP->nb + ipmP->me;
430:       cind[i - sstart]    = i;
431:     }
432:     ISCreateGeneral(comm, send - sstart, cind, PETSC_COPY_VALUES, &is1);
433:     ISCreateGeneral(comm, send - sstart, stepind, PETSC_COPY_VALUES, &sis);
434:     VecScatterCreate(ipmP->bigstep, sis, ipmP->s, is1, &ipmP->step4);
435:     VecScatterCreate(ipmP->s, is1, ipmP->bigrhs, sis, &ipmP->rhs4);
436:     ISDestroy(&sis);
437:     ISDestroy(&is1);
438:   }

440:   PetscFree(stepind);
441:   PetscFree(cind);
442:   PetscFree(ucind);
443:   PetscFree(uceind);
444:   PetscFree(xind);
445:   return 0;
446: }

448: static PetscErrorCode TaoDestroy_IPM(Tao tao)
449: {
450:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

452:   VecDestroy(&ipmP->rd);
453:   VecDestroy(&ipmP->rpe);
454:   VecDestroy(&ipmP->rpi);
455:   VecDestroy(&ipmP->work);
456:   VecDestroy(&ipmP->lambdae);
457:   VecDestroy(&ipmP->lambdai);
458:   VecDestroy(&ipmP->s);
459:   VecDestroy(&ipmP->ds);
460:   VecDestroy(&ipmP->ci);

462:   VecDestroy(&ipmP->rhs_x);
463:   VecDestroy(&ipmP->rhs_lambdae);
464:   VecDestroy(&ipmP->rhs_lambdai);
465:   VecDestroy(&ipmP->rhs_s);

467:   VecDestroy(&ipmP->save_x);
468:   VecDestroy(&ipmP->save_lambdae);
469:   VecDestroy(&ipmP->save_lambdai);
470:   VecDestroy(&ipmP->save_s);

472:   VecScatterDestroy(&ipmP->step1);
473:   VecScatterDestroy(&ipmP->step2);
474:   VecScatterDestroy(&ipmP->step3);
475:   VecScatterDestroy(&ipmP->step4);

477:   VecScatterDestroy(&ipmP->rhs1);
478:   VecScatterDestroy(&ipmP->rhs2);
479:   VecScatterDestroy(&ipmP->rhs3);
480:   VecScatterDestroy(&ipmP->rhs4);

482:   VecScatterDestroy(&ipmP->ci_scat);
483:   VecScatterDestroy(&ipmP->xl_scat);
484:   VecScatterDestroy(&ipmP->xu_scat);

486:   VecDestroy(&ipmP->dlambdai);
487:   VecDestroy(&ipmP->dlambdae);
488:   VecDestroy(&ipmP->Zero_nb);
489:   VecDestroy(&ipmP->One_nb);
490:   VecDestroy(&ipmP->Inf_nb);
491:   VecDestroy(&ipmP->complementarity);

493:   VecDestroy(&ipmP->bigrhs);
494:   VecDestroy(&ipmP->bigstep);
495:   MatDestroy(&ipmP->Ai);
496:   MatDestroy(&ipmP->K);
497:   ISDestroy(&ipmP->isxu);
498:   ISDestroy(&ipmP->isxl);
499:   KSPDestroy(&tao->ksp);
500:   PetscFree(tao->data);
501:   return 0;
502: }

504: static PetscErrorCode TaoSetFromOptions_IPM(Tao tao, PetscOptionItems *PetscOptionsObject)
505: {
506:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

508:   PetscOptionsHeadBegin(PetscOptionsObject, "IPM method for constrained optimization");
509:   PetscOptionsBool("-tao_ipm_monitorkkt", "monitor kkt status", NULL, ipmP->monitorkkt, &ipmP->monitorkkt, NULL);
510:   PetscOptionsReal("-tao_ipm_pushs", "parameter to push initial slack variables away from bounds", NULL, ipmP->pushs, &ipmP->pushs, NULL);
511:   PetscOptionsReal("-tao_ipm_pushnu", "parameter to push initial (inequality) dual variables away from bounds", NULL, ipmP->pushnu, &ipmP->pushnu, NULL);
512:   PetscOptionsHeadEnd();
513:   KSPSetFromOptions(tao->ksp);
514:   return 0;
515: }

517: static PetscErrorCode TaoView_IPM(Tao tao, PetscViewer viewer)
518: {
519:   return 0;
520: }

522: /* IPMObjectiveAndGradient()
523:    f = d'x + 0.5 * x' * H * x
524:    rd = H*x + d + Ae'*lame - Ai'*lami
525:    rpe = Ae*x - be
526:    rpi = Ai*x - yi - bi
527:    mu = yi' * lami/mi;
528:    com = yi.*lami

530:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
531: */
532: /*
533: static PetscErrorCode IPMObjective(TaoLineSearch ls, Vec X, PetscReal *f, void *tptr)
534: {
535:   Tao tao = (Tao)tptr;
536:   TAO_IPM *ipmP = (TAO_IPM*)tao->data;
537:   IPMComputeKKT(tao);
538:   *f = ipmP->phi;
539:   return 0;
540: }
541: */

543: /*
544:    f = d'x + 0.5 * x' * H * x
545:    rd = H*x + d + Ae'*lame - Ai'*lami
546:        Ai =   jac_ineq
547:                I (w/lb)
548:               -I (w/ub)

550:    rpe = ce
551:    rpi = ci - s;
552:    com = s.*lami
553:    mu = yi' * lami/mi;

555:    phi = ||rd|| + ||rpe|| + ||rpi|| + ||com||
556: */
557: static PetscErrorCode IPMComputeKKT(Tao tao)
558: {
559:   TAO_IPM    *ipmP = (TAO_IPM *)tao->data;
560:   PetscScalar norm;

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

564:   if (ipmP->me > 0) {
565:     /* rd = gradient + Ae'*lambdae */
566:     MatMultTranspose(tao->jacobian_equality, ipmP->lambdae, ipmP->work);
567:     VecAXPY(ipmP->rd, 1.0, ipmP->work);

569:     /* rpe = ce(x) */
570:     VecCopy(tao->constraints_equality, ipmP->rpe);
571:   }
572:   if (ipmP->nb > 0) {
573:     /* rd = rd - Ai'*lambdai */
574:     MatMultTranspose(ipmP->Ai, ipmP->lambdai, ipmP->work);
575:     VecAXPY(ipmP->rd, -1.0, ipmP->work);

577:     /* rpi = cin - s */
578:     VecCopy(ipmP->ci, ipmP->rpi);
579:     VecAXPY(ipmP->rpi, -1.0, ipmP->s);

581:     /* com = s .* lami */
582:     VecPointwiseMult(ipmP->complementarity, ipmP->s, ipmP->lambdai);
583:   }
584:   /* phi = ||rd; rpe; rpi; com|| */
585:   VecDot(ipmP->rd, ipmP->rd, &norm);
586:   ipmP->phi = norm;
587:   if (ipmP->me > 0) {
588:     VecDot(ipmP->rpe, ipmP->rpe, &norm);
589:     ipmP->phi += norm;
590:   }
591:   if (ipmP->nb > 0) {
592:     VecDot(ipmP->rpi, ipmP->rpi, &norm);
593:     ipmP->phi += norm;
594:     VecDot(ipmP->complementarity, ipmP->complementarity, &norm);
595:     ipmP->phi += norm;
596:     /* mu = s'*lami/nb */
597:     VecDot(ipmP->s, ipmP->lambdai, &ipmP->mu);
598:     ipmP->mu /= ipmP->nb;
599:   } else {
600:     ipmP->mu = 1.0;
601:   }

603:   ipmP->phi = PetscSqrtScalar(ipmP->phi);
604:   return 0;
605: }

607: /* evaluate user info at current point */
608: PetscErrorCode IPMEvaluate(Tao tao)
609: {
610:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

612:   TaoComputeObjectiveAndGradient(tao, tao->solution, &ipmP->kkt_f, tao->gradient);
613:   TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
614:   if (ipmP->me > 0) {
615:     TaoComputeEqualityConstraints(tao, tao->solution, tao->constraints_equality);
616:     TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre);
617:   }
618:   if (ipmP->mi > 0) {
619:     TaoComputeInequalityConstraints(tao, tao->solution, tao->constraints_inequality);
620:     TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre);
621:   }
622:   if (ipmP->nb > 0) {
623:     /* Ai' =   jac_ineq | I (w/lb) | -I (w/ub)  */
624:     IPMUpdateAi(tao);
625:   }
626:   return 0;
627: }

629: /* Push initial point away from bounds */
630: PetscErrorCode IPMPushInitialPoint(Tao tao)
631: {
632:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

634:   TaoComputeVariableBounds(tao);
635:   VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
636:   if (ipmP->nb > 0) {
637:     VecSet(ipmP->s, ipmP->pushs);
638:     VecSet(ipmP->lambdai, ipmP->pushnu);
639:     if (ipmP->mi > 0) VecSet(tao->DI, ipmP->pushnu);
640:   }
641:   if (ipmP->me > 0) {
642:     VecSet(tao->DE, 1.0);
643:     VecSet(ipmP->lambdae, 1.0);
644:   }
645:   return 0;
646: }

648: PetscErrorCode IPMUpdateAi(Tao tao)
649: {
650:   /* Ai =     Ji
651:               I (w/lb)
652:              -I (w/ub) */

654:   /* Ci =    user->ci
655:              Xi - lb (w/lb)
656:              -Xi + ub (w/ub)  */

658:   TAO_IPM           *ipmP = (TAO_IPM *)tao->data;
659:   MPI_Comm           comm;
660:   PetscInt           i;
661:   PetscScalar        newval;
662:   PetscInt           newrow, newcol, ncols;
663:   const PetscScalar *vals;
664:   const PetscInt    *cols;
665:   PetscInt           astart, aend, jstart, jend;
666:   PetscInt          *nonzeros;
667:   PetscInt           r2, r3, r4;
668:   PetscMPIInt        size;
669:   Vec                solu;
670:   PetscInt           nloc;

672:   r2 = ipmP->mi;
673:   r3 = r2 + ipmP->nxlb;
674:   r4 = r3 + ipmP->nxub;

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

678:   /* Create Ai matrix if it doesn't exist yet */
679:   if (!ipmP->Ai) {
680:     comm = ((PetscObject)(tao->solution))->comm;
681:     MPI_Comm_size(comm, &size);
682:     if (size == 1) {
683:       PetscMalloc1(ipmP->nb, &nonzeros);
684:       for (i = 0; i < ipmP->mi; i++) {
685:         MatGetRow(tao->jacobian_inequality, i, &ncols, NULL, NULL);
686:         nonzeros[i] = ncols;
687:         MatRestoreRow(tao->jacobian_inequality, i, &ncols, NULL, NULL);
688:       }
689:       for (i = r2; i < r4; i++) nonzeros[i] = 1;
690:     }
691:     MatCreate(comm, &ipmP->Ai);
692:     MatSetType(ipmP->Ai, MATAIJ);

694:     TaoGetSolution(tao, &solu);
695:     VecGetLocalSize(solu, &nloc);
696:     MatSetSizes(ipmP->Ai, PETSC_DECIDE, nloc, ipmP->nb, PETSC_DECIDE);
697:     MatSetFromOptions(ipmP->Ai);
698:     MatMPIAIJSetPreallocation(ipmP->Ai, ipmP->nb, NULL, ipmP->nb, NULL);
699:     MatSeqAIJSetPreallocation(ipmP->Ai, PETSC_DEFAULT, nonzeros);
700:     if (size == 1) PetscFree(nonzeros);
701:   }

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

706:   /* Ai w/lb */
707:   if (ipmP->mi) {
708:     MatZeroEntries(ipmP->Ai);
709:     MatGetOwnershipRange(tao->jacobian_inequality, &jstart, &jend);
710:     for (i = jstart; i < jend; i++) {
711:       MatGetRow(tao->jacobian_inequality, i, &ncols, &cols, &vals);
712:       newrow = i;
713:       MatSetValues(ipmP->Ai, 1, &newrow, ncols, cols, vals, INSERT_VALUES);
714:       MatRestoreRow(tao->jacobian_inequality, i, &ncols, &cols, &vals);
715:     }
716:   }

718:   /* I w/ xlb */
719:   if (ipmP->nxlb) {
720:     for (i = 0; i < ipmP->nxlb; i++) {
721:       if (i >= astart && i < aend) {
722:         newrow = i + r2;
723:         newcol = i;
724:         newval = 1.0;
725:         MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES);
726:       }
727:     }
728:   }
729:   if (ipmP->nxub) {
730:     /* I w/ xub */
731:     for (i = 0; i < ipmP->nxub; i++) {
732:       if (i >= astart && i < aend) {
733:         newrow = i + r3;
734:         newcol = i;
735:         newval = -1.0;
736:         MatSetValues(ipmP->Ai, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES);
737:       }
738:     }
739:   }

741:   MatAssemblyBegin(ipmP->Ai, MAT_FINAL_ASSEMBLY);
742:   MatAssemblyEnd(ipmP->Ai, MAT_FINAL_ASSEMBLY);
743:   CHKMEMQ;

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

747:   /* user ci */
748:   if (ipmP->mi > 0) {
749:     VecScatterBegin(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
750:     VecScatterEnd(ipmP->ci_scat, tao->constraints_inequality, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
751:   }
752:   if (!ipmP->work) VecDuplicate(tao->solution, &ipmP->work);
753:   VecCopy(tao->solution, ipmP->work);
754:   if (tao->XL) {
755:     VecAXPY(ipmP->work, -1.0, tao->XL);

757:     /* lower bounds on variables */
758:     if (ipmP->nxlb > 0) {
759:       VecScatterBegin(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
760:       VecScatterEnd(ipmP->xl_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
761:     }
762:   }
763:   if (tao->XU) {
764:     /* upper bounds on variables */
765:     VecCopy(tao->solution, ipmP->work);
766:     VecScale(ipmP->work, -1.0);
767:     VecAXPY(ipmP->work, 1.0, tao->XU);
768:     if (ipmP->nxub > 0) {
769:       VecScatterBegin(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
770:       VecScatterEnd(ipmP->xu_scat, ipmP->work, ipmP->ci, INSERT_VALUES, SCATTER_FORWARD);
771:     }
772:   }
773:   return 0;
774: }

776: /* create K = [ Hlag , 0 , Ae', -Ai'];
777:               [Ae , 0,   0  , 0];
778:               [Ai ,-I,   0 ,  0];
779:               [ 0 , S ,  0,   Y ];  */
780: PetscErrorCode IPMUpdateK(Tao tao)
781: {
782:   TAO_IPM         *ipmP = (TAO_IPM *)tao->data;
783:   MPI_Comm         comm;
784:   PetscMPIInt      size;
785:   PetscInt         i, j, row;
786:   PetscInt         ncols, newcol, newcols[2], newrow;
787:   const PetscInt  *cols;
788:   const PetscReal *vals;
789:   const PetscReal *l, *y;
790:   PetscReal       *newvals;
791:   PetscReal        newval;
792:   PetscInt         subsize;
793:   const PetscInt  *indices;
794:   PetscInt        *nonzeros, *d_nonzeros, *o_nonzeros;
795:   PetscInt         bigsize;
796:   PetscInt         r1, r2, r3;
797:   PetscInt         c1, c2, c3;
798:   PetscInt         klocalsize;
799:   PetscInt         hstart, hend, kstart, kend;
800:   PetscInt         aistart, aiend, aestart, aeend;
801:   PetscInt         sstart, send;

803:   comm = ((PetscObject)(tao->solution))->comm;
804:   MPI_Comm_size(comm, &size);
805:   IPMUpdateAi(tao);

807:   /* allocate workspace */
808:   subsize = PetscMax(ipmP->n, ipmP->nb);
809:   subsize = PetscMax(ipmP->me, subsize);
810:   subsize = PetscMax(2, subsize);
811:   PetscMalloc1(subsize, (PetscInt **)&indices);
812:   PetscMalloc1(subsize, &newvals);

814:   r1 = c1 = ipmP->n;
815:   r2      = r1 + ipmP->me;
816:   c2      = c1 + ipmP->nb;
817:   r3 = c3 = r2 + ipmP->nb;

819:   bigsize = ipmP->n + 2 * ipmP->nb + ipmP->me;
820:   VecGetOwnershipRange(ipmP->bigrhs, &kstart, &kend);
821:   MatGetOwnershipRange(tao->hessian, &hstart, &hend);
822:   klocalsize = kend - kstart;
823:   if (!ipmP->K) {
824:     if (size == 1) {
825:       PetscMalloc1(kend - kstart, &nonzeros);
826:       for (i = 0; i < bigsize; i++) {
827:         if (i < r1) {
828:           MatGetRow(tao->hessian, i, &ncols, NULL, NULL);
829:           nonzeros[i] = ncols;
830:           MatRestoreRow(tao->hessian, i, &ncols, NULL, NULL);
831:           nonzeros[i] += ipmP->me + ipmP->nb;
832:         } else if (i < r2) {
833:           nonzeros[i - kstart] = ipmP->n;
834:         } else if (i < r3) {
835:           nonzeros[i - kstart] = ipmP->n + 1;
836:         } else if (i < bigsize) {
837:           nonzeros[i - kstart] = 2;
838:         }
839:       }
840:       MatCreate(comm, &ipmP->K);
841:       MatSetType(ipmP->K, MATSEQAIJ);
842:       MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE);
843:       MatSeqAIJSetPreallocation(ipmP->K, 0, nonzeros);
844:       MatSetFromOptions(ipmP->K);
845:       PetscFree(nonzeros);
846:     } else {
847:       PetscMalloc1(kend - kstart, &d_nonzeros);
848:       PetscMalloc1(kend - kstart, &o_nonzeros);
849:       for (i = kstart; i < kend; i++) {
850:         if (i < r1) {
851:           /* TODO fix preallocation for mpi mats */
852:           d_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, kend - kstart);
853:           o_nonzeros[i - kstart] = PetscMin(ipmP->n + ipmP->me + ipmP->nb, bigsize - (kend - kstart));
854:         } else if (i < r2) {
855:           d_nonzeros[i - kstart] = PetscMin(ipmP->n, kend - kstart);
856:           o_nonzeros[i - kstart] = PetscMin(ipmP->n, bigsize - (kend - kstart));
857:         } else if (i < r3) {
858:           d_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, kend - kstart);
859:           o_nonzeros[i - kstart] = PetscMin(ipmP->n + 2, bigsize - (kend - kstart));
860:         } else {
861:           d_nonzeros[i - kstart] = PetscMin(2, kend - kstart);
862:           o_nonzeros[i - kstart] = PetscMin(2, bigsize - (kend - kstart));
863:         }
864:       }
865:       MatCreate(comm, &ipmP->K);
866:       MatSetType(ipmP->K, MATMPIAIJ);
867:       MatSetSizes(ipmP->K, klocalsize, klocalsize, PETSC_DETERMINE, PETSC_DETERMINE);
868:       MatMPIAIJSetPreallocation(ipmP->K, 0, d_nonzeros, 0, o_nonzeros);
869:       PetscFree(d_nonzeros);
870:       PetscFree(o_nonzeros);
871:       MatSetFromOptions(ipmP->K);
872:     }
873:   }

875:   MatZeroEntries(ipmP->K);
876:   /* Copy H */
877:   for (i = hstart; i < hend; i++) {
878:     MatGetRow(tao->hessian, i, &ncols, &cols, &vals);
879:     if (ncols > 0) MatSetValues(ipmP->K, 1, &i, ncols, cols, vals, INSERT_VALUES);
880:     MatRestoreRow(tao->hessian, i, &ncols, &cols, &vals);
881:   }

883:   /* Copy Ae and Ae' */
884:   if (ipmP->me > 0) {
885:     MatGetOwnershipRange(tao->jacobian_equality, &aestart, &aeend);
886:     for (i = aestart; i < aeend; i++) {
887:       MatGetRow(tao->jacobian_equality, i, &ncols, &cols, &vals);
888:       if (ncols > 0) {
889:         /*Ae*/
890:         row = i + r1;
891:         MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES);
892:         /*Ae'*/
893:         for (j = 0; j < ncols; j++) {
894:           newcol = i + c2;
895:           newrow = cols[j];
896:           newval = vals[j];
897:           MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES);
898:         }
899:       }
900:       MatRestoreRow(tao->jacobian_equality, i, &ncols, &cols, &vals);
901:     }
902:   }

904:   if (ipmP->nb > 0) {
905:     MatGetOwnershipRange(ipmP->Ai, &aistart, &aiend);
906:     /* Copy Ai,and Ai' */
907:     for (i = aistart; i < aiend; i++) {
908:       row = i + r2;
909:       MatGetRow(ipmP->Ai, i, &ncols, &cols, &vals);
910:       if (ncols > 0) {
911:         /*Ai*/
912:         MatSetValues(ipmP->K, 1, &row, ncols, cols, vals, INSERT_VALUES);
913:         /*-Ai'*/
914:         for (j = 0; j < ncols; j++) {
915:           newcol = i + c3;
916:           newrow = cols[j];
917:           newval = -vals[j];
918:           MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES);
919:         }
920:       }
921:       MatRestoreRow(ipmP->Ai, i, &ncols, &cols, &vals);
922:     }

924:     /* -I */
925:     for (i = kstart; i < kend; i++) {
926:       if (i >= r2 && i < r3) {
927:         newrow = i;
928:         newcol = i - r2 + c1;
929:         newval = -1.0;
930:         MatSetValues(ipmP->K, 1, &newrow, 1, &newcol, &newval, INSERT_VALUES);
931:       }
932:     }

934:     /* Copy L,Y */
935:     VecGetOwnershipRange(ipmP->s, &sstart, &send);
936:     VecGetArrayRead(ipmP->lambdai, &l);
937:     VecGetArrayRead(ipmP->s, &y);

939:     for (i = sstart; i < send; i++) {
940:       newcols[0] = c1 + i;
941:       newcols[1] = c3 + i;
942:       newvals[0] = l[i - sstart];
943:       newvals[1] = y[i - sstart];
944:       newrow     = r3 + i;
945:       MatSetValues(ipmP->K, 1, &newrow, 2, newcols, newvals, INSERT_VALUES);
946:     }

948:     VecRestoreArrayRead(ipmP->lambdai, &l);
949:     VecRestoreArrayRead(ipmP->s, &y);
950:   }

952:   PetscFree(indices);
953:   PetscFree(newvals);
954:   MatAssemblyBegin(ipmP->K, MAT_FINAL_ASSEMBLY);
955:   MatAssemblyEnd(ipmP->K, MAT_FINAL_ASSEMBLY);
956:   return 0;
957: }

959: PetscErrorCode IPMGatherRHS(Tao tao, Vec RHS, Vec X1, Vec X2, Vec X3, Vec X4)
960: {
961:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

963:   /* rhs = [x1      (n)
964:             x2     (me)
965:             x3     (nb)
966:             x4     (nb)] */
967:   if (X1) {
968:     VecScatterBegin(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD);
969:     VecScatterEnd(ipmP->rhs1, X1, RHS, INSERT_VALUES, SCATTER_FORWARD);
970:   }
971:   if (ipmP->me > 0 && X2) {
972:     VecScatterBegin(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD);
973:     VecScatterEnd(ipmP->rhs2, X2, RHS, INSERT_VALUES, SCATTER_FORWARD);
974:   }
975:   if (ipmP->nb > 0) {
976:     if (X3) {
977:       VecScatterBegin(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD);
978:       VecScatterEnd(ipmP->rhs3, X3, RHS, INSERT_VALUES, SCATTER_FORWARD);
979:     }
980:     if (X4) {
981:       VecScatterBegin(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD);
982:       VecScatterEnd(ipmP->rhs4, X4, RHS, INSERT_VALUES, SCATTER_FORWARD);
983:     }
984:   }
985:   return 0;
986: }

988: PetscErrorCode IPMScatterStep(Tao tao, Vec STEP, Vec X1, Vec X2, Vec X3, Vec X4)
989: {
990:   TAO_IPM *ipmP = (TAO_IPM *)tao->data;

992:   CHKMEMQ;
993:   /*        [x1    (n)
994:              x2    (nb) may be 0
995:              x3    (me) may be 0
996:              x4    (nb) may be 0 */
997:   if (X1) {
998:     VecScatterBegin(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD);
999:     VecScatterEnd(ipmP->step1, STEP, X1, INSERT_VALUES, SCATTER_FORWARD);
1000:   }
1001:   if (X2 && ipmP->nb > 0) {
1002:     VecScatterBegin(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD);
1003:     VecScatterEnd(ipmP->step2, STEP, X2, INSERT_VALUES, SCATTER_FORWARD);
1004:   }
1005:   if (X3 && ipmP->me > 0) {
1006:     VecScatterBegin(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD);
1007:     VecScatterEnd(ipmP->step3, STEP, X3, INSERT_VALUES, SCATTER_FORWARD);
1008:   }
1009:   if (X4 && ipmP->nb > 0) {
1010:     VecScatterBegin(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD);
1011:     VecScatterEnd(ipmP->step4, STEP, X4, INSERT_VALUES, SCATTER_FORWARD);
1012:   }
1013:   CHKMEMQ;
1014:   return 0;
1015: }

1017: /*MC
1018:   TAOIPM - Interior point algorithm for generally constrained optimization.

1020:   Option Database Keys:
1021: +   -tao_ipm_pushnu - parameter to push initial dual variables away from bounds
1022: -   -tao_ipm_pushs - parameter to push initial slack variables away from bounds

1024:   Notes:
1025:     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.
1026:   Level: beginner

1028: M*/

1030: PETSC_EXTERN PetscErrorCode TaoCreate_IPM(Tao tao)
1031: {
1032:   TAO_IPM *ipmP;

1034:   tao->ops->setup          = TaoSetup_IPM;
1035:   tao->ops->solve          = TaoSolve_IPM;
1036:   tao->ops->view           = TaoView_IPM;
1037:   tao->ops->setfromoptions = TaoSetFromOptions_IPM;
1038:   tao->ops->destroy        = TaoDestroy_IPM;
1039:   /* tao->ops->computedual = TaoComputeDual_IPM; */

1041:   PetscNew(&ipmP);
1042:   tao->data = (void *)ipmP;

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

1048:   ipmP->dec        = 10000; /* line search criteria */
1049:   ipmP->taumin     = 0.995;
1050:   ipmP->monitorkkt = PETSC_FALSE;
1051:   ipmP->pushs      = 100;
1052:   ipmP->pushnu     = 100;
1053:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1054:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
1055:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
1056:   return 0;
1057: }