Actual source code: pdipm.c

petsc-main 2021-04-20
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/constrained/impls/ipm/pdipm.h>
  3: #include <petscsnes.h>
  4: #include <petsc/private/pcimpl.h>
  5: #include <petsc/private/matimpl.h>

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

 10:    Collective on tao

 12:    Input Parameter:
 13: +  tao - solver context
 14: -  x - vector at which all objects to be evaluated

 16:    Level: beginner

 18: .seealso: TaoPDIPMUpdateConstraints(), TaoPDIPMSetUpBounds()
 19: */
 20: static PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)
 21: {
 23:   TAO_PDIPM      *pdipm=(TAO_PDIPM*)tao->data;

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

 29:   /* Equality constraints and Jacobian */
 30:   if (pdipm->Ng) {
 31:     TaoComputeEqualityConstraints(tao,x,tao->constraints_equality);
 32:     TaoComputeJacobianEquality(tao,x,tao->jacobian_equality,tao->jacobian_equality_pre);
 33:   }

 35:   /* Inequality constraints and Jacobian */
 36:   if (pdipm->Nh) {
 37:     TaoComputeInequalityConstraints(tao,x,tao->constraints_inequality);
 38:     TaoComputeJacobianInequality(tao,x,tao->jacobian_inequality,tao->jacobian_inequality_pre);
 39:   }
 40:   return(0);
 41: }

 43: /*
 44:   TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x

 46:   Collective

 48:   Input Parameter:
 49: + tao - Tao context
 50: - x - vector at which constraints to be evaluted

 52:    Level: beginner

 54: .seealso: TaoPDIPMEvaluateFunctionsAndJacobians()
 55: */
 56: static PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao,Vec x)
 57: {
 58:   PetscErrorCode    ierr;
 59:   TAO_PDIPM         *pdipm=(TAO_PDIPM*)tao->data;
 60:   PetscInt          i,offset,offset1,k,xstart;
 61:   PetscScalar       *carr;
 62:   const PetscInt    *ubptr,*lbptr,*bxptr,*fxptr;
 63:   const PetscScalar *xarr,*xuarr,*xlarr,*garr,*harr;

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

 68:   VecGetArrayRead(x,&xarr);
 69:   VecGetArrayRead(tao->XU,&xuarr);
 70:   VecGetArrayRead(tao->XL,&xlarr);

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

 75:   if (pdipm->Ng) {
 76:     /* (1.a) Inserting updated g(x) */
 77:     VecGetArrayRead(tao->constraints_equality,&garr);
 78:     PetscMemcpy(carr,garr,pdipm->ng*sizeof(PetscScalar));
 79:     VecRestoreArrayRead(tao->constraints_equality,&garr);
 80:   }

 82:   /* (1.b) Update xfixed */
 83:   if (pdipm->Nxfixed) {
 84:     offset = pdipm->ng;
 85:     ISGetIndices(pdipm->isxfixed,&fxptr); /* global indices in x */
 86:     for (k=0;k < pdipm->nxfixed;k++){
 87:       i = fxptr[k]-xstart;
 88:       carr[offset + k] = xarr[i] - xuarr[i];
 89:     }
 90:   }
 91:   VecRestoreArrayWrite(pdipm->ce,&carr);

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

 96:   if (pdipm->Nh) {
 97:     /* (2.a) Inserting updated h(x) */
 98:     VecGetArrayRead(tao->constraints_inequality,&harr);
 99:     PetscMemcpy(carr,harr,pdipm->nh*sizeof(PetscScalar));
100:     VecRestoreArrayRead(tao->constraints_inequality,&harr);
101:   }

103:   /* (2.b) Update xub */
104:   offset = pdipm->nh;
105:   if (pdipm->Nxub) {
106:     ISGetIndices(pdipm->isxub,&ubptr);
107:     for (k=0; k<pdipm->nxub; k++){
108:       i = ubptr[k]-xstart;
109:       carr[offset + k] = xuarr[i] - xarr[i];
110:     }
111:   }

113:   if (pdipm->Nxlb) {
114:     /* (2.c) Update xlb */
115:     offset += pdipm->nxub;
116:     ISGetIndices(pdipm->isxlb,&lbptr); /* global indices in x */
117:     for (k=0; k<pdipm->nxlb; k++){
118:       i = lbptr[k]-xstart;
119:       carr[offset + k] = xarr[i] - xlarr[i];
120:     }
121:   }

123:   if (pdipm->Nxbox) {
124:     /* (2.d) Update xbox */
125:     offset += pdipm->nxlb;
126:     offset1 = offset + pdipm->nxbox;
127:     ISGetIndices(pdipm->isxbox,&bxptr); /* global indices in x */
128:     for (k=0; k<pdipm->nxbox; k++){
129:       i = bxptr[k]-xstart; /* local indices in x */
130:       carr[offset+k]  = xuarr[i] - xarr[i];
131:       carr[offset1+k] = xarr[i]  - xlarr[i];
132:     }
133:   }
134:   VecRestoreArrayWrite(pdipm->ci,&carr);

136:   /* Restoring Vectors */
137:   VecRestoreArrayRead(x,&xarr);
138:   VecRestoreArrayRead(tao->XU,&xuarr);
139:   VecRestoreArrayRead(tao->XL,&xlarr);
140:   return(0);
141: }

143: /*
144:    TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x

146:    Collective

148:    Input Parameter:
149: .  tao - holds pdipm and XL & XU

151:    Level: beginner

153: .seealso: TaoPDIPMUpdateConstraints
154: */
155: static PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
156: {
157:   PetscErrorCode    ierr;
158:   TAO_PDIPM         *pdipm=(TAO_PDIPM*)tao->data;
159:   const PetscScalar *xl,*xu;
160:   PetscInt          n,*ixlb,*ixub,*ixfixed,*ixfree,*ixbox,i,low,high,idx;
161:   MPI_Comm          comm;
162:   PetscInt          sendbuf[5],recvbuf[5];

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

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

171:   VecGetOwnershipRange(tao->XL,&low,&high);
172:   VecGetArrayRead(tao->XL,&xl);
173:   VecGetArrayRead(tao->XU,&xu);
174:   for (i=0; i<n; i++) {
175:     idx = low + i;
176:     if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
177:       if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
178:         ixfixed[pdipm->nxfixed++]  = idx;
179:       } else ixbox[pdipm->nxbox++] = idx;
180:     } else {
181:       if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
182:         ixlb[pdipm->nxlb++] = idx;
183:       } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
184:         ixub[pdipm->nxlb++] = idx;
185:       } else  ixfree[pdipm->nxfree++] = idx;
186:     }
187:   }
188:   VecRestoreArrayRead(tao->XL,&xl);
189:   VecRestoreArrayRead(tao->XU,&xu);

