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 on tao

  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) {
196:     ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);
197:   }
198:   if (pdipm->Nxub) {
199:     ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);
200:   }
201:   if (pdipm->Nxfixed) {
202:     ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);
203:   }
204:   if (pdipm->Nxbox) {
205:     ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);
206:   }
207:   if (pdipm->Nxfree) {
208:     ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);
209:   }
210:   PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);
211:   return 0;
212: }

214: /*
215:    TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
216:    X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
217:      four subvectors need to be initialized and its values copied over to X. Instead
218:      of copying, we use VecPlace/ResetArray functions to share the memory locations for
219:      X and the subvectors

221:    Collective

223:    Input Parameter:
224: .  tao - Tao context

226:    Level: beginner
227: */
228: static PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
229: {
230:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
231:   PetscScalar       *Xarr,*z,*lambdai;
232:   PetscInt          i;
233:   const PetscScalar *xarr,*h;

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

237:   /* Set Initialize X.x = tao->solution */
238:   VecGetArrayRead(tao->solution,&xarr);
239:   PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));
240:   VecRestoreArrayRead(tao->solution,&xarr);

242:   /* Initialize X.lambdae = 0.0 */
243:   if (pdipm->lambdae) {
244:     VecSet(pdipm->lambdae,0.0);
245:   }

247:   /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
248:   if (pdipm->Nci) {
249:     VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
250:     VecSet(pdipm->z,pdipm->push_init_slack);

252:     /* Additional modification for X.lambdai and X.z */
253:     VecGetArrayWrite(pdipm->lambdai,&lambdai);
254:     VecGetArrayWrite(pdipm->z,&z);
255:     if (pdipm->Nh) {
256:       VecGetArrayRead(tao->constraints_inequality,&h);
257:       for (i=0; i < pdipm->nh; i++) {
258:         if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
259:         if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
260:       }
261:       VecRestoreArrayRead(tao->constraints_inequality,&h);
262:     }
263:     VecRestoreArrayWrite(pdipm->lambdai,&lambdai);
264:     VecRestoreArrayWrite(pdipm->z,&z);
265:   }

267:   VecRestoreArrayWrite(pdipm->X,&Xarr);
268:   return 0;
269: }

271: /*
272:    TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X

274:    Input Parameter:
275:    snes - SNES context
276:    X - KKT Vector
277:    *ctx - pdipm context

279:    Output Parameter:
280:    J - Hessian matrix
281:    Jpre - Preconditioner
282: */
283: static PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
284: {
285:   Tao               tao=(Tao)ctx;
286:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
287:   PetscInt          i,row,cols[2],Jrstart,rjstart,nc,j;
288:   const PetscInt    *aj,*ranges,*Jranges,*rranges,*cranges;
289:   const PetscScalar *Xarr,*aa;
290:   PetscScalar       vals[2];
291:   PetscInt          proc,nx_all,*nce_all=pdipm->nce_all;
292:   MPI_Comm          comm;
293:   PetscMPIInt       rank,size;
294:   Mat               jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;

296:   PetscObjectGetComm((PetscObject)snes,&comm);
297:   MPI_Comm_rank(comm,&rank);
298:   MPI_Comm_rank(comm,&size);

300:   MatGetOwnershipRanges(Jpre,&Jranges);
301:   MatGetOwnershipRange(Jpre,&Jrstart,NULL);
302:   MatGetOwnershipRangesColumn(tao->hessian,&rranges);
303:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

305:   VecGetArrayRead(X,&Xarr);

307:   /* (1) insert Z and Ci to the 4th block of Jpre -- overwrite existing values */
308:   if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
309:     vals[0] = 1.0;
310:     for (i=0; i < pdipm->nci; i++) {
311:         row     = Jrstart + pdipm->off_z + i;
312:         cols[0] = Jrstart + pdipm->off_lambdai + i;
313:         cols[1] = row;
314:         vals[1] = Xarr[pdipm->off_lambdai + i]/Xarr[pdipm->off_z + i];
315:         MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
316:     }
317:   } else {
318:     for (i=0; i < pdipm->nci; i++) {
319:       row     = Jrstart + pdipm->off_z + i;
320:       cols[0] = Jrstart + pdipm->off_lambdai + i;
321:       cols[1] = row;
322:       vals[0] = Xarr[pdipm->off_z + i];
323:       vals[1] = Xarr[pdipm->off_lambdai + i];
324:       MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
325:     }
326:   }

328:   /* (2) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
329:   if (pdipm->Ng) {
330:     MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
331:     for (i=0; i<pdipm->ng; i++) {
332:       row = Jrstart + pdipm->off_lambdae + i;

334:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
335:       proc = 0;
336:       for (j=0; j < nc; j++) {
337:         while (aj[j] >= cranges[proc+1]) proc++;
338:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
339:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
340:       }
341:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
342:       if (pdipm->kkt_pd) {
343:         /* add shift \delta_c */
344:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
345:       }
346:     }
347:   }

349:   /* (3) insert 3rd row block of Jpre: [ -grad h, 0, deltac, I] */
350:   if (pdipm->Nh) {
351:     MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
352:     for (i=0; i < pdipm->nh; i++) {
353:       row = Jrstart + pdipm->off_lambdai + i;
354:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
355:       proc = 0;
356:       for (j=0; j < nc; j++) {
357:         while (aj[j] >= cranges[proc+1]) proc++;
358:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
359:         MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
360:       }
361:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
362:       if (pdipm->kkt_pd) {
363:         /* add shift \delta_c */
364:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
365:       }
366:     }
367:   }

