Actual source code: pdipm.c

petsc-master 2020-08-10
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: {
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:   VecSet(pdipm->lambdae,0.0);

255:   /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
256:   VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
257:   VecSet(pdipm->z,pdipm->push_init_slack);

259:   /* Additional modification for X.lambdai and X.z */
260:   VecGetArray(pdipm->lambdai,&lambdai);
261:   VecGetArray(pdipm->z,&z);
262:   if(pdipm->Nh) {
263:     VecGetArrayRead(tao->constraints_inequality,&h);
264:     for (i=0; i < pdipm->nh; i++) {
265:       if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
266:       if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
267:     }
268:     VecRestoreArrayRead(tao->constraints_inequality,&h);
269:   }
270:   VecRestoreArray(pdipm->lambdai,&lambdai);
271:   VecRestoreArray(pdipm->z,&z);

273:   VecRestoreArray(pdipm->X,&Xarr);
274:   return(0);
275: }

277: /*
278:    TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X

280:    Input Parameter:
281:    snes - SNES context
282:    X - KKT Vector
283:    *ctx - pdipm context

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

304:   PetscObjectGetComm((PetscObject)snes,&comm);
305:   MPI_Comm_rank(comm,&rank);
306:   MPI_Comm_rank(comm,&size);

308:   MatGetOwnershipRanges(Jpre,&Jranges);
309:   MatGetOwnershipRange(Jpre,&Jrstart,NULL);
310:   MatGetOwnershipRangesColumn(tao->hessian,&rranges);
311:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

313:   VecGetArrayRead(X,&Xarr);

315:   /* (2) insert Z and Ci to Jpre -- overwrite existing values */
316:   for (i=0; i < pdipm->nci; i++) {
317:     row     = Jrstart + pdipm->off_z + i;
318:     cols[0] = Jrstart + pdipm->off_lambdai + i;
319:     cols[1] = row;
320:     vals[0] = Xarr[pdipm->off_z + i];
321:     vals[1] = Xarr[pdipm->off_lambdai + i];
322:     MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
323:   }

325:   /* (3) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
326:   if(pdipm->Ng) {
327:     MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
328:     for (i=0; i<pdipm->ng; i++){
329:       row = Jrstart + pdipm->off_lambdae + i;
330:       
331:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
332:       proc = 0;
333:       for (j=0; j < nc; j++) {
334:         while (aj[j] >= cranges[proc+1]) proc++;
335:         cols[0] = aj[j] - cranges[proc] + Jranges[proc];
336:         MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
337:       }
338:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
339:     }
340:   }

342:   if(pdipm->Nh) {
343:     /* (4) insert 3nd row block of Jpre: [ grad h, 0, 0, 0] */
344:     MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
345:     for (i=0; i < pdipm->nh; i++){
346:       row = Jrstart + pdipm->off_lambdai + i;
347:       
348:       MatGetRow(tao->jacobian_inequality,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_inequality,i+rjstart,&nc,&aj,&aa);
356:     }
357:   }

359:   /* (5) insert Wxx, grad g' and -grad h' to Jpre */
360:   if(pdipm->Ng) {
361:     MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
362:   }
363:   if(pdipm->Nh) {
364:     MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
365:   }

367:   VecPlaceArray(pdipm->x,Xarr);
368:   TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
369:   VecResetArray(pdipm->x);

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

375:     /* insert Wxx */
376:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
377:     proc = 0;
378:     for (j=0; j < nc; j++) {
379:       while (aj[j] >= cranges[proc+1]) proc++;
380:       cols[0] = aj[j] - cranges[proc] + Jranges[proc];
381:       MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
382:     }
383:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);

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

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

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

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

428: /*
429:    TaoSnesFunction_PDIPM - Evaluate KKT function at X

431:    Input Parameter:
432:    snes - SNES context
433:    X - KKT Vector
434:    *ctx - pdipm

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

451:   VecSet(F,0.0);

453:   VecGetArrayRead(X,&Xarr);
454:   VecGetArray(F,&Farr);

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

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

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

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

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

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

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

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

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

527:   VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
528:   VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
529:   VecResetArray(pdipm->ci);

531:   /* note: pdipm->z is not changed below */
532:   VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
533:   VecNorm(pdipm->z,NORM_2,&res[1]);
534:   VecResetArray(pdipm->z);

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