191:   PetscObjectGetComm((PetscObject)tao,&comm);
192:   sendbuf[0] = pdipm->nxlb;
193:   sendbuf[1] = pdipm->nxub;
194:   sendbuf[2] = pdipm->nxfixed;
195:   sendbuf[3] = pdipm->nxbox;
196:   sendbuf[4] = pdipm->nxfree;

198:   MPI_Allreduce(sendbuf,recvbuf,5,MPIU_INT,MPI_SUM,comm);
199:   pdipm->Nxlb    = recvbuf[0];
200:   pdipm->Nxub    = recvbuf[1];
201:   pdipm->Nxfixed = recvbuf[2];
202:   pdipm->Nxbox   = recvbuf[3];
203:   pdipm->Nxfree  = recvbuf[4];

205:   if (pdipm->Nxlb) {
206:     ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);
207:   }
208:   if (pdipm->Nxub) {
209:     ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);
210:   }
211:   if (pdipm->Nxfixed) {
212:     ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);
213:   }
214:   if (pdipm->Nxbox) {
215:     ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);
216:   }
217:   if (pdipm->Nxfree) {
218:     ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);
219:   }
220:   PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);
221:   return(0);
222: }

224: /*
225:    TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
226:    X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
227:      four subvectors need to be initialized and its values copied over to X. Instead
228:      of copying, we use VecPlace/ResetArray functions to share the memory locations for
229:      X and the subvectors

231:    Collective

233:    Input Parameter:
234: .  tao - Tao context

236:    Level: beginner
237: */
238: static PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
239: {
240:   PetscErrorCode    ierr;
241:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
242:   PetscScalar       *Xarr,*z,*lambdai;
243:   PetscInt          i;
244:   const PetscScalar *xarr,*h;

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

249:   /* Set Initialize X.x = tao->solution */
250:   VecGetArrayRead(tao->solution,&xarr);
251:   PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));
252:   VecRestoreArrayRead(tao->solution,&xarr);

254:   /* Initialize X.lambdae = 0.0 */
255:   if (pdipm->lambdae) {
256:     VecSet(pdipm->lambdae,0.0);
257:   }

259:   /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
260:   if (pdipm->Nci) {
261:     VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
262:     VecSet(pdipm->z,pdipm->push_init_slack);

264:     /* Additional modification for X.lambdai and X.z */
265:     VecGetArrayWrite(pdipm->lambdai,&lambdai);
266:     VecGetArrayWrite(pdipm->z,&z);
267:     if (pdipm->Nh) {
268:       VecGetArrayRead(tao->constraints_inequality,&h);
269:       for (i=0; i < pdipm->nh; i++) {
270:         if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
271:         if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
272:       }
273:       VecRestoreArrayRead(tao->constraints_inequality,&h);
274:     }
275:     VecRestoreArrayWrite(pdipm->lambdai,&lambdai);
276:     VecRestoreArrayWrite(pdipm->z,&z);
277:   }

279:   VecRestoreArrayWrite(pdipm->X,&Xarr);
280:   return(0);
281: }

283: /*
284:    TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X

286:    Input Parameter:
287:    snes - SNES context
288:    X - KKT Vector
289:    *ctx - pdipm context

291:    Output Parameter:
292:    J - Hessian matrix
293:    Jpre - Preconditioner
294: */
295: static PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
296: {
297:   PetscErrorCode    ierr;
298:   Tao               tao=(Tao)ctx;
299:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
300:   PetscInt          i,row,cols[2],Jrstart,rjstart,nc,j;
301:   const PetscInt    *aj,*ranges,*Jranges,*rranges,*cranges;
302:   const PetscScalar *Xarr,*aa;
303:   PetscScalar       vals[2];
304:   PetscInt          proc,nx_all,*nce_all=pdipm->nce_all;
305:   MPI_Comm          comm;
306:   PetscMPIInt       rank,size;
307:   Mat               jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;

310:   PetscObjectGetComm((PetscObject)snes,&comm);
311:   MPI_Comm_rank(comm,&rank);
312:   MPI_Comm_rank(comm,&size);

314:   MatGetOwnershipRanges(Jpre,&Jranges);
315:   MatGetOwnershipRange(Jpre,&Jrstart,NULL);
316:   MatGetOwnershipRangesColumn(tao->hessian,&rranges);
317:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

319:   VecGetArrayRead(X,&Xarr);

321:   /* (1) insert Z and Ci to the 4th block of Jpre -- overwrite existing values */
322:   if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
323:     vals[0] = 1.0;
324:     for (i=0; i < pdipm->nci; i++) {
325:         row     = Jrstart + pdipm->off_z + i;
326:         cols[0] = Jrstart + pdipm->off_lambdai + i;
327:         cols[1] = row;
328:         vals[1] = Xarr[pdipm->off_lambdai + i]/Xarr[pdipm->off_z + i];
329:         MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
330:     }
331:   } else {
332:     for (i=0; i < pdipm->nci; i++) {
333:       row     = Jrstart + pdipm->off_z + i;
334:       cols[0] = Jrstart + pdipm->off_lambdai + i;
335:       cols[1] = row;
336:       vals[0] = Xarr[pdipm->off_z + i];
337:       vals[1] = Xarr[pdipm->off_lambdai + i];
338:       MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
339:     }
340:   }

342:   /* (2) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
343:   if (pdipm->Ng) {
344:     MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
345:     for (i=0; i<pdipm->ng; i++){
346:       row = Jrstart + pdipm->off_lambdae + i;

348:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
349:       proc = 0;
350:       for (j=0; j < nc; j++) {
351:         while (aj[j] >= cranges[proc+1]) proc++;
352:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
353:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
354:       }
355:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
356:       if (pdipm->kkt_pd) {
357:         /* add shift \delta_c */
358:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
359:       }
360:     }
361:   }

363:   /* (3) insert 3nd row block of Jpre: [ -grad h, 0, deltac, I] */
364:   if (pdipm->Nh) {
365:     MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
366:     for (i=0; i < pdipm->nh; i++){
367:       row = Jrstart + pdipm->off_lambdai + i;
368:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
369:       proc = 0;
370:       for (j=0; j < nc; j++) {
371:         while (aj[j] >= cranges[proc+1]) proc++;
372:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
373:         MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
374:       }
375:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
376:       if (pdipm->kkt_pd) {
377:         /* add shift \delta_c */
378:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
379:       }
380:     }
381:   }

383:   /* (4) insert 1st row block of Jpre: [Wxx, grad g', -grad h', 0] */
384:   if (pdipm->Ng) { /* grad g' */
385:     MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
386:   }
387:   if (pdipm->Nh) { /* grad h' */
388:     MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
389:   }