369:   /* (4) insert 1st row block of Jpre: [Wxx, grad g', -grad h', 0] */
370:   if (pdipm->Ng) { /* grad g' */
371:     MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
372:   }
373:   if (pdipm->Nh) { /* grad h' */
374:     MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
375:   }

377:   VecPlaceArray(pdipm->x,Xarr);
378:   TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
379:   VecResetArray(pdipm->x);

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

385:     /* insert Wxx = fxx + ... -- provided by user */
386:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
387:     proc = 0;
388:     for (j=0; j < nc; j++) {
389:       while (aj[j] >= cranges[proc+1]) proc++;
390:       cols[0] = aj[j] - cranges[proc] + Jranges[proc];
391:       if (row == cols[0] && pdipm->kkt_pd) {
392:         /* add shift deltaw to Wxx component */
393:         MatSetValue(Jpre,row,cols[0],aa[j]+pdipm->deltaw,INSERT_VALUES);
394:       } else {
395:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
396:       }
397:     }
398:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);

400:     /* insert grad g' */
401:     if (pdipm->ng) {
402:       MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
403:       MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
404:       proc = 0;
405:       for (j=0; j < nc; j++) {
406:         /* find row ownership of */
407:         while (aj[j] >= ranges[proc+1]) proc++;
408:         nx_all = rranges[proc+1] - rranges[proc];
409:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
410:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
411:       }
412:       MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
413:     }

415:     /* insert -grad h' */
416:     if (pdipm->nh) {
417:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
418:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
419:       proc = 0;
420:       for (j=0; j < nc; j++) {
421:         /* find row ownership of */
422:         while (aj[j] >= ranges[proc+1]) proc++;
423:         nx_all = rranges[proc+1] - rranges[proc];
424:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
425:         MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
426:       }
427:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
428:     }
429:   }
430:   VecRestoreArrayRead(X,&Xarr);

432:   /* (6) assemble Jpre and J */
433:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
434:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);

436:   if (J != Jpre) {
437:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
438:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
439:   }
440:   return 0;
441: }

443: /*
444:    TaoSnesFunction_PDIPM - Evaluate KKT function at X

446:    Input Parameter:
447:    snes - SNES context
448:    X - KKT Vector
449:    *ctx - pdipm

451:    Output Parameter:
452:    F - Updated Lagrangian vector
453: */
454: static PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
455: {
456:   Tao               tao=(Tao)ctx;
457:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
458:   PetscScalar       *Farr;
459:   Vec               x,L1;
460:   PetscInt          i;
461:   const PetscScalar *Xarr,*carr,*zarr,*larr;

463:   VecSet(F,0.0);

465:   VecGetArrayRead(X,&Xarr);
466:   VecGetArrayWrite(F,&Farr);

468:   /* (0) Evaluate f, fx, gradG, gradH at X.x Note: pdipm->x is not changed below */
469:   x = pdipm->x;
470:   VecPlaceArray(x,Xarr);
471:   TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);

473:   /* Update ce, ci, and Jci at X.x */
474:   TaoPDIPMUpdateConstraints(tao,x);
475:   VecResetArray(x);

477:   /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
478:   L1 = pdipm->x;
479:   VecPlaceArray(L1,Farr); /* L1 = 0.0 */
480:   if (pdipm->Nci) {
481:     if (pdipm->Nh) {
482:       /* L1 += gradH'*DI. Note: tao->DI is not changed below */
483:       VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
484:       MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
485:       VecResetArray(tao->DI);
486:     }

488:     /* L1 += Jci_xb'*lambdai_xb */
489:     VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
490:     MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
491:     VecResetArray(pdipm->lambdai_xb);

493:     /* L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
494:     VecScale(L1,-1.0);
495:   }

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

500:   if (pdipm->Nce) {
501:     if (pdipm->Ng) {
502:       /* L1 += gradG'*DE. Note: tao->DE is not changed below */
503:       VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
504:       MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
505:       VecResetArray(tao->DE);
506:     }
507:     if (pdipm->Nxfixed) {
508:       /* L1 += Jce_xfixed'*lambdae_xfixed */
509:       VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
510:       MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
511:       VecResetArray(pdipm->lambdae_xfixed);
512:     }
513:   }
514:   VecResetArray(L1);

516:   /* (2) L2 = ce(x) */
517:   if (pdipm->Nce) {
518:     VecGetArrayRead(pdipm->ce,&carr);
519:     for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
520:     VecRestoreArrayRead(pdipm->ce,&carr);
521:   }

523:   if (pdipm->Nci) {
524:     if (pdipm->solve_symmetric_kkt) {
525:       /* (3) L3 = z - ci(x);
526:          (4) L4 = Lambdai * e - mu/z *e  */
527:       VecGetArrayRead(pdipm->ci,&carr);
528:       larr = Xarr+pdipm->off_lambdai;
529:       zarr = Xarr+pdipm->off_z;
530:       for (i=0; i<pdipm->nci; i++) {
531:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
532:         Farr[pdipm->off_z       + i] = larr[i] - pdipm->mu/zarr[i];
533:       }
534:       VecRestoreArrayRead(pdipm->ci,&carr);
535:     } else {
536:       /* (3) L3 = z - ci(x);
537:          (4) L4 = Z * Lambdai * e - mu * e  */
538:       VecGetArrayRead(pdipm->ci,&carr);
539:       larr = Xarr+pdipm->off_lambdai;
540:       zarr = Xarr+pdipm->off_z;
541:       for (i=0; i<pdipm->nci; i++) {
542:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
543:         Farr[pdipm->off_z       + i] = zarr[i]*larr[i] - pdipm->mu;
544:       }
545:       VecRestoreArrayRead(pdipm->ci,&carr);
546:     }
547:   }

