Actual source code: pdipm.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/constrained/impls/ipm/pdipm.h>
  3: #include <petscsnes.h>

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

  8:    Collective on tao

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

 14:    Level: beginner

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

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

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

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

 41: /*
 42:   TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x

 44:   Collective

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

 50:    Level: beginner

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

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

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

 70:   /* (1) Update ce vector */
 71:   VecGetArray(pdipm->ce,&carr);

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

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

 91:   /* (2) Update ci vector */
 92:   VecGetArray(pdipm->ci,&carr);

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

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

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

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

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

141: /*
142:    TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x

144:    Collective

146:    Input Parameter:
147: .  tao - holds pdipm and XL & XU

149:    Level: beginner

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

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

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

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

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

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

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

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

229:    Collective

231:    Input Parameter:
232: .  tao - Tao context

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

245:   VecGetArray(pdipm->X,&Xarr);

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

252:   /* Initialize X.lambdae = 0.0 */
253:   if (pdipm->lambdae) {
254:     VecSet(pdipm->lambdae,0.0);
255:   }
256:   /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
257:   if (pdipm->lambdai) {
258:     VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
259:   }
260:   if (pdipm->z) {
261:     VecSet(pdipm->z,pdipm->push_init_slack);
262:   }

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

286:   VecRestoreArray(pdipm->X,&Xarr);
287:   return(0);
288: }

290: /*
291:    TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X

293:    Input Parameter:
294:    snes - SNES context
295:    X - KKT Vector
296:    *ctx - pdipm context

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

317:   PetscObjectGetComm((PetscObject)snes,&comm);
318:   MPI_Comm_rank(comm,&rank);
319:   MPI_Comm_rank(comm,&size);

321:   MatGetOwnershipRanges(Jpre,&Jranges);
322:   MatGetOwnershipRange(Jpre,&Jrstart,NULL);
323:   MatGetOwnershipRangesColumn(tao->hessian,&rranges);
324:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

326:   VecGetArrayRead(X,&Xarr);

328:   if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
329:     vals[0] = 1.0;
330:     for (i=0; i < pdipm->nci; i++) {
331:         row     = Jrstart + pdipm->off_z + i;
332:         cols[0] = Jrstart + pdipm->off_lambdai + i;
333:         cols[1] = row;
334:         vals[1] = Xarr[pdipm->off_lambdai + i]/Xarr[pdipm->off_z + i];
335:         MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
336:     }
337:   } else {
338:     /* (2) insert Z and Ci to Jpre -- overwrite existing values */
339:     for (i=0; i < pdipm->nci; i++) {
340:       row     = Jrstart + pdipm->off_z + i;
341:       cols[0] = Jrstart + pdipm->off_lambdai + i;
342:       cols[1] = row;
343:       vals[0] = Xarr[pdipm->off_z + i];
344:       vals[1] = Xarr[pdipm->off_lambdai + i];
345:       MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
346:     }
347:   }

349:   /* (3) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
350:   if (pdipm->Ng) {
351:     MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
352:     for (i=0; i<pdipm->ng; i++){
353:       row = Jrstart + pdipm->off_lambdae + i;

355:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
356:       proc = 0;
357:       for (j=0; j < nc; j++) {
358:         while (aj[j] >= cranges[proc+1]) proc++;
359:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
360:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
361:       }
362:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
363:       if (pdipm->kkt_pd) {
364:         /* (3) insert 2nd row block of Jpre: [ grad g, \delta_c*I, 0, 0] */
365:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
366:       }
367:     }
368:   }

370:   if (pdipm->Nh) {
371:     /* (4) insert 3nd row block of Jpre: [ -grad h, 0, deltac, I] */
372:     MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
373:     for (i=0; i < pdipm->nh; i++){
374:       row = Jrstart + pdipm->off_lambdai + i;
375:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
376:       proc = 0;
377:       for (j=0; j < nc; j++) {
378:         while (aj[j] >= cranges[proc+1]) proc++;
379:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
380:         MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
381:       }
382:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
383:       if (pdipm->kkt_pd) {
384:         MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
385:       }
386:     }
387:   }

389:   /* (5) insert Wxx, grad g' and -grad h' to Jpre */
390:   if (pdipm->Ng) {
391:     MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
392:   }
393:   if (pdipm->Nh) {
394:     MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
395:   }

397:   VecPlaceArray(pdipm->x,Xarr);
398:   TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
399:   VecResetArray(pdipm->x);

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

405:     /* insert Wxx */
406:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
407:     proc = 0;
408:     for (j=0; j < nc; j++) {
409:       while (aj[j] >= cranges[proc+1]) proc++;
410:       cols[0] = aj[j] - cranges[proc] + Jranges[proc];
411:       if (row == cols[0] && pdipm->kkt_pd) {
412:         /* add diagonal shift to Wxx component */
413:         MatSetValue(Jpre,row,cols[0],aa[j]+pdipm->deltaw,INSERT_VALUES);
414:       } else {
415:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
416:       }
417:     }
418:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);

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

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