391:   VecPlaceArray(pdipm->x,Xarr);
392:   TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
393:   VecResetArray(pdipm->x);

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

399:     /* insert Wxx = fxx + ... -- provided by user */
400:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
401:     proc = 0;
402:     for (j=0; j < nc; j++) {
403:       while (aj[j] >= cranges[proc+1]) proc++;
404:       cols[0] = aj[j] - cranges[proc] + Jranges[proc];
405:       if (row == cols[0] && pdipm->kkt_pd) {
406:         /* add shift deltaw to Wxx component */
407:         MatSetValue(Jpre,row,cols[0],aa[j]+pdipm->deltaw,INSERT_VALUES);
408:       } else {
409:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
410:       }
411:     }
412:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);

414:     /* insert grad g' */
415:     if (pdipm->ng) {
416:       MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
417:       MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
418:       proc = 0;
419:       for (j=0; j < nc; j++) {
420:         /* find row ownership of */
421:         while (aj[j] >= ranges[proc+1]) proc++;
422:         nx_all = rranges[proc+1] - rranges[proc];
423:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
424:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
425:       }
426:       MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
427:     }

429:     /* insert -grad h' */
430:     if (pdipm->nh) {
431:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
432:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
433:       proc = 0;
434:       for (j=0; j < nc; j++) {
435:         /* find row ownership of */
436:         while (aj[j] >= ranges[proc+1]) proc++;
437:         nx_all = rranges[proc+1] - rranges[proc];
438:         cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
439:         MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
440:       }
441:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
442:     }
443:   }
444:   VecRestoreArrayRead(X,&Xarr);

446:   /* (6) assemble Jpre and J */
447:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
448:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);

450:   if (J != Jpre) {
451:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
452:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
453:   }
454:   return(0);
455: }

457: /*
458:    TaoSnesFunction_PDIPM - Evaluate KKT function at X

460:    Input Parameter:
461:    snes - SNES context
462:    X - KKT Vector
463:    *ctx - pdipm

465:    Output Parameter:
466:    F - Updated Lagrangian vector
467: */
468: static PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
469: {
470:   PetscErrorCode    ierr;
471:   Tao               tao=(Tao)ctx;
472:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
473:   PetscScalar       *Farr;
474:   Vec               x,L1;
475:   PetscInt          i;
476:   const PetscScalar *Xarr,*carr,*zarr,*larr;

479:   VecSet(F,0.0);

481:   VecGetArrayRead(X,&Xarr);
482:   VecGetArrayWrite(F,&Farr);

484:   /* (0) Evaluate f, fx, gradG, gradH at X.x Note: pdipm->x is not changed below */
485:   x = pdipm->x;
486:   VecPlaceArray(x,Xarr);
487:   TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);

489:   /* Update ce, ci, and Jci at X.x */
490:   TaoPDIPMUpdateConstraints(tao,x);
491:   VecResetArray(x);

493:   /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
494:   L1 = pdipm->x;
495:   VecPlaceArray(L1,Farr); /* L1 = 0.0 */
496:   if (pdipm->Nci) {
497:     if (pdipm->Nh) {
498:       /* L1 += gradH'*DI. Note: tao->DI is not changed below */
499:       VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
500:       MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
501:       VecResetArray(tao->DI);
502:     }

504:     /* L1 += Jci_xb'*lambdai_xb */
505:     VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
506:     MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
507:     VecResetArray(pdipm->lambdai_xb);

509:     /* L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
510:     VecScale(L1,-1.0);
511:   }

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

516:   if (pdipm->Nce) {
517:     if (pdipm->Ng) {
518:       /* L1 += gradG'*DE. Note: tao->DE is not changed below */
519:       VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
520:       MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
521:       VecResetArray(tao->DE);
522:     }
523:     if (pdipm->Nxfixed) {
524:       /* L1 += Jce_xfixed'*lambdae_xfixed */
525:       VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
526:       MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
527:       VecResetArray(pdipm->lambdae_xfixed);
528:     }
529:   }
530:   VecResetArray(L1);

532:   /* (2) L2 = ce(x) */
533:   if (pdipm->Nce) {
534:     VecGetArrayRead(pdipm->ce,&carr);
535:     for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
536:     VecRestoreArrayRead(pdipm->ce,&carr);
537:   }

539:   if (pdipm->Nci) {
540:     if (pdipm->solve_symmetric_kkt) {
541:       /* (3) L3 = z - ci(x);
542:          (4) L4 = Lambdai * e - mu/z *e  */
543:       VecGetArrayRead(pdipm->ci,&carr);
544:       larr = Xarr+pdipm->off_lambdai;
545:       zarr = Xarr+pdipm->off_z;
546:       for (i=0; i<pdipm->nci; i++) {
547:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
548:         Farr[pdipm->off_z       + i] = larr[i] - pdipm->mu/zarr[i];
549:       }
550:       VecRestoreArrayRead(pdipm->ci,&carr);
551:     } else {
552:       /* (3) L3 = z - ci(x);
553:          (4) L4 = Z * Lambdai * e - mu * e  */
554:       VecGetArrayRead(pdipm->ci,&carr);
555:       larr = Xarr+pdipm->off_lambdai;
556:       zarr = Xarr+pdipm->off_z;
557:       for (i=0; i<pdipm->nci; i++) {
558:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
559:         Farr[pdipm->off_z       + i] = zarr[i]*larr[i] - pdipm->mu;
560:       }
561:       VecRestoreArrayRead(pdipm->ci,&carr);
562:     }
563:   }

565:   VecRestoreArrayRead(X,&Xarr);
566:   VecRestoreArrayWrite(F,&Farr);
567:   return(0);
568: }

570: /*
571:  Evaluate F(X), then compute tao->residual = norm2(F_x,F_z) and tao->cnorm = norm2(F_ce,F_ci)
572: */
573: static PetscErrorCode TaoSNESFunction_PDIPM_residual(SNES snes,Vec X,Vec F,void *ctx)
574: {
575:   PetscErrorCode    ierr;
576:   Tao               tao=(Tao)ctx;
577:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
578:   PetscScalar       *Farr,*tmparr;
579:   Vec               L1;
580:   PetscInt          i;
581:   PetscReal         res[2],cnorm[2];
582:   const PetscScalar *Xarr=NULL;

585:   TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);
586:   VecGetArrayWrite(F,&Farr);
587:   VecGetArrayRead(X,&Xarr);

589:   /* compute norm2(F_x), norm2(F_z) */
590:   L1 = pdipm->x;
591:   VecPlaceArray(L1,Farr);
592:   VecNorm(L1,NORM_2,&res[0]);
593:   VecResetArray(L1);