539:   VecRestoreArrayRead(X,&Xarr);
540:   VecRestoreArray(F,&Farr);
541:   return(0);
542: }

544: /*
545:    PDIPMLineSearch - Custom line search used with PDIPM.

547:    Collective on TAO

549:    Notes:
550:    PDIPMLineSearch employs a simple backtracking line-search to keep
551:    the slack variables (z) and inequality constraints lagrange multipliers
552:    (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
553:    alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the 
554:    slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
555:    are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
556:    is also updated as mu = mu + z'lambdai/Nci
557: */
558: PetscErrorCode PDIPMLineSearch(SNESLineSearch linesearch,void *ctx)
559: {
561:   Tao            tao=(Tao)ctx;
562:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;
563:   SNES           snes;
564:   Vec            X,F,Y,W,G;
565:   PetscInt       i,iter;
566:   PetscReal      alpha_p=1.0,alpha_d=1.0,alpha[4];
567:   PetscScalar    *Xarr,*z,*lambdai,dot;
568:   const PetscScalar *dXarr,*dz,*dlambdai;
569:   PetscScalar    *taosolarr;

572:   SNESLineSearchGetSNES(linesearch,&snes);
573:   SNESGetIterationNumber(snes,&iter);

575:   SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
576:   SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);

578:   VecGetArray(X,&Xarr);
579:   VecGetArrayRead(Y,&dXarr);
580:   z  = Xarr + pdipm->off_z;
581:   dz = dXarr + pdipm->off_z;
582:   for (i=0; i < pdipm->nci; i++) {
583:     if (z[i] - dz[i] < 0.0) {
584:       alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
585:     }
586:   }

588:   lambdai  = Xarr + pdipm->off_lambdai;
589:   dlambdai = dXarr + pdipm->off_lambdai;

591:   for (i=0; i<pdipm->nci; i++) {
592:     if (lambdai[i] - dlambdai[i] < 0.0) {
593:       alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
594:     }
595:   }

597:   alpha[0] = alpha_p;
598:   alpha[1] = alpha_d;
599:   VecRestoreArrayRead(Y,&dXarr);
600:   VecRestoreArray(X,&Xarr);

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

605:   alpha_p = alpha[2];
606:   alpha_d = alpha[3];

608:   VecGetArray(X,&Xarr);
609:   VecGetArrayRead(Y,&dXarr);
610:   for (i=0; i<pdipm->nx; i++) {
611:     Xarr[i] = Xarr[i] - alpha_p * dXarr[i];
612:   }

614:   for (i=0; i<pdipm->nce; i++) {
615:     Xarr[i+pdipm->off_lambdae] = Xarr[i+pdipm->off_lambdae] - alpha_d * dXarr[i+pdipm->off_lambdae];
616:   }

618:   for (i=0; i<pdipm->nci; i++) {
619:     Xarr[i+pdipm->off_lambdai] = Xarr[i+pdipm->off_lambdai] - alpha_d * dXarr[i+pdipm->off_lambdai];
620:   }

622:   for (i=0; i<pdipm->nci; i++) {
623:     Xarr[i+pdipm->off_z] = Xarr[i+pdipm->off_z] - alpha_p * dXarr[i+pdipm->off_z];
624:   }

626:   VecGetArray(tao->solution,&taosolarr);
627:   PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
628:   VecRestoreArray(tao->solution,&taosolarr);


631:   VecRestoreArray(X,&Xarr);
632:   VecRestoreArrayRead(Y,&dXarr);

634:   /* Evaluate F at X */
635:   SNESComputeFunction(snes,X,F);
636:   SNESLineSearchComputeNorms(linesearch); /* must call this func, do not know why */

638:   /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
639:   VecDot(pdipm->z,pdipm->lambdai,&dot);

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

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

647:   tao->niter++;
648:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
649:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);

651:   (*tao->ops->convergencetest)(tao,tao->cnvP);
652:   if (tao->reason) {
653:     SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
654:   }
655:   return(0);
656: }

