Actual source code: pdipm.c

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

  3: /*
  4:    TaoPDIPMEvaluateFunctionsAndJacobians - Evaluate the objective function f, gradient fx, constraints, and all the Jacobians at current vector

  6:    Collective

  8:    Input Parameter:
  9: +  tao - solver context
 10: -  x - vector at which all objects to be evaluated

 12:    Level: beginner

 14: .seealso: `TaoPDIPMUpdateConstraints()`, `TaoPDIPMSetUpBounds()`
 15: */
 16: static PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao, Vec x)
 17: {
 18:   TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;

 20:   /* Compute user objective function and gradient */
 21:   TaoComputeObjectiveAndGradient(tao, x, &pdipm->obj, tao->gradient);

 23:   /* Equality constraints and Jacobian */
 24:   if (pdipm->Ng) {
 25:     TaoComputeEqualityConstraints(tao, x, tao->constraints_equality);
 26:     TaoComputeJacobianEquality(tao, x, tao->jacobian_equality, tao->jacobian_equality_pre);
 27:   }

 29:   /* Inequality constraints and Jacobian */
 30:   if (pdipm->Nh) {
 31:     TaoComputeInequalityConstraints(tao, x, tao->constraints_inequality);
 32:     TaoComputeJacobianInequality(tao, x, tao->jacobian_inequality, tao->jacobian_inequality_pre);
 33:   }
 34:   return 0;
 35: }

 37: /*
 38:   TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x

 40:   Collective

 42:   Input Parameter:
 43: + tao - Tao context
 44: - x - vector at which constraints to be evaluated

 46:    Level: beginner

 48: .seealso: `TaoPDIPMEvaluateFunctionsAndJacobians()`
 49: */
 50: static PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao, Vec x)
 51: {
 52:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
 53:   PetscInt           i, offset, offset1, k, xstart;
 54:   PetscScalar       *carr;
 55:   const PetscInt    *ubptr, *lbptr, *bxptr, *fxptr;
 56:   const PetscScalar *xarr, *xuarr, *xlarr, *garr, *harr;

 58:   VecGetOwnershipRange(x, &xstart, NULL);

 60:   VecGetArrayRead(x, &xarr);
 61:   VecGetArrayRead(tao->XU, &xuarr);
 62:   VecGetArrayRead(tao->XL, &xlarr);

 64:   /* (1) Update ce vector */
 65:   VecGetArrayWrite(pdipm->ce, &carr);

 67:   if (pdipm->Ng) {
 68:     /* (1.a) Inserting updated g(x) */
 69:     VecGetArrayRead(tao->constraints_equality, &garr);
 70:     PetscMemcpy(carr, garr, pdipm->ng * sizeof(PetscScalar));
 71:     VecRestoreArrayRead(tao->constraints_equality, &garr);
 72:   }

 74:   /* (1.b) Update xfixed */
 75:   if (pdipm->Nxfixed) {
 76:     offset = pdipm->ng;
 77:     ISGetIndices(pdipm->isxfixed, &fxptr); /* global indices in x */
 78:     for (k = 0; k < pdipm->nxfixed; k++) {
 79:       i                = fxptr[k] - xstart;
 80:       carr[offset + k] = xarr[i] - xuarr[i];
 81:     }
 82:   }
 83:   VecRestoreArrayWrite(pdipm->ce, &carr);

 85:   /* (2) Update ci vector */
 86:   VecGetArrayWrite(pdipm->ci, &carr);

 88:   if (pdipm->Nh) {
 89:     /* (2.a) Inserting updated h(x) */
 90:     VecGetArrayRead(tao->constraints_inequality, &harr);
 91:     PetscMemcpy(carr, harr, pdipm->nh * sizeof(PetscScalar));
 92:     VecRestoreArrayRead(tao->constraints_inequality, &harr);
 93:   }

 95:   /* (2.b) Update xub */
 96:   offset = pdipm->nh;
 97:   if (pdipm->Nxub) {
 98:     ISGetIndices(pdipm->isxub, &ubptr);
 99:     for (k = 0; k < pdipm->nxub; k++) {
100:       i                = ubptr[k] - xstart;
101:       carr[offset + k] = xuarr[i] - xarr[i];
102:     }
103:   }

105:   if (pdipm->Nxlb) {
106:     /* (2.c) Update xlb */
107:     offset += pdipm->nxub;
108:     ISGetIndices(pdipm->isxlb, &lbptr); /* global indices in x */
109:     for (k = 0; k < pdipm->nxlb; k++) {
110:       i                = lbptr[k] - xstart;
111:       carr[offset + k] = xarr[i] - xlarr[i];
112:     }
113:   }

115:   if (pdipm->Nxbox) {
116:     /* (2.d) Update xbox */
117:     offset += pdipm->nxlb;
118:     offset1 = offset + pdipm->nxbox;
119:     ISGetIndices(pdipm->isxbox, &bxptr); /* global indices in x */
120:     for (k = 0; k < pdipm->nxbox; k++) {
121:       i                 = bxptr[k] - xstart; /* local indices in x */
122:       carr[offset + k]  = xuarr[i] - xarr[i];
123:       carr[offset1 + k] = xarr[i] - xlarr[i];
124:     }
125:   }
126:   VecRestoreArrayWrite(pdipm->ci, &carr);

128:   /* Restoring Vectors */
129:   VecRestoreArrayRead(x, &xarr);
130:   VecRestoreArrayRead(tao->XU, &xuarr);
131:   VecRestoreArrayRead(tao->XL, &xlarr);
132:   return 0;
133: }

135: /*
136:    TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x

138:    Collective

140:    Input Parameter:
141: .  tao - holds pdipm and XL & XU

143:    Level: beginner

145: .seealso: `TaoPDIPMUpdateConstraints`
146: */
147: static PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
148: {
149:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
150:   const PetscScalar *xl, *xu;
151:   PetscInt           n, *ixlb, *ixub, *ixfixed, *ixfree, *ixbox, i, low, high, idx;
152:   MPI_Comm           comm;
153:   PetscInt           sendbuf[5], recvbuf[5];

155:   /* Creates upper and lower bounds vectors on x, if not created already */
156:   TaoComputeVariableBounds(tao);

158:   VecGetLocalSize(tao->XL, &n);
159:   PetscMalloc5(n, &ixlb, n, &ixub, n, &ixfree, n, &ixfixed, n, &ixbox);

161:   VecGetOwnershipRange(tao->XL, &low, &high);
162:   VecGetArrayRead(tao->XL, &xl);
163:   VecGetArrayRead(tao->XU, &xu);
164:   for (i = 0; i < n; i++) {
165:     idx = low + i;
166:     if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
167:       if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
168:         ixfixed[pdipm->nxfixed++] = idx;
169:       } else ixbox[pdipm->nxbox++] = idx;
170:     } else {
171:       if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
172:         ixlb[pdipm->nxlb++] = idx;
173:       } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
174:         ixub[pdipm->nxlb++] = idx;
175:       } else ixfree[pdipm->nxfree++] = idx;
176:     }
177:   }
178:   VecRestoreArrayRead(tao->XL, &xl);
179:   VecRestoreArrayRead(tao->XU, &xu);

181:   PetscObjectGetComm((PetscObject)tao, &comm);
182:   sendbuf[0] = pdipm->nxlb;
183:   sendbuf[1] = pdipm->nxub;
184:   sendbuf[2] = pdipm->nxfixed;
185:   sendbuf[3] = pdipm->nxbox;
186:   sendbuf[4] = pdipm->nxfree;

188:   MPI_Allreduce(sendbuf, recvbuf, 5, MPIU_INT, MPI_SUM, comm);
189:   pdipm->Nxlb    = recvbuf[0];
190:   pdipm->Nxub    = recvbuf[1];
191:   pdipm->Nxfixed = recvbuf[2];
192:   pdipm->Nxbox   = recvbuf[3];
193:   pdipm->Nxfree  = recvbuf[4];

195:   if (pdipm->Nxlb) ISCreateGeneral(comm, pdipm->nxlb, ixlb, PETSC_COPY_VALUES, &pdipm->isxlb);
196:   if (pdipm->Nxub) ISCreateGeneral(comm, pdipm->nxub, ixub, PETSC_COPY_VALUES, &pdipm->isxub);
197:   if (pdipm->Nxfixed) ISCreateGeneral(comm, pdipm->nxfixed, ixfixed, PETSC_COPY_VALUES, &pdipm->isxfixed);
198:   if (pdipm->Nxbox) ISCreateGeneral(comm, pdipm->nxbox, ixbox, PETSC_COPY_VALUES, &pdipm->isxbox);
199:   if (pdipm->Nxfree) ISCreateGeneral(comm, pdipm->nxfree, ixfree, PETSC_COPY_VALUES, &pdipm->isxfree);
200:   PetscFree5(ixlb, ixub, ixfixed, ixbox, ixfree);
201:   return 0;
202: }