595:   if (pdipm->z) {
596:     if (pdipm->solve_symmetric_kkt) {
597:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
598:       if (pdipm->Nci) {
599:         VecGetArrayWrite(pdipm->z,&tmparr);
600:         for (i=0; i<pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
601:         VecRestoreArrayWrite(pdipm->z,&tmparr);
602:       }

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

606:       if (pdipm->Nci) {
607:         VecGetArrayWrite(pdipm->z,&tmparr);
608:         for (i=0; i<pdipm->nci; i++) {
609:           tmparr[i] /= Xarr[pdipm->off_z + i];
610:         }
611:         VecRestoreArrayWrite(pdipm->z,&tmparr);
612:       }
613:       VecResetArray(pdipm->z);
614:     } else { /* !solve_symmetric_kkt */
615:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
616:       VecNorm(pdipm->z,NORM_2,&res[1]);
617:       VecResetArray(pdipm->z);
618:     }
619:   } else res[1] = 0.0;

621:   tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);

623:   /* compute norm2(F_ce), norm2(F_ci) */
624:   if (pdipm->Nce) {
625:     VecPlaceArray(pdipm->ce,Farr+pdipm->off_lambdae);
626:     VecNorm(pdipm->ce,NORM_2,&cnorm[0]);
627:     VecResetArray(pdipm->ce);
628:   } else cnorm[0] = 0.0;

630:   VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
631:   VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
632:   VecResetArray(pdipm->ci);

634:   tao->cnorm  = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);

636:   VecRestoreArrayWrite(F,&Farr);
637:   VecRestoreArrayRead(X,&Xarr);
638:   return(0);
639: }

641: /*
642:   KKTAddShifts - Check the inertia of Cholesky factor of KKT matrix.
643:   If it does not match the numbers of prime and dual variables, add shifts to the KKT matrix.
644: */
645: static PetscErrorCode KKTAddShifts(Tao tao,SNES snes,Vec X)
646: {
648:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;
649:   KSP            ksp;
650:   PC             pc;
651:   PCType         ptype;
652:   Mat            Factor;
653:   PetscBool      isCHOL;
654:   PetscInt       nneg,nzero,npos;

657:   /* Get the inertia of Cholesky factor */
658:   SNESGetKSP(snes,&ksp);
659:   KSPGetPC(ksp,&pc);
660:   PCGetType(pc,&ptype);
661:   PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);

663:   if (isCHOL) {
664:     PetscMPIInt       size;
665:     PCFactorGetMatrix(pc,&Factor);
666:     MPI_Comm_size(PetscObjectComm((PetscObject)Factor),&size);
667:     if (Factor->ops->getinertia) {
668: #if defined(PETSC_HAVE_MUMPS)
669:       MatSolverType     stype;
670:       PetscBool         isMUMPS;
671:       PCFactorGetMatSolverType(pc,&stype);
672:       PetscStrcmp(stype, MATSOLVERMUMPS, &isMUMPS);
673:       if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
674:         MatMumpsSetIcntl(Factor,24,1);
675:         if (size > 1) {
676:           MatMumpsSetIcntl(Factor,13,1);
677:         }
678:       }
679: #else
680:       if (size > 1) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires external package MUMPS");
681: #endif
682:       MatGetInertia(Factor,&nneg,&nzero,&npos);

684:       if (npos < pdipm->Nx+pdipm->Nci) {
685:         pdipm->deltaw = PetscMax(pdipm->lastdeltaw/3, 1.e-4*PETSC_MACHINE_EPSILON);
686:         PetscInfo5(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);
687:         TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
688:         PCSetUp(pc);
689:         MatGetInertia(Factor,&nneg,&nzero,&npos);

691:         if (npos < pdipm->Nx+pdipm->Nci) {
692:           pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
693:           while (npos < pdipm->Nx+pdipm->Nci && pdipm->deltaw <= 1.e10) { /* increase deltaw */
694:             PetscInfo5(tao,"  deltaw=%g fails, MatInertia: nneg %D, nzero %D, npos %D(<%D)\n",(double)pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
695:             pdipm->deltaw = PetscMin(8*pdipm->deltaw,PetscPowReal(10,20));
696:             TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
697:             PCSetUp(pc);
698:             MatGetInertia(Factor,&nneg,&nzero,&npos);
699:           }

701:           if (pdipm->deltaw >= 1.e10) {
702:             SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_CONV_FAILED,"Reached maximum delta w will not converge, try different inital x0");
703:           }
704:           PetscInfo1(tao,"Updated deltaw %g\n",(double)pdipm->deltaw);
705:           pdipm->lastdeltaw = pdipm->deltaw;
706:           pdipm->deltaw     = 0.0;
707:         }
708:       }

710:       if (nzero) { /* Jacobian is singular */
711:         if (pdipm->deltac == 0.0) {
712:           pdipm->deltac = PETSC_SQRT_MACHINE_EPSILON;
713:         } else {
714:           pdipm->deltac = pdipm->deltac*PetscPowReal(pdipm->mu,.25);
715:         }
716:         PetscInfo4(tao,"Updated deltac=%g, MatInertia: nneg %D, nzero %D(!=0), npos %D\n",(double)pdipm->deltac,nneg,nzero,npos);
717:         TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
718:         PCSetUp(pc);
719:         MatGetInertia(Factor,&nneg,&nzero,&npos);
720:       }
721:     } else SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires an external package that supports MatGetInertia()");
722:   }
723:   return(0);
724: }

726: /*
727:   PCPreSolve_PDIPM -- called betwee MatFactorNumeric() and MatSolve()
728: */
729: PetscErrorCode PCPreSolve_PDIPM(PC pc,KSP ksp,Vec rhs,Vec x)
730: {
732:   Tao            tao;
733:   TAO_PDIPM      *pdipm;

736:   KSPGetApplicationContext(ksp,&tao);
737:   pdipm = (TAO_PDIPM*)tao->data;
738:   KKTAddShifts(tao,pdipm->snes,pdipm->X);
739:   return(0);
740: }

742: /*
743:    SNESLineSearch_PDIPM - Custom line search used with PDIPM.

745:    Collective on TAO

747:    Notes:
748:    This routine employs a simple backtracking line-search to keep
749:    the slack variables (z) and inequality constraints Lagrange multipliers
750:    (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
751:    alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
752:    slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
753:    are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
754:    is also updated as mu = mu + z'lambdai/Nci
755: */
756: static PetscErrorCode SNESLineSearch_PDIPM(SNESLineSearch linesearch,void *ctx)
757: {
758:   PetscErrorCode    ierr;
759:   Tao               tao=(Tao)ctx;
760:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
761:   SNES              snes;
762:   Vec               X,F,Y,W,G;
763:   PetscInt          i,iter;
764:   PetscReal         alpha_p=1.0,alpha_d=1.0,alpha[4];
765:   PetscScalar       *Xarr,*z,*lambdai,dot,*taosolarr;
766:   const PetscScalar *dXarr,*dz,*dlambdai;

769:   SNESLineSearchGetSNES(linesearch,&snes);
770:   SNESGetIterationNumber(snes,&iter);

772:   SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
773:   SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);