658: /*
659:    TaoSolve_PDIPM

661:    Input Parameter:
662:    tao - TAO context

664:    Output Parameter:
665:    tao - TAO context
666: */
667: PetscErrorCode TaoSolve_PDIPM(Tao tao)
668: {
669:   PetscErrorCode     ierr;
670:   TAO_PDIPM          *pdipm = (TAO_PDIPM*)tao->data;
671:   SNESLineSearch     linesearch;  /* SNESLineSearch context */
672:   Vec                dummy;

675:   /* Initialize all variables */
676:   TaoPDIPMInitializeSolution(tao);

678:   /* Set linesearch */
679:   SNESGetLineSearch(pdipm->snes,&linesearch);
680:   SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
681:   SNESLineSearchShellSetUserFunc(linesearch,PDIPMLineSearch,tao);
682:   SNESLineSearchSetFromOptions(linesearch);

684:   tao->reason = TAO_CONTINUE_ITERATING;

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

690:   TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
691:   TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
692:   VecDestroy(&dummy);
693:   (*tao->ops->convergencetest)(tao,tao->cnvP);
694:   if (tao->reason) {
695:     SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
696:   }

698:   while (tao->reason == TAO_CONTINUE_ITERATING) {
699:     SNESConvergedReason reason;
700:     SNESSolve(pdipm->snes,NULL,pdipm->X);

702:     /* Check SNES convergence */
703:     SNESGetConvergedReason(pdipm->snes,&reason);
704:     if (reason < 0) {
705:       PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
706:     }

708:     /* Check TAO convergence */
709:     if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
710:   }
711:   return(0);
712: }

714: /*
715:    TaoSetup_PDIPM - Sets up tao and pdipm

717:    Input Parameter:
718:    tao - TAO object

720:    Output:   pdipm - initialized object
721: */
722: PetscErrorCode TaoSetup_PDIPM(Tao tao)
723: {
724:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;
726:   MPI_Comm       comm;
727:   PetscMPIInt    rank,size;
728:   PetscInt       row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
729:   PetscInt       offset,*xa,*xb,i,j,rstart,rend;
730:   PetscScalar    one=1.0,neg_one=-1.0,*Xarr;
731:   const PetscInt    *cols,*rranges,*cranges,*aj,*ranges;
732:   const PetscScalar *aa;
733:   Mat            J,jac_equality_trans,jac_inequality_trans;
734:   Mat            Jce_xfixed_trans,Jci_xb_trans;
735:   PetscInt       *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];

738:   PetscObjectGetComm((PetscObject)tao,&comm);
739:   MPI_Comm_rank(comm,&rank);
740:   MPI_Comm_size(comm,&size);

742:   /* (1) Setup Bounds and create Tao vectors */
743:   TaoPDIPMSetUpBounds(tao);

745:   if (!tao->gradient) {
746:     VecDuplicate(tao->solution,&tao->gradient);
747:     VecDuplicate(tao->solution,&tao->stepdirection);
748:   }

750:   /* (2) Get sizes */
751:   /* Size of vector x - This is set by TaoSetInitia√lVector */
752:   VecGetSize(tao->solution,&pdipm->Nx);
753:   VecGetLocalSize(tao->solution,&pdipm->nx);

755:   /* Size of equality constraints and vectors */
756:   if (tao->constraints_equality) {
757:     VecGetSize(tao->constraints_equality,&pdipm->Ng);
758:     VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
759:   } else {
760:     pdipm->ng = pdipm->Ng = 0;
761:   }

763:   pdipm->nce = pdipm->ng + pdipm->nxfixed;
764:   pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;

766:   /* Size of inequality constraints and vectors */
767:   if (tao->constraints_inequality) {
768:     VecGetSize(tao->constraints_inequality,&pdipm->Nh);
769:     VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
770:   } else {
771:     pdipm->nh = pdipm->Nh = 0;
772:   }

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

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

781:   /* list below to TaoView_PDIPM()? */
782:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] nce %d = ng %d + nxfixed %d\n",rank,pdipm->nce,pdipm->ng,pdipm->nxfixed); */
783:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] nci %d = nh %d + nxlb %d + nxub %d + 2*nxbox %d\n",rank,pdipm->nci,pdipm->nh,pdipm->nxlb,pdipm->nxub,pdipm->nxbox); */
784:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] n %d = nx %d + nce %d + 2*nci %d\n",rank,pdipm->n,pdipm->nx,pdipm->nce,pdipm->nci); */