549:   VecRestoreArrayRead(X,&Xarr);
550:   VecRestoreArrayWrite(F,&Farr);
551:   return 0;
552: }

554: /*
555:   Evaluate F(X); then update update tao->gnorm0, tao->step = mu,
556:   tao->residual = norm2(F_x,F_z) and tao->cnorm = norm2(F_ce,F_ci).
557: */
558: static PetscErrorCode TaoSNESFunction_PDIPM_residual(SNES snes,Vec X,Vec F,void *ctx)
559: {
560:   Tao               tao=(Tao)ctx;
561:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
562:   PetscScalar       *Farr,*tmparr;
563:   Vec               L1;
564:   PetscInt          i;
565:   PetscReal         res[2],cnorm[2];
566:   const PetscScalar *Xarr=NULL;

568:   TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);
569:   VecGetArrayWrite(F,&Farr);
570:   VecGetArrayRead(X,&Xarr);

572:   /* compute res[0] = norm2(F_x) */
573:   L1 = pdipm->x;
574:   VecPlaceArray(L1,Farr);
575:   VecNorm(L1,NORM_2,&res[0]);
576:   VecResetArray(L1);

578:   /* compute res[1] = norm2(F_z), cnorm[1] = norm2(F_ci) */
579:   if (pdipm->z) {
580:     if (pdipm->solve_symmetric_kkt) {
581:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
582:       if (pdipm->Nci) {
583:         VecGetArrayWrite(pdipm->z,&tmparr);
584:         for (i=0; i<pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
585:         VecRestoreArrayWrite(pdipm->z,&tmparr);
586:       }

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

590:       if (pdipm->Nci) {
591:         VecGetArrayWrite(pdipm->z,&tmparr);
592:         for (i=0; i<pdipm->nci; i++) {
593:           tmparr[i] /= Xarr[pdipm->off_z + i];
594:         }
595:         VecRestoreArrayWrite(pdipm->z,&tmparr);
596:       }
597:       VecResetArray(pdipm->z);
598:     } else { /* !solve_symmetric_kkt */
599:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
600:       VecNorm(pdipm->z,NORM_2,&res[1]);
601:       VecResetArray(pdipm->z);
602:     }

604:     VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
605:     VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
606:     VecResetArray(pdipm->ci);
607:   } else {
608:     res[1] = 0.0; cnorm[1] = 0.0;
609:   }

611:   /* compute cnorm[0] = norm2(F_ce) */
612:   if (pdipm->Nce) {
613:     VecPlaceArray(pdipm->ce,Farr+pdipm->off_lambdae);
614:     VecNorm(pdipm->ce,NORM_2,&cnorm[0]);
615:     VecResetArray(pdipm->ce);
616:   } else cnorm[0] = 0.0;

618:   VecRestoreArrayWrite(F,&Farr);
619:   VecRestoreArrayRead(X,&Xarr);

621:   tao->gnorm0   = tao->residual;
622:   tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);
623:   tao->cnorm    = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);
624:   tao->step     = pdipm->mu;
625:   return 0;
626: }

628: /*
629:   KKTAddShifts - Check the inertia of Cholesky factor of KKT matrix.
630:   If it does not match the numbers of prime and dual variables, add shifts to the KKT matrix.
631: */
632: static PetscErrorCode KKTAddShifts(Tao tao,SNES snes,Vec X)
633: {
634:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;
635:   KSP            ksp;
636:   PC             pc;
637:   Mat            Factor;
638:   PetscBool      isCHOL;
639:   PetscInt       nneg,nzero,npos;

641:   /* Get the inertia of Cholesky factor */
642:   SNESGetKSP(snes,&ksp);
643:   KSPGetPC(ksp,&pc);
644:   PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);
645:   if (!isCHOL) return 0;

647:   PCFactorGetMatrix(pc,&Factor);
648:   MatGetInertia(Factor,&nneg,&nzero,&npos);