204: /*
205:    TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
206:    X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
207:      four subvectors need to be initialized and its values copied over to X. Instead
208:      of copying, we use VecPlace/ResetArray functions to share the memory locations for
209:      X and the subvectors

211:    Collective

213:    Input Parameter:
214: .  tao - Tao context

216:    Level: beginner
217: */
218: static PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
219: {
220:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
221:   PetscScalar       *Xarr, *z, *lambdai;
222:   PetscInt           i;
223:   const PetscScalar *xarr, *h;

225:   VecGetArrayWrite(pdipm->X, &Xarr);

227:   /* Set Initialize X.x = tao->solution */
228:   VecGetArrayRead(tao->solution, &xarr);
229:   PetscMemcpy(Xarr, xarr, pdipm->nx * sizeof(PetscScalar));
230:   VecRestoreArrayRead(tao->solution, &xarr);

232:   /* Initialize X.lambdae = 0.0 */
233:   if (pdipm->lambdae) VecSet(pdipm->lambdae, 0.0);

235:   /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
236:   if (pdipm->Nci) {
237:     VecSet(pdipm->lambdai, pdipm->push_init_lambdai);
238:     VecSet(pdipm->z, pdipm->push_init_slack);

240:     /* Additional modification for X.lambdai and X.z */
241:     VecGetArrayWrite(pdipm->lambdai, &lambdai);
242:     VecGetArrayWrite(pdipm->z, &z);
243:     if (pdipm->Nh) {
244:       VecGetArrayRead(tao->constraints_inequality, &h);
245:       for (i = 0; i < pdipm->nh; i++) {
246:         if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
247:         if (pdipm->mu / z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu / z[i];
248:       }
249:       VecRestoreArrayRead(tao->constraints_inequality, &h);
250:     }
251:     VecRestoreArrayWrite(pdipm->lambdai, &lambdai);
252:     VecRestoreArrayWrite(pdipm->z, &z);
253:   }

255:   VecRestoreArrayWrite(pdipm->X, &Xarr);
256:   return 0;
257: }

259: /*
260:    TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X

262:    Input Parameter:
263:    snes - SNES context
264:    X - KKT Vector
265:    *ctx - pdipm context

267:    Output Parameter:
268:    J - Hessian matrix
269:    Jpre - Preconditioner
270: */
271: static PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx)
272: {
273:   Tao                tao   = (Tao)ctx;
274:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
275:   PetscInt           i, row, cols[2], Jrstart, rjstart, nc, j;
276:   const PetscInt    *aj, *ranges, *Jranges, *rranges, *cranges;
277:   const PetscScalar *Xarr, *aa;
278:   PetscScalar        vals[2];
279:   PetscInt           proc, nx_all, *nce_all = pdipm->nce_all;
280:   MPI_Comm           comm;
281:   PetscMPIInt        rank, size;

283:   PetscObjectGetComm((PetscObject)snes, &comm);
284:   MPI_Comm_rank(comm, &rank);
285:   MPI_Comm_rank(comm, &size);

287:   MatGetOwnershipRanges(Jpre, &Jranges);
288:   MatGetOwnershipRange(Jpre, &Jrstart, NULL);
289:   MatGetOwnershipRangesColumn(tao->hessian, &rranges);
290:   MatGetOwnershipRangesColumn(tao->hessian, &cranges);

292:   VecGetArrayRead(X, &Xarr);

294:   /* (1) insert Z and Ci to the 4th block of Jpre -- overwrite existing values */
295:   if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
296:     vals[0] = 1.0;
297:     for (i = 0; i < pdipm->nci; i++) {
298:       row     = Jrstart + pdipm->off_z + i;
299:       cols[0] = Jrstart + pdipm->off_lambdai + i;
300:       cols[1] = row;
301:       vals[1] = Xarr[pdipm->off_lambdai + i] / Xarr[pdipm->off_z + i];
302:       MatSetValues(Jpre, 1, &row, 2, cols, vals, INSERT_VALUES);
303:     }
304:   } else {
305:     for (i = 0; i < pdipm->nci; i++) {
306:       row     = Jrstart + pdipm->off_z + i;
307:       cols[0] = Jrstart + pdipm->off_lambdai + i;
308:       cols[1] = row;
309:       vals[0] = Xarr[pdipm->off_z + i];
310:       vals[1] = Xarr[pdipm->off_lambdai + i];
311:       MatSetValues(Jpre, 1, &row, 2, cols, vals, INSERT_VALUES);
312:     }
313:   }

315:   /* (2) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
316:   if (pdipm->Ng) {
317:     MatGetOwnershipRange(tao->jacobian_equality, &rjstart, NULL);
318:     for (i = 0; i < pdipm->ng; i++) {
319:       row = Jrstart + pdipm->off_lambdae + i;

321:       MatGetRow(tao->jacobian_equality, i + rjstart, &nc, &aj, &aa);
322:       proc = 0;
323:       for (j = 0; j < nc; j++) {
324:         while (aj[j] >= cranges[proc + 1]) proc++;
325:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
326:         MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
327:       }
328:       MatRestoreRow(tao->jacobian_equality, i + rjstart, &nc, &aj, &aa);
329:       if (pdipm->kkt_pd) {
330:         /* add shift \delta_c */
331:         MatSetValue(Jpre, row, row, -pdipm->deltac, INSERT_VALUES);
332:       }
333:     }
334:   }

336:   /* (3) insert 3rd row block of Jpre: [ -grad h, 0, deltac, I] */
337:   if (pdipm->Nh) {
338:     MatGetOwnershipRange(tao->jacobian_inequality, &rjstart, NULL);
339:     for (i = 0; i < pdipm->nh; i++) {
340:       row = Jrstart + pdipm->off_lambdai + i;
341:       MatGetRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, &aa);
342:       proc = 0;
343:       for (j = 0; j < nc; j++) {
344:         while (aj[j] >= cranges[proc + 1]) proc++;
345:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
346:         MatSetValue(Jpre, row, cols[0], -aa[j], INSERT_VALUES);
347:       }
348:       MatRestoreRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, &aa);
349:       if (pdipm->kkt_pd) {
350:         /* add shift \delta_c */
351:         MatSetValue(Jpre, row, row, -pdipm->deltac, INSERT_VALUES);
352:       }
353:     }
354:   }

356:   /* (4) insert 1st row block of Jpre: [Wxx, grad g', -grad h', 0] */
357:   if (pdipm->Ng) { /* grad g' */
358:     MatTranspose(tao->jacobian_equality, MAT_REUSE_MATRIX, &pdipm->jac_equality_trans);
359:   }
360:   if (pdipm->Nh) { /* grad h' */
361:     MatTranspose(tao->jacobian_inequality, MAT_REUSE_MATRIX, &pdipm->jac_inequality_trans);
362:   }

364:   VecPlaceArray(pdipm->x, Xarr);
365:   TaoComputeHessian(tao, pdipm->x, tao->hessian, tao->hessian_pre);
366:   VecResetArray(pdipm->x);

368:   MatGetOwnershipRange(tao->hessian, &rjstart, NULL);
369:   for (i = 0; i < pdipm->nx; i++) {
370:     row = Jrstart + i;

372:     /* insert Wxx = fxx + ... -- provided by user */
373:     MatGetRow(tao->hessian, i + rjstart, &nc, &aj, &aa);
374:     proc = 0;
375:     for (j = 0; j < nc; j++) {
376:       while (aj[j] >= cranges[proc + 1]) proc++;
377:       cols[0] = aj[j] - cranges[proc] + Jranges[proc];
378:       if (row == cols[0] && pdipm->kkt_pd) {
379:         /* add shift deltaw to Wxx component */
380:         MatSetValue(Jpre, row, cols[0], aa[j] + pdipm->deltaw, INSERT_VALUES);
381:       } else {
382:         MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
383:       }
384:     }
385:     MatRestoreRow(tao->hessian, i + rjstart, &nc, &aj, &aa);

387:     /* insert grad g' */
388:     if (pdipm->ng) {
389:       MatGetRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, &aa);
390:       MatGetOwnershipRanges(tao->jacobian_equality, &ranges);
391:       proc = 0;
392:       for (j = 0; j < nc; j++) {
393:         /* find row ownership of */
394:         while (aj[j] >= ranges[proc + 1]) proc++;
395:         nx_all  = rranges[proc + 1] - rranges[proc];
396:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
397:         MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
398:       }
399:       MatRestoreRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, &aa);
400:     }