786:   /* (3) Offsets for subvectors */
787:   pdipm->off_lambdae = pdipm->nx;
788:   pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
789:   pdipm->off_z       = pdipm->off_lambdai + pdipm->nci;

791:   /* (4) Create vectors and subvectors */
792:   /* Ce and Ci vectors */
793:   VecCreate(comm,&pdipm->ce);
794:   VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
795:   VecSetFromOptions(pdipm->ce);

797:   VecCreate(comm,&pdipm->ci);
798:   VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
799:   VecSetFromOptions(pdipm->ci);

801:   /* X=[x; lambdae; lambdai; z] for the big KKT system */
802:   VecCreate(comm,&pdipm->X);
803:   VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
804:   VecSetFromOptions(pdipm->X);

806:   /* Subvectors; they share local arrays with X */
807:   VecGetArray(pdipm->X,&Xarr);
808:   /* x shares local array with X.x */
809:   if (pdipm->Nx) {
810:     VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
811:   }

813:   /* lambdae shares local array with X.lambdae */
814:   if (pdipm->Nce) {
815:     VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
816:   }

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

822:     VecCreate(comm,&pdipm->lambdae_xfixed);
823:     VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
824:     VecSetFromOptions(pdipm->lambdae_xfixed);
825:   }

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

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

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

840:   VecCreate(comm,&pdipm->lambdai_xb);
841:   VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
842:   VecSetFromOptions(pdipm->lambdai_xb);

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

846:   /* (5) Create Jacobians Jce_xfixed and Jci */
847:   /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
848:   if (pdipm->Nxfixed) {
849:     /* Create Jce_xfixed */
850:     MatCreate(comm,&pdipm->Jce_xfixed);
851:     MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
852:     MatSetFromOptions(pdipm->Jce_xfixed);
853:     MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
854:     MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);

856:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
857:     ISGetIndices(pdipm->isxfixed,&cols);
858:     k = 0;
859:     for (row = Jcrstart; row < Jcrend; row++) {
860:       MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
861:       k++;
862:     }
863:     ISRestoreIndices(pdipm->isxfixed, &cols);
864:     MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
865:     MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
866:   }

868:   /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
869:   MatCreate(comm,&pdipm->Jci_xb);
870:   MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
871:   MatSetFromOptions(pdipm->Jci_xb);
872:   MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
873:   MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);

875:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
876:   offset = Jcrstart;
877:   if (pdipm->Nxub) {
878:     /* Add xub to Jci_xb */
879:     ISGetIndices(pdipm->isxub,&cols);
880:     k = 0;
881:     for (row = offset; row < offset + pdipm->nxub; row++) {
882:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
883:       k++;
884:     }
885:     ISRestoreIndices(pdipm->isxub, &cols);
886:   }

888:   if (pdipm->Nxlb) {
889:     /* Add xlb to Jci_xb */
890:     ISGetIndices(pdipm->isxlb,&cols);
891:     k = 0;
892:     offset += pdipm->nxub;
893:     for (row = offset; row < offset + pdipm->nxlb; row++) {
894:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
895:       k++;
896:     }
897:     ISRestoreIndices(pdipm->isxlb, &cols);
898:   }

900:   /* Add xbox to Jci_xb */
901:   if (pdipm->Nxbox) {
902:     ISGetIndices(pdipm->isxbox,&cols);
903:     k = 0;
904:     offset += pdipm->nxlb;
905:     for (row = offset; row < offset + pdipm->nxbox; row++) {
906:       MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
907:       tmp = row + pdipm->nxbox;
908:       MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
909:       k++;
910:     }
911:     ISRestoreIndices(pdipm->isxbox, &cols);
912:   }

914:   MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
915:   MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
916:   /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */

918:   /* (6) Set up ISs for PC Fieldsplit */
919:   if (pdipm->solve_reduced_kkt) {
920:     PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
921:     for(i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
922:     for(i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;

924:     ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
925:     ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
926:   }

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

931:   /* Get rstart of KKT matrix */
932:   MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
933:   rstart -= pdipm->n;

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

937:   PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
938:   MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
939:   MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
940:   MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);

942:   MatGetOwnershipRanges(tao->hessian,&rranges);
943:   MatGetOwnershipRangesColumn(tao->hessian,&cranges);

945:   if (pdipm->Ng) {
946:     TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
947:     MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
948:   }
949:   if (pdipm->Nh) {
950:     TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
951:     MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
952:   }

954:   /* Count dnz,onz for preallocation of KKT matrix */
955:   jac_equality_trans   = pdipm->jac_equality_trans;
956:   jac_inequality_trans = pdipm->jac_inequality_trans;
957:   nce_all = pdipm->nce_all;

959:   if (pdipm->Nxfixed) {
960:     MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
961:   }
962:   MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);

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

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

970:   /* Insert tao->hessian */
971:   MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
972:   for (i=0; i<pdipm->nx; i++){
973:     row = rstart + i;

975:     MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
976:     proc = 0;
977:     for (j=0; j < nc; j++) {
978:       while (aj[j] >= cranges[proc+1]) proc++;
979:       col = aj[j] - cranges[proc] + Jranges[proc];
980:       MatPreallocateSet(row,1,&col,dnz,onz);
981:     }
982:     MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);

984:     if(pdipm->ng) {
985:       /* Insert grad g' */
986:       MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
987:       MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
988:       proc = 0;
989:       for (j=0; j < nc; j++) {
990:         /* find row ownership of */
991:         while (aj[j] >= ranges[proc+1]) proc++;
992:         nx_all = rranges[proc+1] - rranges[proc];
993:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
994:         MatPreallocateSet(row,1,&col,dnz,onz);
995:       }
996:       MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
997:     }

999:     /* Insert Jce_xfixed^T' */
1000:     if (pdipm->nxfixed) {
1001:       MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1002:       MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1003:       proc = 0;
1004:       for (j=0; j < nc; j++) {
1005:         /* find row ownership of */
1006:         while (aj[j] >= ranges[proc+1]) proc++;
1007:         nx_all = rranges[proc+1] - rranges[proc];
1008:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1009:         MatPreallocateSet(row,1,&col,dnz,onz);
1010:       }
1011:       MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1012:     }

1014:     if(pdipm->nh) {
1015:       /* Insert -grad h' */
1016:       MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1017:       MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1018:       proc = 0;
1019:       for (j=0; j < nc; j++) {
1020:         /* find row ownership of */
1021:         while (aj[j] >= ranges[proc+1]) proc++;
1022:         nx_all = rranges[proc+1] - rranges[proc];
1023:         col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1024:         MatPreallocateSet(row,1,&col,dnz,onz);
1025:       }
1026:       MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1027:     }

1029:     /* Insert Jci_xb^T' */
1030:     MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1031:     MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1032:     proc = 0;
1033:     for (j=0; j < nc; j++) {
1034:       /* find row ownership of */
1035:       while (aj[j] >= ranges[proc+1]) proc++;
1036:       nx_all = rranges[proc+1] - rranges[proc];
1037:       col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1038:       MatPreallocateSet(row,1,&col,dnz,onz);
1039:     }
1040:     MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1041:   }

1043:   /* 2nd Row block of KKT matrix: [grad Ce, 0, 0, 0] */
1044:   if(pdipm->Ng) {
1045:     MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
1046:     for (i=0; i < pdipm->ng; i++){
1047:       row = rstart + pdipm->off_lambdae + i;

1049:       MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1050:       proc = 0;
1051:       for (j=0; j < nc; j++) {
1052:         while (aj[j] >= cranges[proc+1]) proc++;
1053:         col = aj[j] - cranges[proc] + Jranges[proc];
1054:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1055:       }
1056:       MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1057:     }
1058:   }
1059:   /* Jce_xfixed */
1060:   if (pdipm->Nxfixed) {
1061:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1062:     for (i=0; i < (pdipm->nce - pdipm->ng); i++ ){
1063:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

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

1068:       proc = 0;
1069:       j    = 0;
1070:       while (cols[j] >= cranges[proc+1]) proc++;
1071:       col = cols[j] - cranges[proc] + Jranges[proc];
1072:       MatPreallocateSet(row,1,&col,dnz,onz);
1073:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1074:     }
1075:   }