650:   if (npos < pdipm->Nx+pdipm->Nci) {
651:     pdipm->deltaw = PetscMax(pdipm->lastdeltaw/3, 1.e-4*PETSC_MACHINE_EPSILON);
652:     PetscInfo(tao,"Test reduced deltaw=%g; previous MatInertia: nneg %D, nzero %D, npos %D(<%D)\n",(double)pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
653:     TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
654:     PCSetUp(pc);
655:     MatGetInertia(Factor,&nneg,&nzero,&npos);

657:     if (npos < pdipm->Nx+pdipm->Nci) {
658:       pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
659:       while (npos < pdipm->Nx+pdipm->Nci && pdipm->deltaw <= 1./PETSC_SMALL) { /* increase deltaw */
660:         PetscInfo(tao,"  deltaw=%g fails, MatInertia: nneg %D, nzero %D, npos %D(<%D)\n",(double)pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
661:         pdipm->deltaw = PetscMin(8*pdipm->deltaw,PetscPowReal(10,20));
662:         TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
663:         PCSetUp(pc);
664:         MatGetInertia(Factor,&nneg,&nzero,&npos);
665:       }


669:       PetscInfo(tao,"Updated deltaw %g\n",(double)pdipm->deltaw);
670:       pdipm->lastdeltaw = pdipm->deltaw;
671:       pdipm->deltaw     = 0.0;
672:     }
673:   }

675:   if (nzero) { /* Jacobian is singular */
676:     if (pdipm->deltac == 0.0) {
677:       pdipm->deltac = PETSC_SQRT_MACHINE_EPSILON;
678:     } else {
679:       pdipm->deltac = pdipm->deltac*PetscPowReal(pdipm->mu,.25);
680:     }
681:     PetscInfo(tao,"Updated deltac=%g, MatInertia: nneg %D, nzero %D(!=0), npos %D\n",(double)pdipm->deltac,nneg,nzero,npos);
682:     TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
683:     PCSetUp(pc);
684:     MatGetInertia(Factor,&nneg,&nzero,&npos);
685:   }
686:   return 0;
687: }

689: /*
690:   PCPreSolve_PDIPM -- called betwee MatFactorNumeric() and MatSolve()
691: */
692: PetscErrorCode PCPreSolve_PDIPM(PC pc,KSP ksp)
693: {
694:   Tao            tao;
695:   TAO_PDIPM      *pdipm;

697:   KSPGetApplicationContext(ksp,&tao);
698:   pdipm = (TAO_PDIPM*)tao->data;
699:   KKTAddShifts(tao,pdipm->snes,pdipm->X);
700:   return 0;
701: }

703: /*
704:    SNESLineSearch_PDIPM - Custom line search used with PDIPM.

706:    Collective on TAO

708:    Notes:
709:    This routine employs a simple backtracking line-search to keep
710:    the slack variables (z) and inequality constraints Lagrange multipliers
711:    (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
712:    alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
713:    slack variables are updated as X = X - alpha_d*dx. The constraint multipliers
714:    are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
715:    is also updated as mu = mu + z'lambdai/Nci
716: */
717: static PetscErrorCode SNESLineSearch_PDIPM(SNESLineSearch linesearch,void *ctx)
718: {
719:   Tao               tao=(Tao)ctx;
720:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
721:   SNES              snes;
722:   Vec               X,F,Y;
723:   PetscInt          i,iter;
724:   PetscReal         alpha_p=1.0,alpha_d=1.0,alpha[4];
725:   PetscScalar       *Xarr,*z,*lambdai,dot,*taosolarr;
726:   const PetscScalar *dXarr,*dz,*dlambdai;

728:   SNESLineSearchGetSNES(linesearch,&snes);
729:   SNESGetIterationNumber(snes,&iter);

731:   SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
732:   SNESLineSearchGetVecs(linesearch,&X,&F,&Y,NULL,NULL);

734:   VecGetArrayWrite(X,&Xarr);
735:   VecGetArrayRead(Y,&dXarr);
736:   z  = Xarr + pdipm->off_z;
737:   dz = dXarr + pdipm->off_z;
738:   for (i=0; i < pdipm->nci; i++) {
739:     if (z[i] - dz[i] < 0.0) alpha_p = PetscMin(alpha_p, 0.9999*z[i]/dz[i]);
740:   }

742:   lambdai  = Xarr + pdipm->off_lambdai;
743:   dlambdai = dXarr + pdipm->off_lambdai;

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

749:   alpha[0] = alpha_p;
750:   alpha[1] = alpha_d;
751:   VecRestoreArrayRead(Y,&dXarr);
752:   VecRestoreArrayWrite(X,&Xarr);

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

757:   alpha_p = alpha[2];
758:   alpha_d = alpha[3];

760:   /* X = X - alpha * Y */
761:   VecGetArrayWrite(X,&Xarr);
762:   VecGetArrayRead(Y,&dXarr);
763:   for (i=0; i<pdipm->nx; i++) Xarr[i] -= alpha_p * dXarr[i];
764:   for (i=0; i<pdipm->nce; i++) Xarr[i+pdipm->off_lambdae] -= alpha_d * dXarr[i+pdipm->off_lambdae];

766:   for (i=0; i<pdipm->nci; i++) {
767:     Xarr[i+pdipm->off_lambdai] -= alpha_d * dXarr[i+pdipm->off_lambdai];
768:     Xarr[i+pdipm->off_z]       -= alpha_p * dXarr[i+pdipm->off_z];
769:   }
770:   VecGetArrayWrite(tao->solution,&taosolarr);
771:   PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
772:   VecRestoreArrayWrite(tao->solution,&taosolarr);

774:   VecRestoreArrayWrite(X,&Xarr);
775:   VecRestoreArrayRead(Y,&dXarr);

777:   /* Update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
778:   if (pdipm->z) {
779:     VecDot(pdipm->z,pdipm->lambdai,&dot);
780:   } else dot = 0.0;

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

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

788:   tao->niter++;
789:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
790:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);

792:   (*tao->ops->convergencetest)(tao,tao->cnvP);
793:   if (tao->reason) {
794:     SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
795:   }
796:   return 0;
797: }

799: /*
800:    TaoSolve_PDIPM

802:    Input Parameter:
803:    tao - TAO context

805:    Output Parameter:
806:    tao - TAO context
807: */
808: PetscErrorCode TaoSolve_PDIPM(Tao tao)
809: {
810:   TAO_PDIPM          *pdipm = (TAO_PDIPM*)tao->data;
811:   SNESLineSearch     linesearch; /* SNESLineSearch context */
812:   Vec                dummy;


816:   /* Initialize all variables */
817:   TaoPDIPMInitializeSolution(tao);

819:   /* Set linesearch */
820:   SNESGetLineSearch(pdipm->snes,&linesearch);
821:   SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
822:   SNESLineSearchShellSetUserFunc(linesearch,SNESLineSearch_PDIPM,tao);
823:   SNESLineSearchSetFromOptions(linesearch);

825:   tao->reason = TAO_CONTINUE_ITERATING;

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

831:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
832:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
833:   VecDestroy(&dummy);
834:   (*tao->ops->convergencetest)(tao,tao->cnvP);
835:   if (tao->reason) {
836:     SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
837:   }

839:   while (tao->reason == TAO_CONTINUE_ITERATING) {
840:     SNESConvergedReason reason;
841:     SNESSolve(pdipm->snes,NULL,pdipm->X);

843:     /* Check SNES convergence */
844:     SNESGetConvergedReason(pdipm->snes,&reason);
845:     if (reason < 0) {
846:       PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
847:     }

849:     /* Check TAO convergence */
851:   }
852:   return 0;
853: }

855: /*
856:   TaoView_PDIPM - View PDIPM

858:    Input Parameter:
859:     tao - TAO object
860:     viewer - PetscViewer

862:    Output:
863: */
864: PetscErrorCode TaoView_PDIPM(Tao tao,PetscViewer viewer)
865: {
866:   TAO_PDIPM      *pdipm = (TAO_PDIPM *)tao->data;

868:   tao->constrained = PETSC_TRUE;
869:   PetscViewerASCIIPushTab(viewer);
870:   PetscViewerASCIIPrintf(viewer,"Number of prime=%D, Number of dual=%D\n",pdipm->Nx+pdipm->Nci,pdipm->Nce + pdipm->Nci);
871:   if (pdipm->kkt_pd) {
872:     PetscViewerASCIIPrintf(viewer,"KKT shifts deltaw=%g, deltac=%g\n",(double)pdipm->deltaw,(double)pdipm->deltac);
873:   }
874:   PetscViewerASCIIPopTab(viewer);
875:   return 0;
876: }

878: /*
879:    TaoSetup_PDIPM - Sets up tao and pdipm

881:    Input Parameter:
882:    tao - TAO object

884:    Output:   pdipm - initialized object
885: */
886: PetscErrorCode TaoSetup_PDIPM(Tao tao)
887: {
888:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
889:   PetscErrorCode    ierr;
890:   MPI_Comm          comm;
891:   PetscMPIInt       size;
892:   PetscInt          row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
893:   PetscInt          offset,*xa,*xb,i,j,rstart,rend;
894:   PetscScalar       one=1.0,neg_one=-1.0;
895:   const PetscInt    *cols,*rranges,*cranges,*aj,*ranges;
896:   const PetscScalar *aa,*Xarr;
897:   Mat               J,jac_equality_trans,jac_inequality_trans;
898:   Mat               Jce_xfixed_trans,Jci_xb_trans;
899:   PetscInt          *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];

901:   PetscObjectGetComm((PetscObject)tao,&comm);
902:   MPI_Comm_size(comm,&size);

904:   /* (1) Setup Bounds and create Tao vectors */
905:   TaoPDIPMSetUpBounds(tao);

907:   if (!tao->gradient) {
908:     VecDuplicate(tao->solution,&tao->gradient);
909:     VecDuplicate(tao->solution,&tao->stepdirection);
910:   }

912:   /* (2) Get sizes */
913:   /* Size of vector x - This is set by TaoSetSolution */
914:   VecGetSize(tao->solution,&pdipm->Nx);
915:   VecGetLocalSize(tao->solution,&pdipm->nx);

917:   /* Size of equality constraints and vectors */
918:   if (tao->constraints_equality) {
919:     VecGetSize(tao->constraints_equality,&pdipm->Ng);
920:     VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
921:   } else {
922:     pdipm->ng = pdipm->Ng = 0;
923:   }

925:   pdipm->nce = pdipm->ng + pdipm->nxfixed;
926:   pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;

928:   /* Size of inequality constraints and vectors */
929:   if (tao->constraints_inequality) {
930:     VecGetSize(tao->constraints_inequality,&pdipm->Nh);
931:     VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
932:   } else {
933:     pdipm->nh = pdipm->Nh = 0;
934:   }

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

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

943:   /* (3) Offsets for subvectors */
944:   pdipm->off_lambdae = pdipm->nx;
945:   pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
946:   pdipm->off_z       = pdipm->off_lambdai + pdipm->nci;

948:   /* (4) Create vectors and subvectors */
949:   /* Ce and Ci vectors */
950:   VecCreate(comm,&pdipm->ce);
951:   VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
952:   VecSetFromOptions(pdipm->ce);

954:   VecCreate(comm,&pdipm->ci);
955:   VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
956:   VecSetFromOptions(pdipm->ci);

958:   /* X=[x; lambdae; lambdai; z] for the big KKT system */
959:   VecCreate(comm,&pdipm->X);
960:   VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
961:   VecSetFromOptions(pdipm->X);

963:   /* Subvectors; they share local arrays with X */
964:   VecGetArrayRead(pdipm->X,&Xarr);
965:   /* x shares local array with X.x */
966:   if (pdipm->Nx) {
967:     VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
968:   }

970:   /* lambdae shares local array with X.lambdae */
971:   if (pdipm->Nce) {
972:     VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
973:   }

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

979:     VecCreate(comm,&pdipm->lambdae_xfixed);
980:     VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
981:     VecSetFromOptions(pdipm->lambdae_xfixed);
982:   }

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

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

992:   /* tao->DI which shares local array with X.lambdai_h */
993:   if (pdipm->Nh) {
994:     VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
995:   }
996:   VecCreate(comm,&pdipm->lambdai_xb);
997:   VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
998:   VecSetFromOptions(pdipm->lambdai_xb);

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

1002:   /* (5) Create Jacobians Jce_xfixed and Jci */
1003:   /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
1004:   if (pdipm->Nxfixed) {
1005:     /* Create Jce_xfixed */
1006:     MatCreate(comm,&pdipm->Jce_xfixed);
1007:     MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1008:     MatSetFromOptions(pdipm->Jce_xfixed);
1009:     MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
1010:     MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);

1012:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
1013:     ISGetIndices(pdipm->isxfixed,&cols);
1014:     k = 0;
1015:     for (row = Jcrstart; row < Jcrend; row++) {
1016:       MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
1017:       k++;
1018:     }
1019:     ISRestoreIndices(pdipm->isxfixed, &cols);
1020:     MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1021:     MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1022:   }