452:   /* (6) assemble Jpre and J */
453:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
454:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);

456:   if (J != Jpre) {
457:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
458:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
459:   }
460:   return(0);
461: }

463: /*
464:    TaoSnesFunction_PDIPM - Evaluate KKT function at X

466:    Input Parameter:
467:    snes - SNES context
468:    X - KKT Vector
469:    *ctx - pdipm

471:    Output Parameter:
472:    F - Updated Lagrangian vector
473: */
474: PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
475: {
476:   PetscErrorCode    ierr;
477:   Tao               tao=(Tao)ctx;
478:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
479:   PetscScalar       *Farr,*tmparr;
480:   Vec               x,L1;
481:   PetscInt          i;
482:   PetscReal         res[2],cnorm[2];
483:   const PetscScalar *Xarr,*carr,*zarr,*larr;

486:   VecSet(F,0.0);

488:   VecGetArrayRead(X,&Xarr);
489:   VecGetArray(F,&Farr);

491:   /* (0) Evaluate f, fx, Gx, Hx at X.x Note: pdipm->x is not changed below */
492:   x = pdipm->x;
493:   VecPlaceArray(x,Xarr);
494:   TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);

496:   /* Update ce, ci, and Jci at X.x */
497:   TaoPDIPMUpdateConstraints(tao,x);
498:   VecResetArray(x);

500:   /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
501:   L1 = pdipm->x;
502:   VecPlaceArray(L1,Farr);
503:   if (pdipm->Nci) {
504:     if (pdipm->Nh) {
505:       /* L1 += gradH'*DI. Note: tao->DI is not changed below */
506:       VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
507:       MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
508:       VecResetArray(tao->DI);
509:     }

511:     /* L1 += Jci_xb'*lambdai_xb */
512:     VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
513:     MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
514:     VecResetArray(pdipm->lambdai_xb);

516:     /* (1.4) L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
517:     VecScale(L1,-1.0);
518:   }

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

523:   if (pdipm->Nce) {
524:     if (pdipm->Ng) {
525:       /* L1 += gradG'*DE. Note: tao->DE is not changed below */
526:       VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
527:       MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
528:       VecResetArray(tao->DE);
529:     }
530:     if (pdipm->Nxfixed) {
531:       /* L1 += Jce_xfixed'*lambdae_xfixed */
532:       VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
533:       MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
534:       VecResetArray(pdipm->lambdae_xfixed);
535:     }
536:   }
537:   VecNorm(L1,NORM_2,&res[0]);
538:   VecResetArray(L1);

540:   /* (2) L2 = ce(x) */
541:   if (pdipm->Nce) {
542:     VecGetArrayRead(pdipm->ce,&carr);
543:     for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
544:     VecRestoreArrayRead(pdipm->ce,&carr);
545:   }
546:   VecNorm(pdipm->ce,NORM_2,&cnorm[0]);

548:   if (pdipm->Nci) {
549:     if (pdipm->solve_symmetric_kkt) {
550:       /* (3) L3 = ci(x) - z;
551:          (4) L4 = Lambdai * e - mu/z *e
552:       */
553:       VecGetArrayRead(pdipm->ci,&carr);
554:       larr = Xarr+pdipm->off_lambdai;
555:       zarr = Xarr+pdipm->off_z;
556:       for (i=0; i<pdipm->nci; i++) {
557:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
558:         Farr[pdipm->off_z       + i] = larr[i] - pdipm->mu/zarr[i];
559:       }
560:       VecRestoreArrayRead(pdipm->ci,&carr);
561:     } else {
562:       /* (3) L3 = ci(x) - z;
563:          (4) L4 = Z * Lambdai * e - mu * e
564:       */
565:       VecGetArrayRead(pdipm->ci,&carr);
566:       larr = Xarr+pdipm->off_lambdai;
567:       zarr = Xarr+pdipm->off_z;
568:       for (i=0; i<pdipm->nci; i++) {
569:         Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
570:         Farr[pdipm->off_z       + i] = zarr[i]*larr[i] - pdipm->mu;
571:       }
572:       VecRestoreArrayRead(pdipm->ci,&carr);
573:     }
574:   }

576:   VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
577:   VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
578:   VecResetArray(pdipm->ci);