1077:   /* 3rd Row block of KKT matrix: [ gradCi, 0, 0, -I] */
1078:   if(pdipm->Nh) {
1079:     MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
1080:     for (i=0; i < pdipm->nh; i++){
1081:       row = rstart + pdipm->off_lambdai + i;

1083:       MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1084:       proc = 0;
1085:       for (j=0; j < nc; j++) {
1086:         while (aj[j] >= cranges[proc+1]) proc++;
1087:         col = aj[j] - cranges[proc] + Jranges[proc];
1088:         MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1089:       }
1090:       MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1091:     }
1092:     /* -I */
1093:     for (i=0; i < pdipm->nh; i++){
1094:       row = rstart + pdipm->off_lambdai + i;
1095:       col = rstart + pdipm->off_z + i;
1096:       MatPreallocateSet(row,1,&col,dnz,onz);
1097:     }
1098:   }

1100:   /* Jci_xb */
1101:   MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1102:   for (i=0; i < (pdipm->nci - pdipm->nh); i++ ){
1103:     row = rstart + pdipm->off_lambdai + pdipm->nh + i;

1105:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1106:     if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1107:     proc = 0;
1108:     for (j=0; j < nc; j++) {
1109:       while (cols[j] >= cranges[proc+1]) proc++;
1110:       col = cols[j] - cranges[proc] + Jranges[proc];
1111:       MatPreallocateSet(row,1,&col,dnz,onz);
1112:     }
1113:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1114:     /* -I */
1115:     col = rstart + pdipm->off_z + pdipm->nh + i;
1116:     MatPreallocateSet(row,1,&col,dnz,onz);
1117:   }

1119:   /* 4-th Row block of KKT matrix: Z and Ci */
1120:   for (i=0; i < pdipm->nci; i++) {
1121:     row     = rstart + pdipm->off_z + i;
1122:     cols1[0] = rstart + pdipm->off_lambdai + i;
1123:     cols1[1] = row;
1124:     MatPreallocateSet(row,2,cols1,dnz,onz);
1125:   }

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

1130:   /* Create KKT matrix */
1131:   MatCreate(comm,&J);
1132:   MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1133:   MatSetFromOptions(J);
1134:   MatSeqAIJSetPreallocation(J,0,dnz);
1135:   MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1136:   /* MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); */
1137:   MatPreallocateFinalize(dnz,onz);
1138:   pdipm->K = J;

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

1144:   if (pdipm->solve_reduced_kkt) {
1145:     PC pc;
1146:     KSPGetPC(tao->ksp,&pc);
1147:     PCSetType(pc,PCFIELDSPLIT);
1148:     PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1149:     PCFieldSplitSetIS(pc,"2",pdipm->is2);
1150:     PCFieldSplitSetIS(pc,"1",pdipm->is1);
1151:   }
1152:   SNESSetFromOptions(pdipm->snes);

1154:   /* (9) Insert constant entries to  K */
1155:   /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1156:   MatGetOwnershipRange(J,&rstart,&rend);
1157:   for (i=rstart; i<rend; i++){
1158:     MatSetValue(J,i,i,0.0,INSERT_VALUES);
1159:   }

1161:   /* Row block of K: [ grad Ce, 0, 0, 0] */
1162:   if (pdipm->Nxfixed) {
1163:     MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1164:     for (i=0; i < (pdipm->nce - pdipm->ng); i++ ){
1165:       row = rstart + pdipm->off_lambdae + pdipm->ng + i;

1167:       MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1168:       proc = 0;
1169:       for (j=0; j < nc; j++) {
1170:         while (cols[j] >= cranges[proc+1]) proc++;
1171:         col = cols[j] - cranges[proc] + Jranges[proc];
1172:         MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1173:         MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1174:       }
1175:       MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1176:     }
1177:   }

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

1184:     MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1185:     proc = 0;
1186:     for (j=0; j < nc; j++) {
1187:       while (cols[j] >= cranges[proc+1]) proc++;
1188:       col = cols[j] - cranges[proc] + Jranges[proc];
1189:       MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1190:       MatSetValue(J,row,col,aa[j],INSERT_VALUES);
1191:     }
1192:     MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);

1194:     col = rstart + pdipm->off_z + pdipm->nh + i;
1195:     MatSetValue(J,row,col,-1,INSERT_VALUES);
1196:   }