775:   VecGetArrayWrite(X,&Xarr);
776:   VecGetArrayRead(Y,&dXarr);
777:   z  = Xarr + pdipm->off_z;
778:   dz = dXarr + pdipm->off_z;
779:   for (i=0; i < pdipm->nci; i++) {
780:     if (z[i] - dz[i] < 0.0) {
781:       alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
782:     }
783:   }

785:   lambdai  = Xarr + pdipm->off_lambdai;
786:   dlambdai = dXarr + pdipm->off_lambdai;

788:   for (i=0; i<pdipm->nci; i++) {
789:     if (lambdai[i] - dlambdai[i] < 0.0) {
790:       alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
791:     }
792:   }

794:   alpha[0] = alpha_p;
795:   alpha[1] = alpha_d;
796:   VecRestoreArrayRead(Y,&dXarr);
797:   VecRestoreArrayWrite(X,&Xarr);

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

802:   alpha_p = alpha[2];
803:   alpha_d = alpha[3];

805:   VecGetArrayWrite(X,&Xarr);
806:   VecGetArrayRead(Y,&dXarr);
807:   for (i=0; i<pdipm->nx; i++) Xarr[i] -= alpha_p * dXarr[i];

809:   for (i=0; i<pdipm->nce; i++) {
810:     Xarr[i+pdipm->off_lambdae] -= alpha_d * dXarr[i+pdipm->off_lambdae];
811:   }

813:   for (i=0; i<pdipm->nci; i++) {
814:     Xarr[i+pdipm->off_lambdai] -= alpha_d * dXarr[i+pdipm->off_lambdai];
815:     Xarr[i+pdipm->off_z]       -= alpha_p * dXarr[i+pdipm->off_z];
816:   }
817:   VecGetArrayWrite(tao->solution,&taosolarr);
818:   PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
819:   VecRestoreArrayWrite(tao->solution,&taosolarr);

821:   VecRestoreArrayWrite(X,&Xarr);
822:   VecRestoreArrayRead(Y,&dXarr);

824:   /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
825:   if (pdipm->z) {
826:     VecDot(pdipm->z,pdipm->lambdai,&dot);
827:   } else dot = 0.0;

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

832:   /* Update F; get tao->residual and tao->cnorm */
833:   TaoSNESFunction_PDIPM_residual(snes,X,F,(void*)tao);
834:   SNESLineSearchComputeNorms(linesearch); /* must call this func, do not know why */

836:   tao->niter++;
837:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
838:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);

840:   (*tao->ops->convergencetest)(tao,tao->cnvP);
841:   if (tao->reason) {
842:     SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
843:   }
844:   return(0);
845: }

847: /*
848:    TaoSolve_PDIPM

850:    Input Parameter:
851:    tao - TAO context

853:    Output Parameter:
854:    tao - TAO context
855: */
856: PetscErrorCode TaoSolve_PDIPM(Tao tao)
857: {
858:   PetscErrorCode     ierr;
859:   TAO_PDIPM          *pdipm = (TAO_PDIPM*)tao->data;
860:   SNESLineSearch     linesearch; /* SNESLineSearch context */
861:   Vec                dummy;

864:   if (!tao->constraints_equality && !tao->constraints_inequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_NULL,"Equality and inequality constraints are not set. Either set them or switch to a different algorithm");

866:   /* Initialize all variables */
867:   TaoPDIPMInitializeSolution(tao);

869:   /* Set linesearch */
870:   SNESGetLineSearch(pdipm->snes,&linesearch);
871:   SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
872:   SNESLineSearchShellSetUserFunc(linesearch,SNESLineSearch_PDIPM,tao);
873:   SNESLineSearchSetFromOptions(linesearch);

875:   tao->reason = TAO_CONTINUE_ITERATING;

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

881:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
882:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
883:   VecDestroy(&dummy);
884:   (*tao->ops->convergencetest)(tao,tao->cnvP);
885:   if (tao->reason) {
886:     SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
887:   }

889:   while (tao->reason == TAO_CONTINUE_ITERATING) {
890:     SNESConvergedReason reason;
891:     SNESSolve(pdipm->snes,NULL,pdipm->X);

893:     /* Check SNES convergence */
894:     SNESGetConvergedReason(pdipm->snes,&reason);
895:     if (reason < 0) {
896:       PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
897:     }

899:     /* Check TAO convergence */
900:     if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
901:   }
902:   return(0);
903: }

905: /*
906:   TaoView_PDIPM - View PDIPM

908:    Input Parameter:
909:     tao - TAO object
910:     viewer - PetscViewer

912:    Output:
913: */
914: PetscErrorCode TaoView_PDIPM(Tao tao,PetscViewer viewer)
915: {
916:   TAO_PDIPM      *pdipm = (TAO_PDIPM *)tao->data;

920:   tao->constrained = PETSC_TRUE;
921:   PetscViewerASCIIPushTab(viewer);
922:   PetscViewerASCIIPrintf(viewer,"Number of prime=%D, Number of dual=%D\n",pdipm->Nx+pdipm->Nci,pdipm->Nce + pdipm->Nci);
923:   if (pdipm->kkt_pd) {
924:     PetscViewerASCIIPrintf(viewer,"KKT shifts deltaw=%g, deltac=%g\n",(double)pdipm->deltaw,(double)pdipm->deltac);
925:   }
926:   PetscViewerASCIIPopTab(viewer);
927:   return(0);
928: }

930: /*
931:    TaoSetup_PDIPM - Sets up tao and pdipm

933:    Input Parameter:
934:    tao - TAO object

936:    Output:   pdipm - initialized object
937: */
938: PetscErrorCode TaoSetup_PDIPM(Tao tao)
939: {
940:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
941:   PetscErrorCode    ierr;
942:   MPI_Comm          comm;
943:   PetscMPIInt       rank,size;
944:   PetscInt          row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
945:   PetscInt          offset,*xa,*xb,i,j,rstart,rend;
946:   PetscScalar       one=1.0,neg_one=-1.0;
947:   const PetscInt    *cols,*rranges,*cranges,*aj,*ranges;
948:   const PetscScalar *aa,*Xarr;
949:   Mat               J,jac_equality_trans,jac_inequality_trans;
950:   Mat               Jce_xfixed_trans,Jci_xb_trans;
951:   PetscInt          *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];