580:   /* note: pdipm->z is not changed below */
581:   if (pdipm->z) {
582:     if (pdipm->solve_symmetric_kkt) {
583:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
584:       if (pdipm->Nci) {
585:         zarr = Xarr+pdipm->off_z;
586:         VecGetArrayWrite(pdipm->z,&tmparr);
587:         for (i=0; i<pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
588:         VecRestoreArrayWrite(pdipm->z,&tmparr);
589:       }

591:       VecNorm(pdipm->z,NORM_2,&res[1]);
592:       if (pdipm->Nci) {
593:         zarr = Xarr+pdipm->off_z;
594:         VecGetArray(pdipm->z,&tmparr);
595:         for (i=0; i<pdipm->nci; i++) {
596:           tmparr[i] /= Xarr[pdipm->off_z + i];
597:         }
598:         VecRestoreArray(pdipm->z,&tmparr);
599:       }
600:       VecResetArray(pdipm->z);
601:     } else {
602:       VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
603:       VecNorm(pdipm->z,NORM_2,&res[1]);
604:       VecResetArray(pdipm->z);
605:     }
606:   } else res[1] = 0.0;

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

611:   VecRestoreArrayRead(X,&Xarr);
612:   VecRestoreArray(F,&Farr);
613:   return(0);
614: }

616: #include <petsc/private/matimpl.h>
617: /*
618:    PDIPMLineSearch - Custom line search used with PDIPM.

620:    Collective on TAO

622:    Notes:
623:    PDIPMLineSearch employs a simple backtracking line-search to keep
624:    the slack variables (z) and inequality constraints lagrange multipliers
625:    (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
626:    alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
627:    slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
628:    are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
629:    is also updated as mu = mu + z'lambdai/Nci
630: */
631: PetscErrorCode PDIPMLineSearch(SNESLineSearch linesearch,void *ctx)
632: {
633:   PetscErrorCode    ierr;
634:   Tao               tao=(Tao)ctx;
635:   TAO_PDIPM         *pdipm = (TAO_PDIPM*)tao->data;
636:   SNES              snes;
637:   KSP               ksp;
638:   PC                pc;
639:   PCType            ptype;
640:   Mat               Factor;
641:   Vec               X,F,Y,W,G;
642:   PetscInt          i,iter,nneg,nzero,npos;
643:   PetscReal         alpha_p=1.0,alpha_d=1.0,alpha[4];
644:   PetscScalar       *Xarr,*z,*lambdai,dot,*taosolarr;
645:   const PetscScalar *dXarr,*dz,*dlambdai;
646:   PetscBool         isCHOL;

649:   SNESLineSearchGetSNES(linesearch,&snes);
650:   SNESGetIterationNumber(snes,&iter);

652:   SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
653:   SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);

655:   VecGetArray(X,&Xarr);
656:   VecGetArrayRead(Y,&dXarr);
657:   z  = Xarr + pdipm->off_z;
658:   dz = dXarr + pdipm->off_z;
659:   for (i=0; i < pdipm->nci; i++) {
660:     if (z[i] - dz[i] < 0.0) {
661:       alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
662:     }
663:   }

665:   lambdai  = Xarr + pdipm->off_lambdai;
666:   dlambdai = dXarr + pdipm->off_lambdai;

668:   for (i=0; i<pdipm->nci; i++) {
669:     if (lambdai[i] - dlambdai[i] < 0.0) {
670:       alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
671:     }
672:   }

674:   alpha[0] = alpha_p;
675:   alpha[1] = alpha_d;
676:   VecRestoreArrayRead(Y,&dXarr);
677:   VecRestoreArray(X,&Xarr);

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

682:   alpha_p = alpha[2];
683:   alpha_d = alpha[3];

685:   VecGetArray(X,&Xarr);
686:   VecGetArrayRead(Y,&dXarr);
687:   for (i=0; i<pdipm->nx; i++) {
688:     Xarr[i] = Xarr[i] - alpha_p * dXarr[i];
689:   }

691:   for (i=0; i<pdipm->nce; i++) {
692:     Xarr[i+pdipm->off_lambdae] = Xarr[i+pdipm->off_lambdae] - alpha_d * dXarr[i+pdipm->off_lambdae];
693:   }

695:   for (i=0; i<pdipm->nci; i++) {
696:     Xarr[i+pdipm->off_lambdai] = Xarr[i+pdipm->off_lambdai] - alpha_d * dXarr[i+pdipm->off_lambdai];
697:   }

699:   for (i=0; i<pdipm->nci; i++) {
700:     Xarr[i+pdipm->off_z] = Xarr[i+pdipm->off_z] - alpha_p * dXarr[i+pdipm->off_z];
701:   }

703:   VecGetArray(tao->solution,&taosolarr);
704:   PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
705:   VecRestoreArray(tao->solution,&taosolarr);


708:   VecRestoreArray(X,&Xarr);
709:   VecRestoreArrayRead(Y,&dXarr);

711:   /* Evaluate F at X */
712:   SNESComputeFunction(snes,X,F);
713:   SNESLineSearchComputeNorms(linesearch); /* must call this func, do not know why */

715:   /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
716:   if (pdipm->z) {
717:     VecDot(pdipm->z,pdipm->lambdai,&dot);
718:   } else dot = 0.0;

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

723:   /* Update F; get tao->residual and tao->cnorm */
724:   TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);

726:   tao->niter++;
727:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
728:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);

730:   (*tao->ops->convergencetest)(tao,tao->cnvP);
731:   if (tao->reason) {
732:     SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
733:   }