1024:   /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
1025:   MatCreate(comm,&pdipm->Jci_xb);
1026:   MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1027:   MatSetFromOptions(pdipm->Jci_xb);
1028:   MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
1029:   MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);

1031:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
1032:   offset = Jcrstart;
1033:   if (pdipm->Nxub) {
1034:     /* Add xub to Jci_xb */
1035:     ISGetIndices(pdipm->isxub,&cols);
1036:     k = 0;
1037:     for (row = offset; row < offset + pdipm->nxub; row++) {
1038:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1039:       k++;
1040:     }
1041:     ISRestoreIndices(pdipm->isxub, &cols);
1042:   }

1044:   if (pdipm->Nxlb) {
1045:     /* Add xlb to Jci_xb */
1046:     ISGetIndices(pdipm->isxlb,&cols);
1047:     k = 0;
1048:     offset += pdipm->nxub;
1049:     for (row = offset; row < offset + pdipm->nxlb; row++) {
1050:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
1051:       k++;
1052:     }
1053:     ISRestoreIndices(pdipm->isxlb, &cols);
1054:   }

1056:   /* Add xbox to Jci_xb */
1057:   if (pdipm->Nxbox) {
1058:     ISGetIndices(pdipm->isxbox,&cols);
1059:     k = 0;
1060:     offset += pdipm->nxlb;
1061:     for (row = offset; row < offset + pdipm->nxbox; row++) {
1062:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1063:       tmp = row + pdipm->nxbox;
1064:       MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
1065:       k++;
1066:     }
1067:     ISRestoreIndices(pdipm->isxbox, &cols);
1068:   }