954:   PetscObjectGetComm((PetscObject)tao,&comm);
955:   MPI_Comm_rank(comm,&rank);
956:   MPI_Comm_size(comm,&size);

958:   /* (1) Setup Bounds and create Tao vectors */
959:   TaoPDIPMSetUpBounds(tao);

961:   if (!tao->gradient) {
962:     VecDuplicate(tao->solution,&tao->gradient);
963:     VecDuplicate(tao->solution,&tao->stepdirection);
964:   }

966:   /* (2) Get sizes */
967:   /* Size of vector x - This is set by TaoSetInitialVector */
968:   VecGetSize(tao->solution,&pdipm->Nx);
969:   VecGetLocalSize(tao->solution,&pdipm->nx);

971:   /* Size of equality constraints and vectors */
972:   if (tao->constraints_equality) {
973:     VecGetSize(tao->constraints_equality,&pdipm->Ng);
974:     VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
975:   } else {
976:     pdipm->ng = pdipm->Ng = 0;
977:   }

979:   pdipm->nce = pdipm->ng + pdipm->nxfixed;
980:   pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;

982:   /* Size of inequality constraints and vectors */
983:   if (tao->constraints_inequality) {
984:     VecGetSize(tao->constraints_inequality,&pdipm->Nh);
985:     VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
986:   } else {
987:     pdipm->nh = pdipm->Nh = 0;
988:   }

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

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

997:   /* (3) Offsets for subvectors */
998:   pdipm->off_lambdae = pdipm->nx;
999:   pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
1000:   pdipm->off_z       = pdipm->off_lambdai + pdipm->nci;

1002:   /* (4) Create vectors and subvectors */
1003:   /* Ce and Ci vectors */
1004:   VecCreate(comm,&pdipm->ce);
1005:   VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
1006:   VecSetFromOptions(pdipm->ce);

1008:   VecCreate(comm,&pdipm->ci);
1009:   VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
1010:   VecSetFromOptions(pdipm->ci);

1012:   /* X=[x; lambdae; lambdai; z] for the big KKT system */
1013:   VecCreate(comm,&pdipm->X);
1014:   VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
1015:   VecSetFromOptions(pdipm->X);

1017:   /* Subvectors; they share local arrays with X */
1018:   VecGetArrayRead(pdipm->X,&Xarr);
1019:   /* x shares local array with X.x */
1020:   if (pdipm->Nx) {
1021:     VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
1022:   }

1024:   /* lambdae shares local array with X.lambdae */
1025:   if (pdipm->Nce) {
1026:     VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
1027:   }

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

1033:     VecCreate(comm,&pdipm->lambdae_xfixed);
1034:     VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
1035:     VecSetFromOptions(pdipm->lambdae_xfixed);
1036:   }

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

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

1046:   /* tao->DI which shares local array with X.lambdai_h */
1047:   if (pdipm->Nh) {
1048:     VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
1049:   }
1050:   VecCreate(comm,&pdipm->lambdai_xb);
1051:   VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
1052:   VecSetFromOptions(pdipm->lambdai_xb);

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

1056:   /* (5) Create Jacobians Jce_xfixed and Jci */
1057:   /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
1058:   if (pdipm->Nxfixed) {
1059:     /* Create Jce_xfixed */
1060:     MatCreate(comm,&pdipm->Jce_xfixed);
1061:     MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1062:     MatSetFromOptions(pdipm->Jce_xfixed);
1063:     MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
1064:     MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);

1066:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
1067:     ISGetIndices(pdipm->isxfixed,&cols);
1068:     k = 0;
1069:     for (row = Jcrstart; row < Jcrend; row++) {
1070:       MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
1071:       k++;
1072:     }
1073:     ISRestoreIndices(pdipm->isxfixed, &cols);
1074:     MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1075:     MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1076:   }

1078:   /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
1079:   MatCreate(comm,&pdipm->Jci_xb);
1080:   MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1081:   MatSetFromOptions(pdipm->Jci_xb);
1082:   MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
1083:   MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);

1085:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
1086:   offset = Jcrstart;
1087:   if (pdipm->Nxub) {
1088:     /* Add xub to Jci_xb */
1089:     ISGetIndices(pdipm->isxub,&cols);
1090:     k = 0;
1091:     for (row = offset; row < offset + pdipm->nxub; row++) {
1092:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1093:       k++;
1094:     }
1095:     ISRestoreIndices(pdipm->isxub, &cols);
1096:   }

1098:   if (pdipm->Nxlb) {
1099:     /* Add xlb to Jci_xb */
1100:     ISGetIndices(pdipm->isxlb,&cols);
1101:     k = 0;
1102:     offset += pdipm->nxub;
1103:     for (row = offset; row < offset + pdipm->nxlb; row++) {
1104:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
1105:       k++;
1106:     }
1107:     ISRestoreIndices(pdipm->isxlb, &cols);
1108:   }

1110:   /* Add xbox to Jci_xb */
1111:   if (pdipm->Nxbox) {
1112:     ISGetIndices(pdipm->isxbox,&cols);
1113:     k = 0;
1114:     offset += pdipm->nxlb;
1115:     for (row = offset; row < offset + pdipm->nxbox; row++) {
1116:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1117:       tmp = row + pdipm->nxbox;
1118:       MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
1119:       k++;
1120:     }
1121:     ISRestoreIndices(pdipm->isxbox, &cols);
1122:   }

1124:   MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1125:   MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1126:   /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */

1128:   /* (6) Set up ISs for PC Fieldsplit */
1129:   if (pdipm->solve_reduced_kkt) {
1130:     PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
1131:     for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1132:     for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;

1134:     ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
1135:     ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
1136:   }

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

1141:   /* Get rstart of KKT matrix */
1142:   MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
1143:   rstart -= pdipm->n;

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

1147:   PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
1148:   MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
1149:   MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
1150:   MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);

1152:   MatGetOwnershipRanges(tao->hessian,&rranges);
1153:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

1155:   if (pdipm->Ng) {
1156:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
1157:     MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
1158:   }
1159:   if (pdipm->Nh) {
1160:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
1161:     MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
1162:   }

1164:   /* Count dnz,onz for preallocation of KKT matrix */
1165:   jac_equality_trans   = pdipm->jac_equality_trans;
1166:   jac_inequality_trans = pdipm->jac_inequality_trans;
1167:   nce_all = pdipm->nce_all;

1169:   if (pdipm->Nxfixed) {
1170:     MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
1171:   }
1172:   MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);

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

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