735:   if (!pdipm->kkt_pd) return(0);

737:   /* Get the inertia of Cholesky factor to set shifts for next SNES interation */
738:   SNESGetKSP(snes,&ksp);
739:   KSPGetPC(ksp,&pc);
740:   PCGetType(pc,&ptype);
741:   PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);

743:   if (isCHOL) {
744:     PetscMPIInt       size;
745:     PCFactorGetMatrix(pc,&Factor);
746:     MPI_Comm_size(PetscObjectComm((PetscObject)Factor),&size);
747:     if (Factor->ops->getinertia) {
748: #if defined(PETSC_HAVE_MUMPS)
749:       MatSolverType     stype;
750:       PetscBool         isMUMPS;
751:       PCFactorGetMatSolverType(pc,&stype);
752:       PetscStrcmp(stype, MATSOLVERMUMPS, &isMUMPS);
753:       if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
754:         MatMumpsSetIcntl(Factor,24,1);
755:         if (size > 1) {
756:           MatMumpsSetIcntl(Factor,13,1);
757:         }
758:       }
759: #else
760:       if (size > 1) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires external package MUMPS");
761: #endif
762:       MatGetInertia(Factor,&nneg,&nzero,&npos);

764:       if (npos < pdipm->Nx+pdipm->Nci) {
765:         pdipm->deltaw = PetscMax(pdipm->lastdeltaw/3, 1.e-4*PETSC_MACHINE_EPSILON);
766:         PetscInfo5(tao,"Test reduced deltaw=%g; previous MatInertia: nneg %d, nzero %d, npos %d(<%d)\n",pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
767:         TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
768:         PCSetUp(pc);
769:         MatGetInertia(Factor,&nneg,&nzero,&npos);

771:         if (npos < pdipm->Nx+pdipm->Nci) {
772:           pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
773:           while (npos < pdipm->Nx+pdipm->Nci && pdipm->deltaw <= 1.e10) { /* increase deltaw */
774:             PetscInfo5(tao,"  deltaw=%g fails, MatInertia: nneg %d, nzero %d, npos %d(<%d)\n",pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
775:             pdipm->deltaw = PetscMin(8*pdipm->deltaw,PetscPowReal(10,20));
776:             TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
777:             PCSetUp(pc);
778:             MatGetInertia(Factor,&nneg,&nzero,&npos);
779:           }

781:           if (pdipm->deltaw >= 1.e10) {
782:             SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_CONV_FAILED,"Reached maximum delta w will not converge, try different inital x0");
783:           }
784:           PetscInfo1(tao,"Updated deltaw %g\n",pdipm->deltaw);
785:           pdipm->lastdeltaw = pdipm->deltaw;
786:           pdipm->deltaw     = 0.0;
787:         }
788:       }

790:       if (nzero) { /* Jacobian is singular */
791:         if (pdipm->deltac == 0.0) {
792:           pdipm->deltac = 1.e8*PETSC_MACHINE_EPSILON;
793:         } else {
794:           pdipm->deltac = pdipm->deltac*PetscPowReal(pdipm->mu,.25);
795:         }
796:         PetscInfo4(tao,"Updated deltac=%g, MatInertia: nneg %D, nzero %D(!=0), npos %D\n",pdipm->deltac,nneg,nzero,npos);
797:       }
798:     } else
799:       SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires an external package that supports MatGetInertia()");
800:   }
801:   return(0);
802: }

804: /*
805:    TaoSolve_PDIPM

807:    Input Parameter:
808:    tao - TAO context

810:    Output Parameter:
811:    tao - TAO context
812: */
813: PetscErrorCode TaoSolve_PDIPM(Tao tao)
814: {
815:   PetscErrorCode     ierr;
816:   TAO_PDIPM          *pdipm = (TAO_PDIPM*)tao->data;
817:   SNESLineSearch     linesearch; /* SNESLineSearch context */
818:   Vec                dummy;

821:   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");

823:   /* Initialize all variables */
824:   TaoPDIPMInitializeSolution(tao);

826:   /* Set linesearch */
827:   SNESGetLineSearch(pdipm->snes,&linesearch);
828:   SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
829:   SNESLineSearchShellSetUserFunc(linesearch,PDIPMLineSearch,tao);
830:   SNESLineSearchSetFromOptions(linesearch);

832:   tao->reason = TAO_CONTINUE_ITERATING;

834:   /* -tao_monitor for iteration 0 and check convergence */
835:   VecDuplicate(pdipm->X,&dummy);
836:   TaoSNESFunction_PDIPM(pdipm->snes,pdipm->X,dummy,(void*)tao);

838:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
839:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
840:   VecDestroy(&dummy);
841:   (*tao->ops->convergencetest)(tao,tao->cnvP);
842:   if (tao->reason) {
843:     SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
844:   }

846:   while (tao->reason == TAO_CONTINUE_ITERATING) {
847:     SNESConvergedReason reason;
848:     SNESSolve(pdipm->snes,NULL,pdipm->X);

850:     /* Check SNES convergence */
851:     SNESGetConvergedReason(pdipm->snes,&reason);
852:     if (reason < 0) {
853:       PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
854:     }

856:     /* Check TAO convergence */
857:     if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
858:   }
859:   return(0);
860: }