1070:   MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1071:   MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1072:   /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */

1074:   /* (6) Set up ISs for PC Fieldsplit */
1075:   if (pdipm->solve_reduced_kkt) {
1076:     PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
1077:     for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1078:     for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;

1080:     ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
1081:     ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
1082:   }

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

1087:   /* Get rstart of KKT matrix */
1088:   MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
1089:   rstart -= pdipm->n;

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

1093:   PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
1094:   MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
1095:   MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
1096:   MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);

1098:   MatGetOwnershipRanges(tao->hessian,&rranges);
1099:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

1101:   if (pdipm->Ng) {
1102:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
1103:     MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
1104:   }
1105:   if (pdipm->Nh) {
1106:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
1107:     MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
1108:   }

1110:   /* Count dnz,onz for preallocation of KKT matrix */
1111:   jac_equality_trans   = pdipm->jac_equality_trans;
1112:   jac_inequality_trans = pdipm->jac_inequality_trans;
1113:   nce_all = pdipm->nce_all;

1115:   if (pdipm->Nxfixed) {
1116:     MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
1117:   }
1118:   MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);

1120:   MatPreallocateInitialize(comm,pdipm->n,pdipm->n,dnz,onz);

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

1126:   /* Insert tao->hessian */
1127:   MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
1128:   for (i=0; i<pdipm->nx; i++) {
1129:     row = rstart + i;

1131:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1132:     proc = 0;
1133:     for (j=0; j < nc; j++) {
1134:       while (aj[j] >= cranges[proc+1]) proc++;
1135:       col = aj[j] - cranges[proc] + Jranges[proc];
1136:       MatPreallocateSet(row,1,&col,dnz,onz);
1137:     }
1138:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);

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

1155:     /* Insert Jce_xfixed^T' */
1156:     if (pdipm->nxfixed) {
1157:       MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1158:       MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1159:       proc = 0;
1160:       for (j=0; j < nc; j++) {
1161:         /* find row ownership of */
1162:         while (aj[j] >= ranges[proc+1]) proc++;
1163:         nx_all = rranges[proc+1] - rranges[proc];
1164:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1165:         MatPreallocateSet(row,1,&col,dnz,onz);
1166:       }
1167:       MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1168:     }

1170:     if (pdipm->nh) {
1171:       /* Insert -grad h' */
1172:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1173:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1174:       proc = 0;
1175:       for (j=0; j < nc; j++) {
1176:         /* find row ownership of */
1177:         while (aj[j] >= ranges[proc+1]) proc++;
1178:         nx_all = rranges[proc+1] - rranges[proc];
1179:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1180:         MatPreallocateSet(row,1,&col,dnz,onz);
1181:       }
1182:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1183:     }

1185:     /* Insert Jci_xb^T' */
1186:     MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1187:     MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1188:     proc = 0;
1189:     for (j=0; j < nc; j++) {
1190:       /* find row ownership of */
1191:       while (aj[j] >= ranges[proc+1]) proc++;
1192:       nx_all = rranges[proc+1] - rranges[proc];
1193:       col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1194:       MatPreallocateSet(row,1,&col,dnz,onz);
1195:     }
1196:     MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1197:   }

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