402:     /* insert -grad h' */
403:     if (pdipm->nh) {
404:       MatGetRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, &aa);
405:       MatGetOwnershipRanges(tao->jacobian_inequality, &ranges);
406:       proc = 0;
407:       for (j = 0; j < nc; j++) {
408:         /* find row ownership of */
409:         while (aj[j] >= ranges[proc + 1]) proc++;
410:         nx_all  = rranges[proc + 1] - rranges[proc];
411:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
412:         MatSetValue(Jpre, row, cols[0], -aa[j], INSERT_VALUES);
413:       }
414:       MatRestoreRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, &aa);
415:     }
416:   }
417:   VecRestoreArrayRead(X, &Xarr);

419:   /* (6) assemble Jpre and J */
420:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
421:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);

423:   if (J != Jpre) {
424:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
425:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
426:   }
427:   return 0;
428: }

430: /*
431:    TaoSnesFunction_PDIPM - Evaluate KKT function at X

433:    Input Parameter:
434:    snes - SNES context
435:    X - KKT Vector
436:    *ctx - pdipm

438:    Output Parameter:
439:    F - Updated Lagrangian vector
440: */
441: static PetscErrorCode TaoSNESFunction_PDIPM(SNES snes, Vec X, Vec F, void *ctx)
442: {
443:   Tao                tao   = (Tao)ctx;
444:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
445:   PetscScalar       *Farr;
446:   Vec                x, L1;
447:   PetscInt           i;
448:   const PetscScalar *Xarr, *carr, *zarr, *larr;

450:   VecSet(F, 0.0);

452:   VecGetArrayRead(X, &Xarr);
453:   VecGetArrayWrite(F, &Farr);

455:   /* (0) Evaluate f, fx, gradG, gradH at X.x Note: pdipm->x is not changed below */
456:   x = pdipm->x;
457:   VecPlaceArray(x, Xarr);
458:   TaoPDIPMEvaluateFunctionsAndJacobians(tao, x);

460:   /* Update ce, ci, and Jci at X.x */
461:   TaoPDIPMUpdateConstraints(tao, x);
462:   VecResetArray(x);

464:   /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
465:   L1 = pdipm->x;
466:   VecPlaceArray(L1, Farr); /* L1 = 0.0 */
467:   if (pdipm->Nci) {
468:     if (pdipm->Nh) {
469:       /* L1 += gradH'*DI. Note: tao->DI is not changed below */
470:       VecPlaceArray(tao->DI, Xarr + pdipm->off_lambdai);
471:       MatMultTransposeAdd(tao->jacobian_inequality, tao->DI, L1, L1);
472:       VecResetArray(tao->DI);
473:     }

475:     /* L1 += Jci_xb'*lambdai_xb */
476:     VecPlaceArray(pdipm->lambdai_xb, Xarr + pdipm->off_lambdai + pdipm->nh);
477:     MatMultTransposeAdd(pdipm->Jci_xb, pdipm->lambdai_xb, L1, L1);
478:     VecResetArray(pdipm->lambdai_xb);

480:     /* L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
481:     VecScale(L1, -1.0);
482:   }

484:   /* L1 += fx */
485:   VecAXPY(L1, 1.0, tao->gradient);

487:   if (pdipm->Nce) {
488:     if (pdipm->Ng) {
489:       /* L1 += gradG'*DE. Note: tao->DE is not changed below */
490:       VecPlaceArray(tao->DE, Xarr + pdipm->off_lambdae);
491:       MatMultTransposeAdd(tao->jacobian_equality, tao->DE, L1, L1);
492:       VecResetArray(tao->DE);
493:     }
494:     if (pdipm->Nxfixed) {
495:       /* L1 += Jce_xfixed'*lambdae_xfixed */
496:       VecPlaceArray(pdipm->lambdae_xfixed, Xarr + pdipm->off_lambdae + pdipm->ng);
497:       MatMultTransposeAdd(pdipm->Jce_xfixed, pdipm->lambdae_xfixed, L1, L1);
498:       VecResetArray(pdipm->lambdae_xfixed);
499:     }
500:   }
501:   VecResetArray(L1);

503:   /* (2) L2 = ce(x) */
504:   if (pdipm->Nce) {
505:     VecGetArrayRead(pdipm->ce, &carr);
506:     for (i = 0; i < pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
507:     VecRestoreArrayRead(pdipm->ce, &carr);
508:   }

510:   if (pdipm->Nci) {
511:     if (pdipm->solve_symmetric_kkt) {
512:       /* (3) L3 = z - ci(x);
513:          (4) L4 = Lambdai * e - mu/z *e  */
514:       VecGetArrayRead(pdipm->ci, &carr);
515:       larr = Xarr + pdipm->off_lambdai;
516:       zarr = Xarr + pdipm->off_z;
517:       for (i = 0; i < pdipm->nci; i++) {
518:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
519:         Farr[pdipm->off_z + i]       = larr[i] - pdipm->mu / zarr[i];
520:       }
521:       VecRestoreArrayRead(pdipm->ci, &carr);
522:     } else {
523:       /* (3) L3 = z - ci(x);
524:          (4) L4 = Z * Lambdai * e - mu * e  */
525:       VecGetArrayRead(pdipm->ci, &carr);
526:       larr = Xarr + pdipm->off_lambdai;
527:       zarr = Xarr + pdipm->off_z;
528:       for (i = 0; i < pdipm->nci; i++) {
529:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
530:         Farr[pdipm->off_z + i]       = zarr[i] * larr[i] - pdipm->mu;
531:       }
532:       VecRestoreArrayRead(pdipm->ci, &carr);
533:     }
534:   }

536:   VecRestoreArrayRead(X, &Xarr);
537:   VecRestoreArrayWrite(F, &Farr);
538:   return 0;
539: }

541: /*
542:   Evaluate F(X); then update update tao->gnorm0, tao->step = mu,
543:   tao->residual = norm2(F_x,F_z) and tao->cnorm = norm2(F_ce,F_ci).
544: */
545: static PetscErrorCode TaoSNESFunction_PDIPM_residual(SNES snes, Vec X, Vec F, void *ctx)
546: {
547:   Tao                tao   = (Tao)ctx;
548:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
549:   PetscScalar       *Farr, *tmparr;
550:   Vec                L1;
551:   PetscInt           i;
552:   PetscReal          res[2], cnorm[2];
553:   const PetscScalar *Xarr = NULL;

555:   TaoSNESFunction_PDIPM(snes, X, F, (void *)tao);
556:   VecGetArrayWrite(F, &Farr);
557:   VecGetArrayRead(X, &Xarr);

559:   /* compute res[0] = norm2(F_x) */
560:   L1 = pdipm->x;
561:   VecPlaceArray(L1, Farr);
562:   VecNorm(L1, NORM_2, &res[0]);
563:   VecResetArray(L1);

565:   /* compute res[1] = norm2(F_z), cnorm[1] = norm2(F_ci) */
566:   if (pdipm->z) {
567:     if (pdipm->solve_symmetric_kkt) {
568:       VecPlaceArray(pdipm->z, Farr + pdipm->off_z);
569:       if (pdipm->Nci) {
570:         VecGetArrayWrite(pdipm->z, &tmparr);
571:         for (i = 0; i < pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
572:         VecRestoreArrayWrite(pdipm->z, &tmparr);
573:       }

575:       VecNorm(pdipm->z, NORM_2, &res[1]);

577:       if (pdipm->Nci) {
578:         VecGetArrayWrite(pdipm->z, &tmparr);
579:         for (i = 0; i < pdipm->nci; i++) tmparr[i] /= Xarr[pdipm->off_z + i];
580:         VecRestoreArrayWrite(pdipm->z, &tmparr);
581:       }
582:       VecResetArray(pdipm->z);
583:     } else { /* !solve_symmetric_kkt */
584:       VecPlaceArray(pdipm->z, Farr + pdipm->off_z);
585:       VecNorm(pdipm->z, NORM_2, &res[1]);
586:       VecResetArray(pdipm->z);
587:     }

589:     VecPlaceArray(pdipm->ci, Farr + pdipm->off_lambdai);
590:     VecNorm(pdipm->ci, NORM_2, &cnorm[1]);
591:     VecResetArray(pdipm->ci);
592:   } else {
593:     res[1]   = 0.0;
594:     cnorm[1] = 0.0;
595:   }

597:   /* compute cnorm[0] = norm2(F_ce) */
598:   if (pdipm->Nce) {
599:     VecPlaceArray(pdipm->ce, Farr + pdipm->off_lambdae);
600:     VecNorm(pdipm->ce, NORM_2, &cnorm[0]);
601:     VecResetArray(pdipm->ce);
602:   } else cnorm[0] = 0.0;

604:   VecRestoreArrayWrite(F, &Farr);
605:   VecRestoreArrayRead(X, &Xarr);

607:   tao->gnorm0   = tao->residual;
608:   tao->residual = PetscSqrtReal(res[0] * res[0] + res[1] * res[1]);
609:   tao->cnorm    = PetscSqrtReal(cnorm[0] * cnorm[0] + cnorm[1] * cnorm[1]);
610:   tao->step     = pdipm->mu;
611:   return 0;
612: }

614: /*
615:   KKTAddShifts - Check the inertia of Cholesky factor of KKT matrix.
616:   If it does not match the numbers of prime and dual variables, add shifts to the KKT matrix.
617: */
618: static PetscErrorCode KKTAddShifts(Tao tao, SNES snes, Vec X)
619: {
620:   TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
621:   KSP        ksp;
622:   PC         pc;
623:   Mat        Factor;
624:   PetscBool  isCHOL;
625:   PetscInt   nneg, nzero, npos;

627:   /* Get the inertia of Cholesky factor */
628:   SNESGetKSP(snes, &ksp);
629:   KSPGetPC(ksp, &pc);
630:   PetscObjectTypeCompare((PetscObject)pc, PCCHOLESKY, &isCHOL);
631:   if (!isCHOL) return 0;

633:   PCFactorGetMatrix(pc, &Factor);
634:   MatGetInertia(Factor, &nneg, &nzero, &npos);

636:   if (npos < pdipm->Nx + pdipm->Nci) {
637:     pdipm->deltaw = PetscMax(pdipm->lastdeltaw / 3, 1.e-4 * PETSC_MACHINE_EPSILON);
638:     PetscInfo(tao, "Test reduced deltaw=%g; previous MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT ", npos %" PetscInt_FMT "(<%" PetscInt_FMT ")\n", (double)pdipm->deltaw, nneg, nzero, npos, pdipm->Nx + pdipm->Nci);
639:     TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
640:     PCSetUp(pc);
641:     MatGetInertia(Factor, &nneg, &nzero, &npos);

643:     if (npos < pdipm->Nx + pdipm->Nci) {
644:       pdipm->deltaw = pdipm->lastdeltaw;                                           /* in case reduction update does not help, this prevents that step from impacting increasing update */
645:       while (npos < pdipm->Nx + pdipm->Nci && pdipm->deltaw <= 1. / PETSC_SMALL) { /* increase deltaw */
646:         PetscInfo(tao, "  deltaw=%g fails, MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT ", npos %" PetscInt_FMT "(<%" PetscInt_FMT ")\n", (double)pdipm->deltaw, nneg, nzero, npos, pdipm->Nx + pdipm->Nci);
647:         pdipm->deltaw = PetscMin(8 * pdipm->deltaw, PetscPowReal(10, 20));
648:         TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
649:         PCSetUp(pc);
650:         MatGetInertia(Factor, &nneg, &nzero, &npos);
651:       }


655:       PetscInfo(tao, "Updated deltaw %g\n", (double)pdipm->deltaw);
656:       pdipm->lastdeltaw = pdipm->deltaw;
657:       pdipm->deltaw     = 0.0;
658:     }
659:   }

661:   if (nzero) { /* Jacobian is singular */
662:     if (pdipm->deltac == 0.0) {
663:       pdipm->deltac = PETSC_SQRT_MACHINE_EPSILON;
664:     } else {
665:       pdipm->deltac = pdipm->deltac * PetscPowReal(pdipm->mu, .25);
666:     }
667:     PetscInfo(tao, "Updated deltac=%g, MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT "(!=0), npos %" PetscInt_FMT "\n", (double)pdipm->deltac, nneg, nzero, npos);
668:     TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
669:     PCSetUp(pc);
670:     MatGetInertia(Factor, &nneg, &nzero, &npos);
671:   }
672:   return 0;
673: }

675: /*
676:   PCPreSolve_PDIPM -- called betwee MatFactorNumeric() and MatSolve()
677: */
678: PetscErrorCode PCPreSolve_PDIPM(PC pc, KSP ksp)
679: {
680:   Tao        tao;
681:   TAO_PDIPM *pdipm;

683:   KSPGetApplicationContext(ksp, &tao);
684:   pdipm = (TAO_PDIPM *)tao->data;
685:   KKTAddShifts(tao, pdipm->snes, pdipm->X);
686:   return 0;
687: }

689: /*
690:    SNESLineSearch_PDIPM - Custom line search used with PDIPM.

692:    Collective

694:    Notes:
695:    This routine employs a simple backtracking line-search to keep
696:    the slack variables (z) and inequality constraints Lagrange multipliers
697:    (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
698:    alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
699:    slack variables are updated as X = X - alpha_d*dx. The constraint multipliers
700:    are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
701:    is also updated as mu = mu + z'lambdai/Nci
702: */
703: static PetscErrorCode SNESLineSearch_PDIPM(SNESLineSearch linesearch, void *ctx)
704: {
705:   Tao                tao   = (Tao)ctx;
706:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
707:   SNES               snes;
708:   Vec                X, F, Y;
709:   PetscInt           i, iter;
710:   PetscReal          alpha_p = 1.0, alpha_d = 1.0, alpha[4];
711:   PetscScalar       *Xarr, *z, *lambdai, dot, *taosolarr;
712:   const PetscScalar *dXarr, *dz, *dlambdai;

714:   SNESLineSearchGetSNES(linesearch, &snes);
715:   SNESGetIterationNumber(snes, &iter);

717:   SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
718:   SNESLineSearchGetVecs(linesearch, &X, &F, &Y, NULL, NULL);

720:   VecGetArrayWrite(X, &Xarr);
721:   VecGetArrayRead(Y, &dXarr);
722:   z  = Xarr + pdipm->off_z;
723:   dz = dXarr + pdipm->off_z;
724:   for (i = 0; i < pdipm->nci; i++) {
725:     if (z[i] - dz[i] < 0.0) alpha_p = PetscMin(alpha_p, 0.9999 * z[i] / dz[i]);
726:   }

728:   lambdai  = Xarr + pdipm->off_lambdai;
729:   dlambdai = dXarr + pdipm->off_lambdai;

731:   for (i = 0; i < pdipm->nci; i++) {
732:     if (lambdai[i] - dlambdai[i] < 0.0) alpha_d = PetscMin(0.9999 * lambdai[i] / dlambdai[i], alpha_d);
733:   }

735:   alpha[0] = alpha_p;
736:   alpha[1] = alpha_d;
737:   VecRestoreArrayRead(Y, &dXarr);
738:   VecRestoreArrayWrite(X, &Xarr);

740:   /* alpha = min(alpha) over all processes */
741:   MPI_Allreduce(alpha, alpha + 2, 2, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)tao));

743:   alpha_p = alpha[2];
744:   alpha_d = alpha[3];

746:   /* X = X - alpha * Y */
747:   VecGetArrayWrite(X, &Xarr);
748:   VecGetArrayRead(Y, &dXarr);
749:   for (i = 0; i < pdipm->nx; i++) Xarr[i] -= alpha_p * dXarr[i];
750:   for (i = 0; i < pdipm->nce; i++) Xarr[i + pdipm->off_lambdae] -= alpha_d * dXarr[i + pdipm->off_lambdae];

752:   for (i = 0; i < pdipm->nci; i++) {
753:     Xarr[i + pdipm->off_lambdai] -= alpha_d * dXarr[i + pdipm->off_lambdai];
754:     Xarr[i + pdipm->off_z] -= alpha_p * dXarr[i + pdipm->off_z];
755:   }
756:   VecGetArrayWrite(tao->solution, &taosolarr);
757:   PetscMemcpy(taosolarr, Xarr, pdipm->nx * sizeof(PetscScalar));
758:   VecRestoreArrayWrite(tao->solution, &taosolarr);

760:   VecRestoreArrayWrite(X, &Xarr);
761:   VecRestoreArrayRead(Y, &dXarr);

763:   /* Update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
764:   if (pdipm->z) {
765:     VecDot(pdipm->z, pdipm->lambdai, &dot);
766:   } else dot = 0.0;

768:   /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu)  */
769:   pdipm->mu = pdipm->mu_update_factor * dot / pdipm->Nci;

771:   /* Update F; get tao->residual and tao->cnorm */
772:   TaoSNESFunction_PDIPM_residual(snes, X, F, (void *)tao);

774:   tao->niter++;
775:   TaoLogConvergenceHistory(tao, pdipm->obj, tao->residual, tao->cnorm, tao->niter);
776:   TaoMonitor(tao, tao->niter, pdipm->obj, tao->residual, tao->cnorm, pdipm->mu);

778:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
779:   if (tao->reason) SNESSetConvergedReason(snes, SNES_CONVERGED_FNORM_ABS);
780:   return 0;
781: }

783: /*
784:    TaoSolve_PDIPM

786:    Input Parameter:
787:    tao - TAO context

789:    Output Parameter:
790:    tao - TAO context
791: */
792: PetscErrorCode TaoSolve_PDIPM(Tao tao)
793: {
794:   TAO_PDIPM     *pdipm = (TAO_PDIPM *)tao->data;
795:   SNESLineSearch linesearch; /* SNESLineSearch context */
796:   Vec            dummy;


800:   /* Initialize all variables */
801:   TaoPDIPMInitializeSolution(tao);

803:   /* Set linesearch */
804:   SNESGetLineSearch(pdipm->snes, &linesearch);
805:   SNESLineSearchSetType(linesearch, SNESLINESEARCHSHELL);
806:   SNESLineSearchShellSetUserFunc(linesearch, SNESLineSearch_PDIPM, tao);
807:   SNESLineSearchSetFromOptions(linesearch);

809:   tao->reason = TAO_CONTINUE_ITERATING;

811:   /* -tao_monitor for iteration 0 and check convergence */
812:   VecDuplicate(pdipm->X, &dummy);
813:   TaoSNESFunction_PDIPM_residual(pdipm->snes, pdipm->X, dummy, (void *)tao);

815:   TaoLogConvergenceHistory(tao, pdipm->obj, tao->residual, tao->cnorm, tao->niter);
816:   TaoMonitor(tao, tao->niter, pdipm->obj, tao->residual, tao->cnorm, pdipm->mu);
817:   VecDestroy(&dummy);
818:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
819:   if (tao->reason) SNESSetConvergedReason(pdipm->snes, SNES_CONVERGED_FNORM_ABS);

821:   while (tao->reason == TAO_CONTINUE_ITERATING) {
822:     SNESConvergedReason reason;
823:     SNESSolve(pdipm->snes, NULL, pdipm->X);

825:     /* Check SNES convergence */
826:     SNESGetConvergedReason(pdipm->snes, &reason);
827:     if (reason < 0) PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes), "SNES solve did not converged due to reason %s\n", SNESConvergedReasons[reason]);

829:     /* Check TAO convergence */
831:   }
832:   return 0;
833: }

835: /*
836:   TaoView_PDIPM - View PDIPM

838:    Input Parameter:
839:     tao - TAO object
840:     viewer - PetscViewer

842:    Output:
843: */
844: PetscErrorCode TaoView_PDIPM(Tao tao, PetscViewer viewer)
845: {
846:   TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;

848:   tao->constrained = PETSC_TRUE;
849:   PetscViewerASCIIPushTab(viewer);
850:   PetscViewerASCIIPrintf(viewer, "Number of prime=%" PetscInt_FMT ", Number of dual=%" PetscInt_FMT "\n", pdipm->Nx + pdipm->Nci, pdipm->Nce + pdipm->Nci);
851:   if (pdipm->kkt_pd) PetscViewerASCIIPrintf(viewer, "KKT shifts deltaw=%g, deltac=%g\n", (double)pdipm->deltaw, (double)pdipm->deltac);
852:   PetscViewerASCIIPopTab(viewer);
853:   return 0;
854: }

856: /*
857:    TaoSetup_PDIPM - Sets up tao and pdipm

859:    Input Parameter:
860:    tao - TAO object

862:    Output:   pdipm - initialized object
863: */
864: PetscErrorCode TaoSetup_PDIPM(Tao tao)
865: {
866:   TAO_PDIPM         *pdipm = (TAO_PDIPM *)tao->data;
867:   MPI_Comm           comm;
868:   PetscMPIInt        size;
869:   PetscInt           row, col, Jcrstart, Jcrend, k, tmp, nc, proc, *nh_all, *ng_all;
870:   PetscInt           offset, *xa, *xb, i, j, rstart, rend;
871:   PetscScalar        one = 1.0, neg_one = -1.0;
872:   const PetscInt    *cols, *rranges, *cranges, *aj, *ranges;
873:   const PetscScalar *aa, *Xarr;
874:   Mat                J;
875:   Mat                Jce_xfixed_trans, Jci_xb_trans;
876:   PetscInt          *dnz, *onz, rjstart, nx_all, *nce_all, *Jranges, cols1[2];

878:   PetscObjectGetComm((PetscObject)tao, &comm);
879:   MPI_Comm_size(comm, &size);

881:   /* (1) Setup Bounds and create Tao vectors */
882:   TaoPDIPMSetUpBounds(tao);

884:   if (!tao->gradient) {
885:     VecDuplicate(tao->solution, &tao->gradient);
886:     VecDuplicate(tao->solution, &tao->stepdirection);
887:   }

889:   /* (2) Get sizes */
890:   /* Size of vector x - This is set by TaoSetSolution */
891:   VecGetSize(tao->solution, &pdipm->Nx);
892:   VecGetLocalSize(tao->solution, &pdipm->nx);

894:   /* Size of equality constraints and vectors */
895:   if (tao->constraints_equality) {
896:     VecGetSize(tao->constraints_equality, &pdipm->Ng);
897:     VecGetLocalSize(tao->constraints_equality, &pdipm->ng);
898:   } else {
899:     pdipm->ng = pdipm->Ng = 0;
900:   }

902:   pdipm->nce = pdipm->ng + pdipm->nxfixed;
903:   pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;

905:   /* Size of inequality constraints and vectors */
906:   if (tao->constraints_inequality) {
907:     VecGetSize(tao->constraints_inequality, &pdipm->Nh);
908:     VecGetLocalSize(tao->constraints_inequality, &pdipm->nh);
909:   } else {
910:     pdipm->nh = pdipm->Nh = 0;
911:   }

913:   pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2 * pdipm->nxbox;
914:   pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2 * pdipm->Nxbox;

916:   /* Full size of the KKT system to be solved */
917:   pdipm->n = pdipm->nx + pdipm->nce + 2 * pdipm->nci;
918:   pdipm->N = pdipm->Nx + pdipm->Nce + 2 * pdipm->Nci;

920:   /* (3) Offsets for subvectors */
921:   pdipm->off_lambdae = pdipm->nx;
922:   pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
923:   pdipm->off_z       = pdipm->off_lambdai + pdipm->nci;

925:   /* (4) Create vectors and subvectors */
926:   /* Ce and Ci vectors */
927:   VecCreate(comm, &pdipm->ce);
928:   VecSetSizes(pdipm->ce, pdipm->nce, pdipm->Nce);
929:   VecSetFromOptions(pdipm->ce);

931:   VecCreate(comm, &pdipm->ci);
932:   VecSetSizes(pdipm->ci, pdipm->nci, pdipm->Nci);
933:   VecSetFromOptions(pdipm->ci);

935:   /* X=[x; lambdae; lambdai; z] for the big KKT system */
936:   VecCreate(comm, &pdipm->X);
937:   VecSetSizes(pdipm->X, pdipm->n, pdipm->N);
938:   VecSetFromOptions(pdipm->X);

940:   /* Subvectors; they share local arrays with X */
941:   VecGetArrayRead(pdipm->X, &Xarr);
942:   /* x shares local array with X.x */
943:   if (pdipm->Nx) VecCreateMPIWithArray(comm, 1, pdipm->nx, pdipm->Nx, Xarr, &pdipm->x);

945:   /* lambdae shares local array with X.lambdae */
946:   if (pdipm->Nce) VecCreateMPIWithArray(comm, 1, pdipm->nce, pdipm->Nce, Xarr + pdipm->off_lambdae, &pdipm->lambdae);

948:   /* tao->DE shares local array with X.lambdae_g */
949:   if (pdipm->Ng) {
950:     VecCreateMPIWithArray(comm, 1, pdipm->ng, pdipm->Ng, Xarr + pdipm->off_lambdae, &tao->DE);

952:     VecCreate(comm, &pdipm->lambdae_xfixed);
953:     VecSetSizes(pdipm->lambdae_xfixed, pdipm->nxfixed, PETSC_DECIDE);
954:     VecSetFromOptions(pdipm->lambdae_xfixed);
955:   }

957:   if (pdipm->Nci) {
958:     /* lambdai shares local array with X.lambdai */
959:     VecCreateMPIWithArray(comm, 1, pdipm->nci, pdipm->Nci, Xarr + pdipm->off_lambdai, &pdipm->lambdai);

961:     /* z for slack variables; it shares local array with X.z */
962:     VecCreateMPIWithArray(comm, 1, pdipm->nci, pdipm->Nci, Xarr + pdipm->off_z, &pdipm->z);
963:   }

965:   /* tao->DI which shares local array with X.lambdai_h */
966:   if (pdipm->Nh) VecCreateMPIWithArray(comm, 1, pdipm->nh, pdipm->Nh, Xarr + pdipm->off_lambdai, &tao->DI);
967:   VecCreate(comm, &pdipm->lambdai_xb);
968:   VecSetSizes(pdipm->lambdai_xb, (pdipm->nci - pdipm->nh), PETSC_DECIDE);
969:   VecSetFromOptions(pdipm->lambdai_xb);

971:   VecRestoreArrayRead(pdipm->X, &Xarr);

973:   /* (5) Create Jacobians Jce_xfixed and Jci */
974:   /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
975:   if (pdipm->Nxfixed) {
976:     /* Create Jce_xfixed */
977:     MatCreate(comm, &pdipm->Jce_xfixed);
978:     MatSetSizes(pdipm->Jce_xfixed, pdipm->nxfixed, pdipm->nx, PETSC_DECIDE, pdipm->Nx);
979:     MatSetFromOptions(pdipm->Jce_xfixed);
980:     MatSeqAIJSetPreallocation(pdipm->Jce_xfixed, 1, NULL);
981:     MatMPIAIJSetPreallocation(pdipm->Jce_xfixed, 1, NULL, 1, NULL);

983:     MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, &Jcrend);
984:     ISGetIndices(pdipm->isxfixed, &cols);
985:     k = 0;
986:     for (row = Jcrstart; row < Jcrend; row++) {
987:       MatSetValues(pdipm->Jce_xfixed, 1, &row, 1, cols + k, &one, INSERT_VALUES);
988:       k++;
989:     }
990:     ISRestoreIndices(pdipm->isxfixed, &cols);
991:     MatAssemblyBegin(pdipm->Jce_xfixed, MAT_FINAL_ASSEMBLY);
992:     MatAssemblyEnd(pdipm->Jce_xfixed, MAT_FINAL_ASSEMBLY);
993:   }

995:   /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
996:   MatCreate(comm, &pdipm->Jci_xb);
997:   MatSetSizes(pdipm->Jci_xb, pdipm->nci - pdipm->nh, pdipm->nx, PETSC_DECIDE, pdipm->Nx);
998:   MatSetFromOptions(pdipm->Jci_xb);
999:   MatSeqAIJSetPreallocation(pdipm->Jci_xb, 1, NULL);
1000:   MatMPIAIJSetPreallocation(pdipm->Jci_xb, 1, NULL, 1, NULL);

1002:   MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, &Jcrend);
1003:   offset = Jcrstart;
1004:   if (pdipm->Nxub) {
1005:     /* Add xub to Jci_xb */
1006:     ISGetIndices(pdipm->isxub, &cols);
1007:     k = 0;
1008:     for (row = offset; row < offset + pdipm->nxub; row++) {
1009:       MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &neg_one, INSERT_VALUES);
1010:       k++;
1011:     }
1012:     ISRestoreIndices(pdipm->isxub, &cols);
1013:   }