862: /*
863:   TaoView_PDIPM - View PDIPM

865:    Input Parameter:
866:     tao - TAO object
867:     viewer - PetscViewer

869:    Output:
870: */
871: PetscErrorCode TaoView_PDIPM(Tao tao,PetscViewer viewer)
872: {
873:   TAO_PDIPM      *pdipm = (TAO_PDIPM *)tao->data;

877:   tao->constrained = PETSC_TRUE;
878:   PetscViewerASCIIPushTab(viewer);
879:   PetscViewerASCIIPrintf(viewer,"Number of prime=%D, Number of dual=%D\n",pdipm->Nx+pdipm->Nci,pdipm->Nce + pdipm->Nci);
880:   if (pdipm->kkt_pd) {
881:     PetscViewerASCIIPrintf(viewer,"KKT shifts deltaw=%g, deltac=%g\n",pdipm->deltaw,pdipm->deltac);
882:   }
883:   PetscViewerASCIIPopTab(viewer);
884:   return(0);
885: }

887: /*
888:    TaoSetup_PDIPM - Sets up tao and pdipm

890:    Input Parameter:
891:    tao - TAO object

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

911:   PetscObjectGetComm((PetscObject)tao,&comm);
912:   MPI_Comm_rank(comm,&rank);
913:   MPI_Comm_size(comm,&size);

915:   /* (1) Setup Bounds and create Tao vectors */
916:   TaoPDIPMSetUpBounds(tao);

918:   if (!tao->gradient) {
919:     VecDuplicate(tao->solution,&tao->gradient);
920:     VecDuplicate(tao->solution,&tao->stepdirection);
921:   }

923:   /* (2) Get sizes */
924:   /* Size of vector x - This is set by TaoSetInitialVector */
925:   VecGetSize(tao->solution,&pdipm->Nx);
926:   VecGetLocalSize(tao->solution,&pdipm->nx);

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

936:   pdipm->nce = pdipm->ng + pdipm->nxfixed;
937:   pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;

939:   /* Size of inequality constraints and vectors */
940:   if (tao->constraints_inequality) {
941:     VecGetSize(tao->constraints_inequality,&pdipm->Nh);
942:     VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
943:   } else {
944:     pdipm->nh = pdipm->Nh = 0;
945:   }

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

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

954:   /* (3) Offsets for subvectors */
955:   pdipm->off_lambdae = pdipm->nx;
956:   pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
957:   pdipm->off_z       = pdipm->off_lambdai + pdipm->nci;

959:   /* (4) Create vectors and subvectors */
960:   /* Ce and Ci vectors */
961:   VecCreate(comm,&pdipm->ce);
962:   VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
963:   VecSetFromOptions(pdipm->ce);

965:   VecCreate(comm,&pdipm->ci);
966:   VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
967:   VecSetFromOptions(pdipm->ci);

969:   /* X=[x; lambdae; lambdai; z] for the big KKT system */
970:   VecCreate(comm,&pdipm->X);
971:   VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
972:   VecSetFromOptions(pdipm->X);

974:   /* Subvectors; they share local arrays with X */
975:   VecGetArray(pdipm->X,&Xarr);
976:   /* x shares local array with X.x */
977:   if (pdipm->Nx) {
978:     VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
979:   }

981:   /* lambdae shares local array with X.lambdae */
982:   if (pdipm->Nce) {
983:     VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
984:   }

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

990:     VecCreate(comm,&pdipm->lambdae_xfixed);
991:     VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
992:     VecSetFromOptions(pdipm->lambdae_xfixed);
993:   }

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

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

1003:   /* tao->DI which shares local array with X.lambdai_h */
1004:   if (pdipm->Nh) {
1005:     VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
1006:   }

1008:   VecCreate(comm,&pdipm->lambdai_xb);
1009:   VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
1010:   VecSetFromOptions(pdipm->lambdai_xb);

1012:   VecRestoreArray(pdipm->X,&Xarr);

1014:   /* (5) Create Jacobians Jce_xfixed and Jci */
1015:   /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
1016:   if (pdipm->Nxfixed) {
1017:     /* Create Jce_xfixed */
1018:     MatCreate(comm,&pdipm->Jce_xfixed);
1019:     MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1020:     MatSetFromOptions(pdipm->Jce_xfixed);
1021:     MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
1022:     MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);

1024:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
1025:     ISGetIndices(pdipm->isxfixed,&cols);
1026:     k = 0;
1027:     for (row = Jcrstart; row < Jcrend; row++) {
1028:       MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
1029:       k++;
1030:     }
1031:     ISRestoreIndices(pdipm->isxfixed, &cols);
1032:     MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1033:     MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1034:   }