1205:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1206:       proc = 0;
1207:       for (j=0; j < nc; j++) {
1208:         while (aj[j] >= cranges[proc+1]) proc++;
1209:         col = aj[j] - cranges[proc] + Jranges[proc];
1210:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1211:       }
1212:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1213:     }
1214:   }
1215:   /* Jce_xfixed */
1216:   if (pdipm->Nxfixed) {
1217:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1218:     for (i=0; i < (pdipm->nce - pdipm->ng); i++) {
1219:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

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

1224:       proc = 0;
1225:       j    = 0;
1226:       while (cols[j] >= cranges[proc+1]) proc++;
1227:       col = cols[j] - cranges[proc] + Jranges[proc];
1228:       MatPreallocateSet(row,1,&col,dnz,onz);
1229:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1230:     }
1231:   }

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

1239:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1240:       proc = 0;
1241:       for (j=0; j < nc; j++) {
1242:         while (aj[j] >= cranges[proc+1]) proc++;
1243:         col = aj[j] - cranges[proc] + Jranges[proc];
1244:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1245:       }
1246:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1247:     }
1248:     /* I */
1249:     for (i=0; i < pdipm->nh; i++) {
1250:       row = rstart + pdipm->off_lambdai + i;
1251:       col = rstart + pdipm->off_z + i;
1252:       MatPreallocateSet(row,1,&col,dnz,onz);
1253:     }
1254:   }

1256:   /* Jci_xb */
1257:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1258:   for (i=0; i < (pdipm->nci - pdipm->nh); i++) {
1259:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1261:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1263:     proc = 0;
1264:     for (j=0; j < nc; j++) {
1265:       while (cols[j] >= cranges[proc+1]) proc++;
1266:       col = cols[j] - cranges[proc] + Jranges[proc];
1267:       MatPreallocateSet(row,1,&col,dnz,onz);
1268:     }
1269:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1270:     /* I */
1271:     col = rstart + pdipm->off_z + pdipm->nh + i;
1272:     MatPreallocateSet(row,1,&col,dnz,onz);
1273:   }

1275:   /* 4-th Row block of KKT matrix: Z and Ci */
1276:   for (i=0; i < pdipm->nci; i++) {
1277:     row     = rstart + pdipm->off_z + i;
1278:     cols1[0] = rstart + pdipm->off_lambdai + i;
1279:     cols1[1] = row;
1280:     MatPreallocateSet(row,2,cols1,dnz,onz);
1281:   }

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

1286:   /* Create KKT matrix */
1287:   MatCreate(comm,&J);
1288:   MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1289:   MatSetFromOptions(J);
1290:   MatSeqAIJSetPreallocation(J,0,dnz);
1291:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1292:   MatPreallocateFinalize(dnz,onz);
1293:   pdipm->K = J;

1295:   /* (8) Insert constant entries to  K */
1296:   /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1297:   MatGetOwnershipRange(J,&rstart,&rend);
1298:   for (i=rstart; i<rend; i++) {
1299:     MatSetValue(J,i,i,0.0,INSERT_VALUES);
1300:   }
1301:   /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1302:   if (pdipm->kkt_pd) {
1303:       for (i=0; i<pdipm->nh; i++) {
1304:         row  = rstart + i;
1305:         MatSetValue(J,row,row,pdipm->deltaw,INSERT_VALUES);
1306:       }
1307:   }

1309:   /* Row block of K: [ grad Ce, 0, 0, 0] */
1310:   if (pdipm->Nxfixed) {
1311:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1312:     for (i=0; i < (pdipm->nce - pdipm->ng); i++) {
1313:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1315:       MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1316:       proc = 0;
1317:       for (j=0; j < nc; j++) {
1318:         while (cols[j] >= cranges[proc+1]) proc++;
1319:         col = cols[j] - cranges[proc] + Jranges[proc];
1320:         MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1321:         MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1322:       }
1323:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1324:     }
1325:   }

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

1332:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1333:     proc = 0;
1334:     for (j=0; j < nc; j++) {
1335:       while (cols[j] >= cranges[proc+1]) proc++;
1336:       col = cols[j] - cranges[proc] + Jranges[proc];
1337:       MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1338:       MatSetValue(J,row,col,-aa[j],INSERT_VALUES);
1339:     }
1340:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);

1342:     col = rstart + pdipm->off_z + pdipm->nh + i;
1343:     MatSetValue(J,row,col,1,INSERT_VALUES);
1344:   }

1346:   for (i=0; i < pdipm->nh; i++) {
1347:     row = rstart + pdipm->off_lambdai + i;
1348:     col = rstart + pdipm->off_z + i;
1349:     MatSetValue(J,row,col,1,INSERT_VALUES);
1350:   }

1352:   /* Row block of K: [ 0, 0, I, ...] */
1353:   for (i=0; i < pdipm->nci; i++) {
1354:     row = rstart + pdipm->off_z + i;
1355:     col = rstart + pdipm->off_lambdai + i;
1356:     MatSetValue(J,row,col,1,INSERT_VALUES);
1357:   }

1359:   if (pdipm->Nxfixed) {
1360:     MatDestroy(&Jce_xfixed_trans);
1361:   }
1362:   MatDestroy(&Jci_xb_trans);
1363:   PetscFree3(ng_all,nh_all,Jranges);

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

1369:   if (pdipm->solve_reduced_kkt) {
1370:     PC pc;
1371:     KSPGetPC(tao->ksp,&pc);
1372:     PCSetType(pc,PCFIELDSPLIT);
1373:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1374:     PCFieldSplitSetIS(pc,"2",pdipm->is2);
1375:     PCFieldSplitSetIS(pc,"1",pdipm->is1);
1376:   }
1377:   SNESSetFromOptions(pdipm->snes);