1180:   /* Insert tao->hessian */
1181:   MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
1182:   for (i=0; i<pdipm->nx; i++){
1183:     row = rstart + i;

1185:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1186:     proc = 0;
1187:     for (j=0; j < nc; j++) {
1188:       while (aj[j] >= cranges[proc+1]) proc++;
1189:       col = aj[j] - cranges[proc] + Jranges[proc];
1190:       MatPreallocateSet(row,1,&col,dnz,onz);
1191:     }
1192:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);

1194:     if (pdipm->ng) {
1195:       /* Insert grad g' */
1196:       MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1197:       MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
1198:       proc = 0;
1199:       for (j=0; j < nc; j++) {
1200:         /* find row ownership of */
1201:         while (aj[j] >= ranges[proc+1]) proc++;
1202:         nx_all = rranges[proc+1] - rranges[proc];
1203:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1204:         MatPreallocateSet(row,1,&col,dnz,onz);
1205:       }
1206:       MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1207:     }

1209:     /* Insert Jce_xfixed^T' */
1210:     if (pdipm->nxfixed) {
1211:       MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1212:       MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1213:       proc = 0;
1214:       for (j=0; j < nc; j++) {
1215:         /* find row ownership of */
1216:         while (aj[j] >= ranges[proc+1]) proc++;
1217:         nx_all = rranges[proc+1] - rranges[proc];
1218:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1219:         MatPreallocateSet(row,1,&col,dnz,onz);
1220:       }
1221:       MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1222:     }

1224:     if (pdipm->nh) {
1225:       /* Insert -grad h' */
1226:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1227:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1228:       proc = 0;
1229:       for (j=0; j < nc; j++) {
1230:         /* find row ownership of */
1231:         while (aj[j] >= ranges[proc+1]) proc++;
1232:         nx_all = rranges[proc+1] - rranges[proc];
1233:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1234:         MatPreallocateSet(row,1,&col,dnz,onz);
1235:       }
1236:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1237:     }

1239:     /* Insert Jci_xb^T' */
1240:     MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1241:     MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1242:     proc = 0;
1243:     for (j=0; j < nc; j++) {
1244:       /* find row ownership of */
1245:       while (aj[j] >= ranges[proc+1]) proc++;
1246:       nx_all = rranges[proc+1] - rranges[proc];
1247:       col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1248:       MatPreallocateSet(row,1,&col,dnz,onz);
1249:     }
1250:     MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1251:   }

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

1259:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1260:       proc = 0;
1261:       for (j=0; j < nc; j++) {
1262:         while (aj[j] >= cranges[proc+1]) proc++;
1263:         col = aj[j] - cranges[proc] + Jranges[proc];
1264:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1265:       }
1266:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1267:     }
1268:   }
1269:   /* Jce_xfixed */
1270:   if (pdipm->Nxfixed) {
1271:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1272:     for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1273:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1275:       MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1276:       if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");

1278:       proc = 0;
1279:       j    = 0;
1280:       while (cols[j] >= cranges[proc+1]) proc++;
1281:       col = cols[j] - cranges[proc] + Jranges[proc];
1282:       MatPreallocateSet(row,1,&col,dnz,onz);
1283:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1284:     }
1285:   }

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

1293:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1294:       proc = 0;
1295:       for (j=0; j < nc; j++) {
1296:         while (aj[j] >= cranges[proc+1]) proc++;
1297:         col = aj[j] - cranges[proc] + Jranges[proc];
1298:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1299:       }
1300:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1301:     }
1302:     /* I */
1303:     for (i=0; i < pdipm->nh; i++){
1304:       row = rstart + pdipm->off_lambdai + i;
1305:       col = rstart + pdipm->off_z + i;
1306:       MatPreallocateSet(row,1,&col,dnz,onz);
1307:     }
1308:   }

1310:   /* Jci_xb */
1311:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1312:   for (i=0; i < (pdipm->nci - pdipm->nh); i++){
1313:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1315:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1316:     if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1317:     proc = 0;
1318:     for (j=0; j < nc; j++) {
1319:       while (cols[j] >= cranges[proc+1]) proc++;
1320:       col = cols[j] - cranges[proc] + Jranges[proc];
1321:       MatPreallocateSet(row,1,&col,dnz,onz);
1322:     }
1323:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1324:     /* I */
1325:     col = rstart + pdipm->off_z + pdipm->nh + i;
1326:     MatPreallocateSet(row,1,&col,dnz,onz);
1327:   }

1329:   /* 4-th Row block of KKT matrix: Z and Ci */
1330:   for (i=0; i < pdipm->nci; i++) {
1331:     row     = rstart + pdipm->off_z + i;
1332:     cols1[0] = rstart + pdipm->off_lambdai + i;
1333:     cols1[1] = row;
1334:     MatPreallocateSet(row,2,cols1,dnz,onz);
1335:   }

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

1340:   /* Create KKT matrix */
1341:   MatCreate(comm,&J);
1342:   MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1343:   MatSetFromOptions(J);
1344:   MatSeqAIJSetPreallocation(J,0,dnz);
1345:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1346:   /* MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); */
1347:   MatPreallocateFinalize(dnz,onz);
1348:   pdipm->K = J;

1350:   /* (8) Set up nonlinear solver SNES */
1351:   SNESSetFunction(pdipm->snes,NULL,TaoSNESFunction_PDIPM,(void*)tao);
1352:   SNESSetJacobian(pdipm->snes,J,J,TaoSNESJacobian_PDIPM,(void*)tao);

1354:   if (pdipm->solve_reduced_kkt) {
1355:     PC pc;
1356:     KSPGetPC(tao->ksp,&pc);
1357:     PCSetType(pc,PCFIELDSPLIT);
1358:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1359:     PCFieldSplitSetIS(pc,"2",pdipm->is2);
1360:     PCFieldSplitSetIS(pc,"1",pdipm->is1);
1361:   }
1362:   SNESSetFromOptions(pdipm->snes);

1364:   /* (9) Insert constant entries to  K */
1365:   /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1366:   MatGetOwnershipRange(J,&rstart,&rend);
1367:   for (i=rstart; i<rend; i++){
1368:     MatSetValue(J,i,i,0.0,INSERT_VALUES);
1369:   }
1370:   /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1371:   if (pdipm->kkt_pd){
1372:       for (i=0; i<pdipm->nh; i++){
1373:         row  = rstart + i;
1374:         MatSetValue(J,row,row,pdipm->deltaw,INSERT_VALUES);
1375:       }
1376:   }