1015:   if (pdipm->Nxlb) {
1016:     /* Add xlb to Jci_xb */
1017:     ISGetIndices(pdipm->isxlb, &cols);
1018:     k = 0;
1019:     offset += pdipm->nxub;
1020:     for (row = offset; row < offset + pdipm->nxlb; row++) {
1021:       MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &one, INSERT_VALUES);
1022:       k++;
1023:     }
1024:     ISRestoreIndices(pdipm->isxlb, &cols);
1025:   }

1027:   /* Add xbox to Jci_xb */
1028:   if (pdipm->Nxbox) {
1029:     ISGetIndices(pdipm->isxbox, &cols);
1030:     k = 0;
1031:     offset += pdipm->nxlb;
1032:     for (row = offset; row < offset + pdipm->nxbox; row++) {
1033:       MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &neg_one, INSERT_VALUES);
1034:       tmp = row + pdipm->nxbox;
1035:       MatSetValues(pdipm->Jci_xb, 1, &tmp, 1, cols + k, &one, INSERT_VALUES);
1036:       k++;
1037:     }
1038:     ISRestoreIndices(pdipm->isxbox, &cols);
1039:   }

1041:   MatAssemblyBegin(pdipm->Jci_xb, MAT_FINAL_ASSEMBLY);
1042:   MatAssemblyEnd(pdipm->Jci_xb, MAT_FINAL_ASSEMBLY);
1043:   /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */

1045:   /* (6) Set up ISs for PC Fieldsplit */
1046:   if (pdipm->solve_reduced_kkt) {
1047:     PetscMalloc2(pdipm->nx + pdipm->nce, &xa, 2 * pdipm->nci, &xb);
1048:     for (i = 0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1049:     for (i = 0; i < 2 * pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;

1051:     ISCreateGeneral(comm, pdipm->nx + pdipm->nce, xa, PETSC_OWN_POINTER, &pdipm->is1);
1052:     ISCreateGeneral(comm, 2 * pdipm->nci, xb, PETSC_OWN_POINTER, &pdipm->is2);
1053:   }

1055:   /* (7) Gather offsets from all processes */
1056:   PetscMalloc1(size, &pdipm->nce_all);

1058:   /* Get rstart of KKT matrix */
1059:   MPI_Scan(&pdipm->n, &rstart, 1, MPIU_INT, MPI_SUM, comm);
1060:   rstart -= pdipm->n;

1062:   MPI_Allgather(&pdipm->nce, 1, MPIU_INT, pdipm->nce_all, 1, MPIU_INT, comm);

1064:   PetscMalloc3(size, &ng_all, size, &nh_all, size, &Jranges);
1065:   MPI_Allgather(&rstart, 1, MPIU_INT, Jranges, 1, MPIU_INT, comm);
1066:   MPI_Allgather(&pdipm->nh, 1, MPIU_INT, nh_all, 1, MPIU_INT, comm);
1067:   MPI_Allgather(&pdipm->ng, 1, MPIU_INT, ng_all, 1, MPIU_INT, comm);

1069:   MatGetOwnershipRanges(tao->hessian, &rranges);
1070:   MatGetOwnershipRangesColumn(tao->hessian, &cranges);

1072:   if (pdipm->Ng) {
1073:     TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre);
1074:     MatTranspose(tao->jacobian_equality, MAT_INITIAL_MATRIX, &pdipm->jac_equality_trans);
1075:   }
1076:   if (pdipm->Nh) {
1077:     TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre);
1078:     MatTranspose(tao->jacobian_inequality, MAT_INITIAL_MATRIX, &pdipm->jac_inequality_trans);
1079:   }