1379:   /* (10) Setup PCPreSolve() for pdipm->solve_symmetric_kkt */
1380:   if (pdipm->solve_symmetric_kkt) {
1381:     KSP       ksp;
1382:     PC        pc;
1383:     PetscBool isCHOL;
1384:     SNESGetKSP(pdipm->snes,&ksp);
1385:     KSPGetPC(ksp,&pc);
1386:     PCSetPreSolve(pc,PCPreSolve_PDIPM);

1388:     PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);
1389:     if (isCHOL) {
1390:       Mat        Factor;
1391:       PetscBool  isMUMPS;
1392:       PCFactorGetMatrix(pc,&Factor);
1393:       PetscObjectTypeCompare((PetscObject)Factor,"mumps",&isMUMPS);
1394:       if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
1395: #if defined(PETSC_HAVE_MUMPS)
1396:         MatMumpsSetIcntl(Factor,24,1); /* detection of null pivot rows */
1397:         if (size > 1) {
1398:           MatMumpsSetIcntl(Factor,13,1); /* parallelism of the root node (enable ScaLAPACK) and its splitting */
1399:         }
1400: #else
1401:         SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires external package MUMPS");
1402: #endif
1403:       }
1404:     }
1405:   }
1406:   return 0;
1407: }

1409: /*
1410:    TaoDestroy_PDIPM - Destroys the pdipm object

1412:    Input:
1413:    full pdipm

1415:    Output:
1416:    Destroyed pdipm
1417: */
1418: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1419: {
1420:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1422:   /* Freeing Vectors assocaiated with KKT (X) */
1423:   VecDestroy(&pdipm->x); /* Solution x */
1424:   VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1425:   VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1426:   VecDestroy(&pdipm->z);       /* Slack variables */
1427:   VecDestroy(&pdipm->X);       /* Big KKT system vector [x; lambdae; lambdai; z] */

1429:   /* work vectors */
1430:   VecDestroy(&pdipm->lambdae_xfixed);
1431:   VecDestroy(&pdipm->lambdai_xb);

1433:   /* Legrangian equality and inequality Vec */
1434:   VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1435:   VecDestroy(&pdipm->ci); /* Vec of inequality constraints */

1437:   /* Matrices */
1438:   MatDestroy(&pdipm->Jce_xfixed);
1439:   MatDestroy(&pdipm->Jci_xb)); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb; J(nxbx)] */
1440:   MatDestroy(&pdipm->K);

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

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

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

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

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

1463:   if (pdipm->solve_reduced_kkt) {
1464:     ISDestroy(&pdipm->is1);
1465:     ISDestroy(&pdipm->is2);
1466:   }

1468:   /* SNES */
1469:   SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1470:   PetscFree(pdipm->nce_all);
1471:   MatDestroy(&pdipm->jac_equality_trans);
1472:   MatDestroy(&pdipm->jac_inequality_trans);

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

1477:   /* Destroy Dual */
1478:   VecDestroy(&tao->DE); /* equality dual */
1479:   VecDestroy(&tao->DI); /* dinequality dual */
1480:   return 0;
1481: }

1483: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1484: {
1485:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1487:   PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1488:   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);
1489:   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);
1490:   PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1491:   PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1492:   PetscOptionsBool("-tao_pdipm_symmetric_kkt","Solve non reduced symmetric KKT system",NULL,pdipm->solve_symmetric_kkt,&pdipm->solve_symmetric_kkt,NULL);
1493:   PetscOptionsBool("-tao_pdipm_kkt_shift_pd","Add shifts to make KKT matrix positive definite",NULL,pdipm->kkt_pd,&pdipm->kkt_pd,NULL);
1494:   PetscOptionsTail();
1495:   return 0;
1496: }

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

1501:   Option Database Keys:
1502: +   -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1503: .   -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1504: .   -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1505: .   -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1506: -   -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite

1508:   Level: beginner
1509: M*/
1510: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1511: {
1512:   TAO_PDIPM      *pdipm;

1514:   tao->ops->setup          = TaoSetup_PDIPM;
1515:   tao->ops->solve          = TaoSolve_PDIPM;
1516:   tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1517:   tao->ops->view           = TaoView_PDIPM;
1518:   tao->ops->destroy        = TaoDestroy_PDIPM;

1520:   PetscNewLog(tao,&pdipm);
1521:   tao->data = (void*)pdipm;

1523:   pdipm->nx      = pdipm->Nx      = 0;
1524:   pdipm->nxfixed = pdipm->Nxfixed = 0;
1525:   pdipm->nxlb    = pdipm->Nxlb    = 0;
1526:   pdipm->nxub    = pdipm->Nxub    = 0;
1527:   pdipm->nxbox   = pdipm->Nxbox   = 0;
1528:   pdipm->nxfree  = pdipm->Nxfree  = 0;

1530:   pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1531:   pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1532:   pdipm->n  = pdipm->N  = 0;
1533:   pdipm->mu = 1.0;
1534:   pdipm->mu_update_factor = 0.1;

1536:   pdipm->deltaw     = 0.0;
1537:   pdipm->lastdeltaw = 3*1.e-4;
1538:   pdipm->deltac     = 0.0;
1539:   pdipm->kkt_pd     = PETSC_FALSE;

1541:   pdipm->push_init_slack     = 1.0;
1542:   pdipm->push_init_lambdai   = 1.0;
1543:   pdipm->solve_reduced_kkt   = PETSC_FALSE;
1544:   pdipm->solve_symmetric_kkt = PETSC_TRUE;

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

1550:   SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1551:   SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1552:   SNESGetKSP(pdipm->snes,&tao->ksp);
1553:   PetscObjectReference((PetscObject)tao->ksp);
1554:   KSPSetApplicationContext(tao->ksp,(void *)tao);
1555:   return 0;
1556: }