1036:   /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
1037:   MatCreate(comm,&pdipm->Jci_xb);
1038:   MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1039:   MatSetFromOptions(pdipm->Jci_xb);
1040:   MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
1041:   MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);

1043:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
1044:   offset = Jcrstart;
1045:   if (pdipm->Nxub) {
1046:     /* Add xub to Jci_xb */
1047:     ISGetIndices(pdipm->isxub,&cols);
1048:     k = 0;
1049:     for (row = offset; row < offset + pdipm->nxub; row++) {
1050:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1051:       k++;
1052:     }
1053:     ISRestoreIndices(pdipm->isxub, &cols);
1054:   }

1056:   if (pdipm->Nxlb) {
1057:     /* Add xlb to Jci_xb */
1058:     ISGetIndices(pdipm->isxlb,&cols);
1059:     k = 0;
1060:     offset += pdipm->nxub;
1061:     for (row = offset; row < offset + pdipm->nxlb; row++) {
1062:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
1063:       k++;
1064:     }
1065:     ISRestoreIndices(pdipm->isxlb, &cols);
1066:   }

1068:   /* Add xbox to Jci_xb */
1069:   if (pdipm->Nxbox) {
1070:     ISGetIndices(pdipm->isxbox,&cols);
1071:     k = 0;
1072:     offset += pdipm->nxlb;
1073:     for (row = offset; row < offset + pdipm->nxbox; row++) {
1074:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1075:       tmp = row + pdipm->nxbox;
1076:       MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
1077:       k++;
1078:     }
1079:     ISRestoreIndices(pdipm->isxbox, &cols);
1080:   }

1082:   MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1083:   MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1084:   /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */

1086:   /* (6) Set up ISs for PC Fieldsplit */
1087:   if (pdipm->solve_reduced_kkt) {
1088:     PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
1089:     for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1090:     for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;

1092:     ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
1093:     ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
1094:   }

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

1099:   /* Get rstart of KKT matrix */
1100:   MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
1101:   rstart -= pdipm->n;

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

1105:   PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
1106:   MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
1107:   MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
1108:   MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);

1110:   MatGetOwnershipRanges(tao->hessian,&rranges);
1111:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

1113:   if (pdipm->Ng) {
1114:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
1115:     MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
1116:   }
1117:   if (pdipm->Nh) {
1118:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
1119:     MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
1120:   }

1122:   /* Count dnz,onz for preallocation of KKT matrix */
1123:   jac_equality_trans   = pdipm->jac_equality_trans;
1124:   jac_inequality_trans = pdipm->jac_inequality_trans;
1125:   nce_all = pdipm->nce_all;

1127:   if (pdipm->Nxfixed) {
1128:     MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
1129:   }
1130:   MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);

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

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

1138:   /* Insert tao->hessian */
1139:   MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
1140:   for (i=0; i<pdipm->nx; i++){
1141:     row = rstart + i;

1143:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1144:     proc = 0;
1145:     for (j=0; j < nc; j++) {
1146:       while (aj[j] >= cranges[proc+1]) proc++;
1147:       col = aj[j] - cranges[proc] + Jranges[proc];
1148:       MatPreallocateSet(row,1,&col,dnz,onz);
1149:     }
1150:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);

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

1167:     /* Insert Jce_xfixed^T' */
1168:     if (pdipm->nxfixed) {
1169:       MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1170:       MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1171:       proc = 0;
1172:       for (j=0; j < nc; j++) {
1173:         /* find row ownership of */
1174:         while (aj[j] >= ranges[proc+1]) proc++;
1175:         nx_all = rranges[proc+1] - rranges[proc];
1176:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1177:         MatPreallocateSet(row,1,&col,dnz,onz);
1178:       }
1179:       MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1180:     }

1182:     if (pdipm->nh) {
1183:       /* Insert -grad h' */
1184:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1185:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1186:       proc = 0;
1187:       for (j=0; j < nc; j++) {
1188:         /* find row ownership of */
1189:         while (aj[j] >= ranges[proc+1]) proc++;
1190:         nx_all = rranges[proc+1] - rranges[proc];
1191:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1192:         MatPreallocateSet(row,1,&col,dnz,onz);
1193:       }
1194:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1195:     }

1197:     /* Insert Jci_xb^T' */
1198:     MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1199:     MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1200:     proc = 0;
1201:     for (j=0; j < nc; j++) {
1202:       /* find row ownership of */
1203:       while (aj[j] >= ranges[proc+1]) proc++;
1204:       nx_all = rranges[proc+1] - rranges[proc];
1205:       col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1206:       MatPreallocateSet(row,1,&col,dnz,onz);
1207:     }
1208:     MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1209:   }

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