1081:   /* Count dnz,onz for preallocation of KKT matrix */
1082:   nce_all = pdipm->nce_all;

1084:   if (pdipm->Nxfixed) MatTranspose(pdipm->Jce_xfixed, MAT_INITIAL_MATRIX, &Jce_xfixed_trans);
1085:   MatTranspose(pdipm->Jci_xb, MAT_INITIAL_MATRIX, &Jci_xb_trans);

1087:   MatPreallocateBegin(comm, pdipm->n, pdipm->n, dnz, onz);

1089:   /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
1090:   TaoPDIPMEvaluateFunctionsAndJacobians(tao, pdipm->x);
1091:   TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);

1093:   /* Insert tao->hessian */
1094:   MatGetOwnershipRange(tao->hessian, &rjstart, NULL);
1095:   for (i = 0; i < pdipm->nx; i++) {
1096:     row = rstart + i;

1098:     MatGetRow(tao->hessian, i + rjstart, &nc, &aj, NULL);
1099:     proc = 0;
1100:     for (j = 0; j < nc; j++) {
1101:       while (aj[j] >= cranges[proc + 1]) proc++;
1102:       col = aj[j] - cranges[proc] + Jranges[proc];
1103:       MatPreallocateSet(row, 1, &col, dnz, onz);
1104:     }
1105:     MatRestoreRow(tao->hessian, i + rjstart, &nc, &aj, NULL);

1107:     if (pdipm->ng) {
1108:       /* Insert grad g' */
1109:       MatGetRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, NULL);
1110:       MatGetOwnershipRanges(tao->jacobian_equality, &ranges);
1111:       proc = 0;
1112:       for (j = 0; j < nc; j++) {
1113:         /* find row ownership of */
1114:         while (aj[j] >= ranges[proc + 1]) proc++;
1115:         nx_all = rranges[proc + 1] - rranges[proc];
1116:         col    = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1117:         MatPreallocateSet(row, 1, &col, dnz, onz);
1118:       }
1119:       MatRestoreRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, NULL);
1120:     }

1122:     /* Insert Jce_xfixed^T' */
1123:     if (pdipm->nxfixed) {
1124:       MatGetRow(Jce_xfixed_trans, i + rjstart, &nc, &aj, NULL);
1125:       MatGetOwnershipRanges(pdipm->Jce_xfixed, &ranges);
1126:       proc = 0;
1127:       for (j = 0; j < nc; j++) {
1128:         /* find row ownership of */
1129:         while (aj[j] >= ranges[proc + 1]) proc++;
1130:         nx_all = rranges[proc + 1] - rranges[proc];
1131:         col    = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1132:         MatPreallocateSet(row, 1, &col, dnz, onz);
1133:       }
1134:       MatRestoreRow(Jce_xfixed_trans, i + rjstart, &nc, &aj, NULL);
1135:     }