1198:   for (i=0; i < pdipm->nh; i++){
1199:     row = rstart + pdipm->off_lambdai + i;
1200:     col = rstart + pdipm->off_z + i;
1201:     MatSetValue(J,row,col,-1,INSERT_VALUES);
1202:   }

1204:   if (pdipm->Nxfixed) {
1205:     MatDestroy(&Jce_xfixed_trans);
1206:   }
1207:   MatDestroy(&Jci_xb_trans);
1208:   PetscFree3(ng_all,nh_all,Jranges);
1209:   return(0);
1210: }

1212: /*
1213:    TaoDestroy_PDIPM - Destroys the pdipm object

1215:    Input:
1216:    full pdipm

1218:    Output:
1219:    Destroyed pdipm
1220: */
1221: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1222: {
1223:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1227:   /* Freeing Vectors assocaiated with KKT (X) */
1228:   VecDestroy(&pdipm->x); /* Solution x */
1229:   VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1230:   VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1231:   VecDestroy(&pdipm->z);       /* Slack variables */
1232:   VecDestroy(&pdipm->X);       /* Big KKT system vector [x; lambdae; lambdai; z] */

1234:   /* work vectors */
1235:   VecDestroy(&pdipm->lambdae_xfixed);
1236:   VecDestroy(&pdipm->lambdai_xb);

1238:   /* Legrangian equality and inequality Vec */
1239:   VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1240:   VecDestroy(&pdipm->ci); /* Vec of inequality constraints */

1242:   /* Matrices */
1243:   MatDestroy(&pdipm->Jce_xfixed);
1244:   MatDestroy(&pdipm->Jci_xb); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1245:   MatDestroy(&pdipm->K);

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

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

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

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

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

1268:   if (pdipm->solve_reduced_kkt) {
1269:     ISDestroy(&pdipm->is1);
1270:     ISDestroy(&pdipm->is2);
1271:   }

1273:   /* SNES */
1274:   SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1275:   PetscFree(pdipm->nce_all);
1276:   MatDestroy(&pdipm->jac_equality_trans);
1277:   MatDestroy(&pdipm->jac_inequality_trans);

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

1282:   /* Destroy Dual */
1283:   VecDestroy(&tao->DE); /* equality dual */
1284:   VecDestroy(&tao->DI); /* dinequality dual */
1285:   return(0);
1286: }

1288: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1289: {
1290:   TAO_PDIPM      *pdipm = (TAO_PDIPM*)tao->data;

1294:   PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1295:   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);
1296:   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);
1297:   PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1298:   PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1299:   PetscOptionsTail();
1300:   return(0);
1301: }

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

1306:   Option Database Keys:
1307: +   -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1308: .   -tao_pdipm_push_init_slack  - parameter to push initial slack variables away from bounds (> 0)
1309: -   -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)

1311:   Level: beginner
1312: M*/
1313: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1314: {
1315:   TAO_PDIPM      *pdipm;

1319:   tao->ops->setup          = TaoSetup_PDIPM;
1320:   tao->ops->solve          = TaoSolve_PDIPM;
1321:   tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1322:   tao->ops->destroy        = TaoDestroy_PDIPM;

1324:   PetscNewLog(tao,&pdipm);
1325:   tao->data = (void*)pdipm;

1327:   pdipm->nx      = pdipm->Nx      = 0;
1328:   pdipm->nxfixed = pdipm->Nxfixed = 0;
1329:   pdipm->nxlb    = pdipm->Nxlb    = 0;
1330:   pdipm->nxub    = pdipm->Nxub    = 0;
1331:   pdipm->nxbox   = pdipm->Nxbox   = 0;
1332:   pdipm->nxfree  = pdipm->Nxfree  = 0;

1334:   pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1335:   pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1336:   pdipm->n  = pdipm->N  = 0;
1337:   pdipm->mu = 1.0;
1338:   pdipm->mu_update_factor = 0.1;

1340:   pdipm->push_init_slack   = 1.0;
1341:   pdipm->push_init_lambdai = 1.0;
1342:   pdipm->solve_reduced_kkt = PETSC_FALSE;

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

1348:   SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1349:   SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1350:   SNESGetKSP(pdipm->snes,&tao->ksp);
1351:   PetscObjectReference((PetscObject)tao->ksp);
1352:   return(0);
1353: }