1217:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1218:       proc = 0;
1219:       for (j=0; j < nc; j++) {
1220:         while (aj[j] >= cranges[proc+1]) proc++;
1221:         col = aj[j] - cranges[proc] + Jranges[proc];
1222:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1223:       }
1224:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1225:     }
1226:   }
1227:   /* Jce_xfixed */
1228:   if (pdipm->Nxfixed) {
1229:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1230:     for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1231:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

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

1236:       proc = 0;
1237:       j    = 0;
1238:       while (cols[j] >= cranges[proc+1]) proc++;
1239:       col = cols[j] - cranges[proc] + Jranges[proc];
1240:       MatPreallocateSet(row,1,&col,dnz,onz);
1241:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1242:     }
1243:   }

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

1251:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1252:       proc = 0;
1253:       for (j=0; j < nc; j++) {
1254:         while (aj[j] >= cranges[proc+1]) proc++;
1255:         col = aj[j] - cranges[proc] + Jranges[proc];
1256:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1257:       }
1258:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1259:     }
1260:     /* I */
1261:     for (i=0; i < pdipm->nh; i++){
1262:       row = rstart + pdipm->off_lambdai + i;
1263:       col = rstart + pdipm->off_z + i;
1264:       MatPreallocateSet(row,1,&col,dnz,onz);
1265:     }
1266:   }

1268:   /* Jci_xb */
1269:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1270:   for (i=0; i < (pdipm->nci - pdipm->nh); i++){
1271:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1273:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1274:     if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1275:     proc = 0;
1276:     for (j=0; j < nc; j++) {
1277:       while (cols[j] >= cranges[proc+1]) proc++;
1278:       col = cols[j] - cranges[proc] + Jranges[proc];
1279:       MatPreallocateSet(row,1,&col,dnz,onz);
1280:     }
1281:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1282:     /* I */
1283:     col = rstart + pdipm->off_z + pdipm->nh + i;
1284:     MatPreallocateSet(row,1,&col,dnz,onz);
1285:   }

1287:   /* 4-th Row block of KKT matrix: Z and Ci */
1288:   for (i=0; i < pdipm->nci; i++) {
1289:     row     = rstart + pdipm->off_z + i;
1290:     cols1[0] = rstart + pdipm->off_lambdai + i;
1291:     cols1[1] = row;
1292:     MatPreallocateSet(row,2,cols1,dnz,onz);
1293:   }

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

1298:   /* Create KKT matrix */
1299:   MatCreate(comm,&J);
1300:   MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1301:   MatSetFromOptions(J);
1302:   MatSeqAIJSetPreallocation(J,0,dnz);
1303:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1304:   /* MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); */
1305:   MatPreallocateFinalize(dnz,onz);
1306:   pdipm->K = J;

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

1312:   if (pdipm->solve_reduced_kkt) {
1313:     PC pc;
1314:     KSPGetPC(tao->ksp,&pc);
1315:     PCSetType(pc,PCFIELDSPLIT);
1316:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1317:     PCFieldSplitSetIS(pc,"2",pdipm->is2);
1318:     PCFieldSplitSetIS(pc,"1",pdipm->is1);
1319:   }
1320:   SNESSetFromOptions(pdipm->snes);

1322:   /* (9) Insert constant entries to  K */
1323:   /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1324:   MatGetOwnershipRange(J,&rstart,&rend);
1325:   for (i=rstart; i<rend; i++){
1326:     MatSetValue(J,i,i,0.0,INSERT_VALUES);
1327:   }
1328:   /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1329:   if (pdipm->kkt_pd){
1330:       for (i=0; i<pdipm->nh; i++){
1331:         row  = rstart + i;
1332:         MatSetValue(J,row,row,pdipm->deltaw,INSERT_VALUES);
1333:       }
1334:   }

1336:   /* Row block of K: [ grad Ce, 0, 0, 0] */
1337:   if (pdipm->Nxfixed) {
1338:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1339:     for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1340:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1342:       MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1343:       proc = 0;
1344:       for (j=0; j < nc; j++) {
1345:         while (cols[j] >= cranges[proc+1]) proc++;
1346:         col = cols[j] - cranges[proc] + Jranges[proc];
1347:         MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1348:         MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1349:       }
1350:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1351:     }
1352:   }

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

1359:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1360:     proc = 0;
1361:     for (j=0; j < nc; j++) {
1362:       while (cols[j] >= cranges[proc+1]) proc++;
1363:       col = cols[j] - cranges[proc] + Jranges[proc];
1364:       MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1365:       MatSetValue(J,row,col,-aa[j],INSERT_VALUES);
1366:     }
1367:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);

1369:     col = rstart + pdipm->off_z + pdipm->nh + i;
1370:     MatSetValue(J,row,col,1,INSERT_VALUES);
1371:   }

1373:   for (i=0; i < pdipm->nh; i++){
1374:     row = rstart + pdipm->off_lambdai + i;
1375:     col = rstart + pdipm->off_z + i;
1376:     MatSetValue(J,row,col,1,INSERT_VALUES);
1377:   }