1137:     if (pdipm->nh) {
1138:       /* Insert -grad h' */
1139:       MatGetRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, NULL);
1140:       MatGetOwnershipRanges(tao->jacobian_inequality, &ranges);
1141:       proc = 0;
1142:       for (j = 0; j < nc; j++) {
1143:         /* find row ownership of */
1144:         while (aj[j] >= ranges[proc + 1]) proc++;
1145:         nx_all = rranges[proc + 1] - rranges[proc];
1146:         col    = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1147:         MatPreallocateSet(row, 1, &col, dnz, onz);
1148:       }
1149:       MatRestoreRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, NULL);
1150:     }

1152:     /* Insert Jci_xb^T' */
1153:     MatGetRow(Jci_xb_trans, i + rjstart, &nc, &aj, NULL);
1154:     MatGetOwnershipRanges(pdipm->Jci_xb, &ranges);
1155:     proc = 0;
1156:     for (j = 0; j < nc; j++) {
1157:       /* find row ownership of */
1158:       while (aj[j] >= ranges[proc + 1]) proc++;
1159:       nx_all = rranges[proc + 1] - rranges[proc];
1160:       col    = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1161:       MatPreallocateSet(row, 1, &col, dnz, onz);
1162:     }
1163:     MatRestoreRow(Jci_xb_trans, i + rjstart, &nc, &aj, NULL);
1164:   }

1166:   /* 2nd Row block of KKT matrix: [grad Ce, deltac*I, 0, 0] */
1167:   if (pdipm->Ng) {
1168:     MatGetOwnershipRange(tao->jacobian_equality, &rjstart, NULL);
1169:     for (i = 0; i < pdipm->ng; i++) {
1170:       row = rstart + pdipm->off_lambdae + i;

1172:       MatGetRow(tao->jacobian_equality, i + rjstart, &nc, &aj, NULL);
1173:       proc = 0;
1174:       for (j = 0; j < nc; j++) {
1175:         while (aj[j] >= cranges[proc + 1]) proc++;
1176:         col = aj[j] - cranges[proc] + Jranges[proc];
1177:         MatPreallocateSet(row, 1, &col, dnz, onz); /* grad g */
1178:       }
1179:       MatRestoreRow(tao->jacobian_equality, i + rjstart, &nc, &aj, NULL);
1180:     }
1181:   }
1182:   /* Jce_xfixed */
1183:   if (pdipm->Nxfixed) {
1184:     MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, NULL);
1185:     for (i = 0; i < (pdipm->nce - pdipm->ng); i++) {
1186:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1188:       MatGetRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, NULL);

1191:       proc = 0;
1192:       j    = 0;
1193:       while (cols[j] >= cranges[proc + 1]) proc++;
1194:       col = cols[j] - cranges[proc] + Jranges[proc];
1195:       MatPreallocateSet(row, 1, &col, dnz, onz);
1196:       MatRestoreRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, NULL);
1197:     }
1198:   }

1200:   /* 3rd Row block of KKT matrix: [ gradCi, 0, deltac*I, -I] */
1201:   if (pdipm->Nh) {
1202:     MatGetOwnershipRange(tao->jacobian_inequality, &rjstart, NULL);
1203:     for (i = 0; i < pdipm->nh; i++) {
1204:       row = rstart + pdipm->off_lambdai + i;

1206:       MatGetRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, NULL);
1207:       proc = 0;
1208:       for (j = 0; j < nc; j++) {
1209:         while (aj[j] >= cranges[proc + 1]) proc++;
1210:         col = aj[j] - cranges[proc] + Jranges[proc];
1211:         MatPreallocateSet(row, 1, &col, dnz, onz); /* grad h */
1212:       }
1213:       MatRestoreRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, NULL);
1214:     }
1215:     /* I */
1216:     for (i = 0; i < pdipm->nh; i++) {
1217:       row = rstart + pdipm->off_lambdai + i;
1218:       col = rstart + pdipm->off_z + i;
1219:       MatPreallocateSet(row, 1, &col, dnz, onz);
1220:     }
1221:   }

1223:   /* Jci_xb */
1224:   MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, NULL);
1225:   for (i = 0; i < (pdipm->nci - pdipm->nh); i++) {
1226:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1228:     MatGetRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, NULL);
1230:     proc = 0;
1231:     for (j = 0; j < nc; j++) {
1232:       while (cols[j] >= cranges[proc + 1]) proc++;
1233:       col = cols[j] - cranges[proc] + Jranges[proc];
1234:       MatPreallocateSet(row, 1, &col, dnz, onz);
1235:     }
1236:     MatRestoreRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, NULL);
1237:     /* I */
1238:     col = rstart + pdipm->off_z + pdipm->nh + i;
1239:     MatPreallocateSet(row, 1, &col, dnz, onz);
1240:   }

1242:   /* 4-th Row block of KKT matrix: Z and Ci */
1243:   for (i = 0; i < pdipm->nci; i++) {
1244:     row      = rstart + pdipm->off_z + i;
1245:     cols1[0] = rstart + pdipm->off_lambdai + i;
1246:     cols1[1] = row;
1247:     MatPreallocateSet(row, 2, cols1, dnz, onz);
1248:   }

1250:   /* diagonal entry */
1251:   for (i = 0; i < pdipm->n; i++) dnz[i]++; /* diagonal entry */

1253:   /* Create KKT matrix */
1254:   MatCreate(comm, &J);
1255:   MatSetSizes(J, pdipm->n, pdipm->n, PETSC_DECIDE, PETSC_DECIDE);
1256:   MatSetFromOptions(J);
1257:   MatSeqAIJSetPreallocation(J, 0, dnz);
1258:   MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz);
1259:   MatPreallocateEnd(dnz, onz);
1260:   pdipm->K = J;

1262:   /* (8) Insert constant entries to  K */
1263:   /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1264:   MatGetOwnershipRange(J, &rstart, &rend);
1265:   for (i = rstart; i < rend; i++) MatSetValue(J, i, i, 0.0, INSERT_VALUES);
1266:   /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1267:   if (pdipm->kkt_pd) {
1268:     for (i = 0; i < pdipm->nh; i++) {
1269:       row = rstart + i;
1270:       MatSetValue(J, row, row, pdipm->deltaw, INSERT_VALUES);
1271:     }
1272:   }

1274:   /* Row block of K: [ grad Ce, 0, 0, 0] */
1275:   if (pdipm->Nxfixed) {
1276:     MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, NULL);
1277:     for (i = 0; i < (pdipm->nce - pdipm->ng); i++) {
1278:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1280:       MatGetRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, &aa);
1281:       proc = 0;
1282:       for (j = 0; j < nc; j++) {
1283:         while (cols[j] >= cranges[proc + 1]) proc++;
1284:         col = cols[j] - cranges[proc] + Jranges[proc];
1285:         MatSetValue(J, row, col, aa[j], INSERT_VALUES); /* grad Ce */
1286:         MatSetValue(J, col, row, aa[j], INSERT_VALUES); /* grad Ce' */
1287:       }
1288:       MatRestoreRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, &aa);
1289:     }
1290:   }

1292:   /* Row block of K: [ -grad Ci, 0, 0, I] */
1293:   MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, NULL);
1294:   for (i = 0; i < pdipm->nci - pdipm->nh; i++) {
1295:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1297:     MatGetRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, &aa);
1298:     proc = 0;
1299:     for (j = 0; j < nc; j++) {
1300:       while (cols[j] >= cranges[proc + 1]) proc++;
1301:       col = cols[j] - cranges[proc] + Jranges[proc];
1302:       MatSetValue(J, col, row, -aa[j], INSERT_VALUES);
1303:       MatSetValue(J, row, col, -aa[j], INSERT_VALUES);
1304:     }
1305:     MatRestoreRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, &aa);

1307:     col = rstart + pdipm->off_z + pdipm->nh + i;
1308:     MatSetValue(J, row, col, 1, INSERT_VALUES);
1309:   }

1311:   for (i = 0; i < pdipm->nh; i++) {
1312:     row = rstart + pdipm->off_lambdai + i;
1313:     col = rstart + pdipm->off_z + i;
1314:     MatSetValue(J, row, col, 1, INSERT_VALUES);
1315:   }

1317:   /* Row block of K: [ 0, 0, I, ...] */
1318:   for (i = 0; i < pdipm->nci; i++) {
1319:     row = rstart + pdipm->off_z + i;
1320:     col = rstart + pdipm->off_lambdai + i;
1321:     MatSetValue(J, row, col, 1, INSERT_VALUES);
1322:   }