1378:   /* Row block of K: [ grad Ce, 0, 0, 0] */
1379:   if (pdipm->Nxfixed) {
1380:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1381:     for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1382:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1384:       MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1385:       proc = 0;
1386:       for (j=0; j < nc; j++) {
1387:         while (cols[j] >= cranges[proc+1]) proc++;
1388:         col = cols[j] - cranges[proc] + Jranges[proc];
1389:         MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1390:         MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1391:       }
1392:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1393:     }
1394:   }

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

1401:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1402:     proc = 0;
1403:     for (j=0; j < nc; j++) {
1404:       while (cols[j] >= cranges[proc+1]) proc++;
1405:       col = cols[j] - cranges[proc] + Jranges[proc];
1406:       MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1407:       MatSetValue(J,row,col,-aa[j],INSERT_VALUES);
1408:     }
1409:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);

1411:     col = rstart + pdipm->off_z + pdipm->nh + i;
1412:     MatSetValue(J,row,col,1,INSERT_VALUES);
1413:   }

1415:   for (i=0; i < pdipm->nh; i++){
1416:     row = rstart + pdipm->off_lambdai + i;
1417:     col = rstart + pdipm->off_z + i;
1418:     MatSetValue(J,row,col,1,INSERT_VALUES);
1419:   }

1421:   /* Row block of K: [ 0, 0, I, ...] */
1422:   for (i=0; i < pdipm->nci; i++){
1423:     row = rstart + pdipm->off_z + i;
1424:     col = rstart + pdipm->off_lambdai + i;
1425:     MatSetValue(J,row,col,1,INSERT_VALUES);
1426:   }

1428:   if (pdipm->Nxfixed) {
1429:     MatDestroy(&Jce_xfixed_trans);
1430:   }
1431:   MatDestroy(&Jci_xb_trans);
1432:   PetscFree3(ng_all,nh_all,Jranges);

1434:   /* (10) Setup PCPreSolve() for pdipm->solve_symmetric_kkt */
1435:   if (pdipm->solve_symmetric_kkt) {
1436:     KSP ksp;
1437:     PC  pc;
1438:     SNESGetKSP(pdipm->snes,&ksp);
1439:     KSPGetPC(ksp,&pc);
1440:     pc->ops->presolve = PCPreSolve_PDIPM;
1441:   }
1442:   return(0);
1443: }

1445: /*
1446:    TaoDestroy_PDIPM - Destroys the pdipm object

1448:    Input:
1449:    full pdipm

1451:    Output:
1452:    Destroyed pdipm
1453: */
1454: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1455: {
1456:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1460:   /* Freeing Vectors assocaiated with KKT (X) */
1461:   VecDestroy(&pdipm->x); /* Solution x */
1462:   VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1463:   VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1464:   VecDestroy(&pdipm->z);       /* Slack variables */
1465:   VecDestroy(&pdipm->X);       /* Big KKT system vector [x; lambdae; lambdai; z] */

1467:   /* work vectors */
1468:   VecDestroy(&pdipm->lambdae_xfixed);
1469:   VecDestroy(&pdipm->lambdai_xb);

1471:   /* Legrangian equality and inequality Vec */
1472:   VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1473:   VecDestroy(&pdipm->ci); /* Vec of inequality constraints */

1475:   /* Matrices */
1476:   MatDestroy(&pdipm->Jce_xfixed);
1477:   MatDestroy(&pdipm->Jci_xb); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1478:   MatDestroy(&pdipm->K);

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

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

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

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

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

1501:   if (pdipm->solve_reduced_kkt) {
1502:     ISDestroy(&pdipm->is1);
1503:     ISDestroy(&pdipm->is2);
1504:   }

1506:   /* SNES */
1507:   SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1508:   PetscFree(pdipm->nce_all);
1509:   MatDestroy(&pdipm->jac_equality_trans);
1510:   MatDestroy(&pdipm->jac_inequality_trans);

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

1515:   /* Destroy Dual */
1516:   VecDestroy(&tao->DE); /* equality dual */
1517:   VecDestroy(&tao->DI); /* dinequality dual */
1518:   return(0);
1519: }

1521: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1522: {
1523:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1527:   PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1528:   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);
1529:   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);
1530:   PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1531:   PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1532:   PetscOptionsBool("-tao_pdipm_symmetric_kkt","Solve non reduced symmetric KKT system",NULL,pdipm->solve_symmetric_kkt,&pdipm->solve_symmetric_kkt,NULL);
1533:   PetscOptionsBool("-tao_pdipm_kkt_shift_pd","Add shifts to make KKT matrix positive definite",NULL,pdipm->kkt_pd,&pdipm->kkt_pd,NULL);
1534:   PetscOptionsTail();
1535:   return(0);
1536: }

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

1541:   Option Database Keys:
1542: +   -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1543: .   -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1544: .   -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1545: .   -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1546: -   -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite

1548:   Level: beginner
1549: M*/
1550: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1551: {
1552:   TAO_PDIPM      *pdipm;

1556:   tao->ops->setup          = TaoSetup_PDIPM;
1557:   tao->ops->solve          = TaoSolve_PDIPM;
1558:   tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1559:   tao->ops->view           = TaoView_PDIPM;
1560:   tao->ops->destroy        = TaoDestroy_PDIPM;

1562:   PetscNewLog(tao,&pdipm);
1563:   tao->data = (void*)pdipm;

1565:   pdipm->nx      = pdipm->Nx      = 0;
1566:   pdipm->nxfixed = pdipm->Nxfixed = 0;
1567:   pdipm->nxlb    = pdipm->Nxlb    = 0;
1568:   pdipm->nxub    = pdipm->Nxub    = 0;
1569:   pdipm->nxbox   = pdipm->Nxbox   = 0;
1570:   pdipm->nxfree  = pdipm->Nxfree  = 0;

1572:   pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1573:   pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1574:   pdipm->n  = pdipm->N  = 0;
1575:   pdipm->mu = 1.0;
1576:   pdipm->mu_update_factor = 0.1;

1578:   pdipm->deltaw     = 0.0;
1579:   pdipm->lastdeltaw = 3*1.e-4;
1580:   pdipm->deltac     = 0.0;
1581:   pdipm->kkt_pd     = PETSC_FALSE;

1583:   pdipm->push_init_slack     = 1.0;
1584:   pdipm->push_init_lambdai   = 1.0;
1585:   pdipm->solve_reduced_kkt   = PETSC_FALSE;
1586:   pdipm->solve_symmetric_kkt = PETSC_TRUE;

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

1592:   SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1593:   SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1594:   SNESGetKSP(pdipm->snes,&tao->ksp);
1595:   PetscObjectReference((PetscObject)tao->ksp);
1596:   KSPSetApplicationContext(tao->ksp,(void *)tao);
1597:   return(0);
1598: }