1379:   /* Row block of K: [ 0, 0, I, ...] */
1380:   for (i=0; i < pdipm->nci; i++){
1381:     row = rstart + pdipm->off_z + i;
1382:     col = rstart + pdipm->off_lambdai + i;
1383:     MatSetValue(J,row,col,1,INSERT_VALUES);
1384:   }

1386:   if (pdipm->Nxfixed) {
1387:     MatDestroy(&Jce_xfixed_trans);
1388:   }
1389:   MatDestroy(&Jci_xb_trans);
1390:   PetscFree3(ng_all,nh_all,Jranges);
1391:   return(0);
1392: }

1394: /*
1395:    TaoDestroy_PDIPM - Destroys the pdipm object

1397:    Input:
1398:    full pdipm

1400:    Output:
1401:    Destroyed pdipm
1402: */
1403: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1404: {
1405:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1409:   /* Freeing Vectors assocaiated with KKT (X) */
1410:   VecDestroy(&pdipm->x); /* Solution x */
1411:   VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1412:   VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1413:   VecDestroy(&pdipm->z);       /* Slack variables */
1414:   VecDestroy(&pdipm->X);       /* Big KKT system vector [x; lambdae; lambdai; z] */

1416:   /* work vectors */
1417:   VecDestroy(&pdipm->lambdae_xfixed);
1418:   VecDestroy(&pdipm->lambdai_xb);

1420:   /* Legrangian equality and inequality Vec */
1421:   VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1422:   VecDestroy(&pdipm->ci); /* Vec of inequality constraints */

1424:   /* Matrices */
1425:   MatDestroy(&pdipm->Jce_xfixed);
1426:   MatDestroy(&pdipm->Jci_xb); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1427:   MatDestroy(&pdipm->K);

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

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

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

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

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

1450:   if (pdipm->solve_reduced_kkt) {
1451:     ISDestroy(&pdipm->is1);
1452:     ISDestroy(&pdipm->is2);
1453:   }

1455:   /* SNES */
1456:   SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1457:   PetscFree(pdipm->nce_all);
1458:   MatDestroy(&pdipm->jac_equality_trans);
1459:   MatDestroy(&pdipm->jac_inequality_trans);

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

1464:   /* Destroy Dual */
1465:   VecDestroy(&tao->DE); /* equality dual */
1466:   VecDestroy(&tao->DI); /* dinequality dual */
1467:   return(0);
1468: }

1470: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1471: {
1472:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1476:   PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1477:   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);
1478:   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);
1479:   PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1480:   PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1481:   PetscOptionsBool("-tao_pdipm_symmetric_kkt","Solve non reduced symmetric KKT system",NULL,pdipm->solve_symmetric_kkt,&pdipm->solve_symmetric_kkt,NULL);
1482:   PetscOptionsBool("-tao_pdipm_kkt_shift_pd","Add shifts to make KKT matrix positive definite",NULL,pdipm->kkt_pd,&pdipm->kkt_pd,NULL);
1483:   PetscOptionsTail();
1484:   return(0);
1485: }

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

1490:   Option Database Keys:
1491: +   -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1492: .   -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1493: .   -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1494: .   -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1495: -   -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite

1497:   Level: beginner
1498: M*/
1499: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1500: {
1501:   TAO_PDIPM      *pdipm;

1505:   tao->ops->setup          = TaoSetup_PDIPM;
1506:   tao->ops->solve          = TaoSolve_PDIPM;
1507:   tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1508:   tao->ops->view           = TaoView_PDIPM;
1509:   tao->ops->destroy        = TaoDestroy_PDIPM;

1511:   PetscNewLog(tao,&pdipm);
1512:   tao->data = (void*)pdipm;

1514:   pdipm->nx      = pdipm->Nx      = 0;
1515:   pdipm->nxfixed = pdipm->Nxfixed = 0;
1516:   pdipm->nxlb    = pdipm->Nxlb    = 0;
1517:   pdipm->nxub    = pdipm->Nxub    = 0;
1518:   pdipm->nxbox   = pdipm->Nxbox   = 0;
1519:   pdipm->nxfree  = pdipm->Nxfree  = 0;

1521:   pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1522:   pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1523:   pdipm->n  = pdipm->N  = 0;
1524:   pdipm->mu = 1.0;
1525:   pdipm->mu_update_factor = 0.1;

1527:   pdipm->deltaw     = 0.0;
1528:   pdipm->lastdeltaw = 3*1.e-4;
1529:   pdipm->deltac     = 0.0;
1530:   pdipm->kkt_pd     = PETSC_FALSE;

1532:   pdipm->push_init_slack     = 1.0;
1533:   pdipm->push_init_lambdai   = 1.0;
1534:   pdipm->solve_reduced_kkt   = PETSC_FALSE;
1535:   pdipm->solve_symmetric_kkt = PETSC_TRUE;

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

1541:   SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1542:   SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1543:   SNESGetKSP(pdipm->snes,&tao->ksp);
1544:   PetscObjectReference((PetscObject)tao->ksp);
1545:   return(0);
1546: }