1324:   if (pdipm->Nxfixed) MatDestroy(&Jce_xfixed_trans);
1325:   MatDestroy(&Jci_xb_trans);
1326:   PetscFree3(ng_all, nh_all, Jranges);

1328:   /* (9) Set up nonlinear solver SNES */
1329:   SNESSetFunction(pdipm->snes, NULL, TaoSNESFunction_PDIPM, (void *)tao);
1330:   SNESSetJacobian(pdipm->snes, J, J, TaoSNESJacobian_PDIPM, (void *)tao);

1332:   if (pdipm->solve_reduced_kkt) {
1333:     PC pc;
1334:     KSPGetPC(tao->ksp, &pc);
1335:     PCSetType(pc, PCFIELDSPLIT);
1336:     PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1337:     PCFieldSplitSetIS(pc, "2", pdipm->is2);
1338:     PCFieldSplitSetIS(pc, "1", pdipm->is1);
1339:   }
1340:   SNESSetFromOptions(pdipm->snes);

1342:   /* (10) Setup PCPreSolve() for pdipm->solve_symmetric_kkt */
1343:   if (pdipm->solve_symmetric_kkt) {
1344:     KSP       ksp;
1345:     PC        pc;
1346:     PetscBool isCHOL;
1347:     SNESGetKSP(pdipm->snes, &ksp);
1348:     KSPGetPC(ksp, &pc);
1349:     PCSetPreSolve(pc, PCPreSolve_PDIPM);

1351:     PetscObjectTypeCompare((PetscObject)pc, PCCHOLESKY, &isCHOL);
1352:     if (isCHOL) {
1353:       Mat       Factor;
1354:       PetscBool isMUMPS;
1355:       PCFactorGetMatrix(pc, &Factor);
1356:       PetscObjectTypeCompare((PetscObject)Factor, "mumps", &isMUMPS);
1357:       if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
1358: #if defined(PETSC_HAVE_MUMPS)
1359:         MatMumpsSetIcntl(Factor, 24, 1); /* detection of null pivot rows */
1360:         if (size > 1) { MatMumpsSetIcntl(Factor, 13, 1); /* parallelism of the root node (enable ScaLAPACK) and its splitting */ }
1361: #else
1362:         SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Requires external package MUMPS");
1363: #endif
1364:       }
1365:     }
1366:   }
1367:   return 0;
1368: }

1370: /*
1371:    TaoDestroy_PDIPM - Destroys the pdipm object

1373:    Input:
1374:    full pdipm

1376:    Output:
1377:    Destroyed pdipm
1378: */
1379: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1380: {
1381:   TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;

1383:   /* Freeing Vectors assocaiated with KKT (X) */
1384:   VecDestroy(&pdipm->x);       /* Solution x */
1385:   VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1386:   VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1387:   VecDestroy(&pdipm->z);       /* Slack variables */
1388:   VecDestroy(&pdipm->X);       /* Big KKT system vector [x; lambdae; lambdai; z] */

1390:   /* work vectors */
1391:   VecDestroy(&pdipm->lambdae_xfixed);
1392:   VecDestroy(&pdipm->lambdai_xb);

1394:   /* Legrangian equality and inequality Vec */
1395:   VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1396:   VecDestroy(&pdipm->ci); /* Vec of inequality constraints */

1398:   /* Matrices */
1399:   MatDestroy(&pdipm->Jce_xfixed);
1400:   MatDestroy(&pdipm->Jci_xb)); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb; J(nxbx)] */
1401:   MatDestroy(&pdipm->K);

1403:   /* Index Sets */
1404:   if (pdipm->Nxub) { ISDestroy(&pdipm->isxub); /* Finite upper bound only -inf < x < ub */ }

1406:   if (pdipm->Nxlb) { ISDestroy(&pdipm->isxlb); /* Finite lower bound only  lb <= x < inf */ }

1408:   if (pdipm->Nxfixed) { ISDestroy(&pdipm->isxfixed); /* Fixed variables         lb =  x = ub */ }

1410:   if (pdipm->Nxbox) { ISDestroy(&pdipm->isxbox); /* Boxed variables         lb <= x <= ub */ }

1412:   if (pdipm->Nxfree) { ISDestroy(&pdipm->isxfree); /* Free variables        -inf <= x <= inf */ }

1414:   if (pdipm->solve_reduced_kkt) {
1415:     ISDestroy(&pdipm->is1);
1416:     ISDestroy(&pdipm->is2);
1417:   }

1419:   /* SNES */
1420:   SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1421:   PetscFree(pdipm->nce_all);
1422:   MatDestroy(&pdipm->jac_equality_trans);
1423:   MatDestroy(&pdipm->jac_inequality_trans);

1425:   /* Destroy pdipm */
1426:   PetscFree(tao->data); /* Holding locations of pdipm */

1428:   /* Destroy Dual */
1429:   VecDestroy(&tao->DE); /* equality dual */
1430:   VecDestroy(&tao->DI); /* dinequality dual */
1431:   return 0;
1432: }

1434: PetscErrorCode TaoSetFromOptions_PDIPM(Tao tao, PetscOptionItems *PetscOptionsObject)
1435: {
1436:   TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;

1438:   PetscOptionsHeadBegin(PetscOptionsObject, "PDIPM method for constrained optimization");
1439:   PetscOptionsReal("-tao_pdipm_push_init_slack", "parameter to push initial slack variables away from bounds", NULL, pdipm->push_init_slack, &pdipm->push_init_slack, NULL);
1440:   PetscOptionsReal("-tao_pdipm_push_init_lambdai", "parameter to push initial (inequality) dual variables away from bounds", NULL, pdipm->push_init_lambdai, &pdipm->push_init_lambdai, NULL);
1441:   PetscOptionsBool("-tao_pdipm_solve_reduced_kkt", "Solve reduced KKT system using Schur-complement", NULL, pdipm->solve_reduced_kkt, &pdipm->solve_reduced_kkt, NULL);
1442:   PetscOptionsReal("-tao_pdipm_mu_update_factor", "Update scalar for barrier parameter (mu) update", NULL, pdipm->mu_update_factor, &pdipm->mu_update_factor, NULL);
1443:   PetscOptionsBool("-tao_pdipm_symmetric_kkt", "Solve non reduced symmetric KKT system", NULL, pdipm->solve_symmetric_kkt, &pdipm->solve_symmetric_kkt, NULL);
1444:   PetscOptionsBool("-tao_pdipm_kkt_shift_pd", "Add shifts to make KKT matrix positive definite", NULL, pdipm->kkt_pd, &pdipm->kkt_pd, NULL);
1445:   PetscOptionsHeadEnd();
1446:   return 0;
1447: }

1449: /*MC
1450:   TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.

1452:   Option Database Keys:
1453: +   -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1454: .   -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1455: .   -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1456: .   -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1457: -   -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite

1459:   Level: beginner
1460: M*/
1461: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1462: {
1463:   TAO_PDIPM *pdipm;

1465:   tao->ops->setup          = TaoSetup_PDIPM;
1466:   tao->ops->solve          = TaoSolve_PDIPM;
1467:   tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1468:   tao->ops->view           = TaoView_PDIPM;
1469:   tao->ops->destroy        = TaoDestroy_PDIPM;

1471:   PetscNew(&pdipm);
1472:   tao->data = (void *)pdipm;

1474:   pdipm->nx = pdipm->Nx = 0;
1475:   pdipm->nxfixed = pdipm->Nxfixed = 0;
1476:   pdipm->nxlb = pdipm->Nxlb = 0;
1477:   pdipm->nxub = pdipm->Nxub = 0;
1478:   pdipm->nxbox = pdipm->Nxbox = 0;
1479:   pdipm->nxfree = pdipm->Nxfree = 0;

1481:   pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1482:   pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1483:   pdipm->n = pdipm->N     = 0;
1484:   pdipm->mu               = 1.0;
1485:   pdipm->mu_update_factor = 0.1;

1487:   pdipm->deltaw     = 0.0;
1488:   pdipm->lastdeltaw = 3 * 1.e-4;
1489:   pdipm->deltac     = 0.0;
1490:   pdipm->kkt_pd     = PETSC_FALSE;

1492:   pdipm->push_init_slack     = 1.0;
1493:   pdipm->push_init_lambdai   = 1.0;
1494:   pdipm->solve_reduced_kkt   = PETSC_FALSE;
1495:   pdipm->solve_symmetric_kkt = PETSC_TRUE;

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

1501:   SNESCreate(((PetscObject)tao)->comm, &pdipm->snes);
1502:   SNESSetOptionsPrefix(pdipm->snes, tao->hdr.prefix);
1503:   SNESGetKSP(pdipm->snes, &tao->ksp);
1504:   PetscObjectReference((PetscObject)tao->ksp);
1505:   KSPSetApplicationContext(tao->ksp, (void *)tao);
1506:   return 0;
1507: }