Actual source code: pdipm.c
1: #include <../src/tao/constrained/impls/ipm/pdipm.h>
3: /*
4: TaoPDIPMEvaluateFunctionsAndJacobians - Evaluate the objective function f, gradient fx, constraints, and all the Jacobians at current vector
6: Collective on tao
8: Input Parameter:
9: + tao - solver context
10: - x - vector at which all objects to be evaluated
12: Level: beginner
14: .seealso: TaoPDIPMUpdateConstraints(), TaoPDIPMSetUpBounds()
15: */
16: static PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)
17: {
18: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
20: /* Compute user objective function and gradient */
21: TaoComputeObjectiveAndGradient(tao,x,&pdipm->obj,tao->gradient);
23: /* Equality constraints and Jacobian */
24: if (pdipm->Ng) {
25: TaoComputeEqualityConstraints(tao,x,tao->constraints_equality);
26: TaoComputeJacobianEquality(tao,x,tao->jacobian_equality,tao->jacobian_equality_pre);
27: }
29: /* Inequality constraints and Jacobian */
30: if (pdipm->Nh) {
31: TaoComputeInequalityConstraints(tao,x,tao->constraints_inequality);
32: TaoComputeJacobianInequality(tao,x,tao->jacobian_inequality,tao->jacobian_inequality_pre);
33: }
34: return 0;
35: }
37: /*
38: TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x
40: Collective
42: Input Parameter:
43: + tao - Tao context
44: - x - vector at which constraints to be evaluated
46: Level: beginner
48: .seealso: TaoPDIPMEvaluateFunctionsAndJacobians()
49: */
50: static PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao,Vec x)
51: {
52: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
53: PetscInt i,offset,offset1,k,xstart;
54: PetscScalar *carr;
55: const PetscInt *ubptr,*lbptr,*bxptr,*fxptr;
56: const PetscScalar *xarr,*xuarr,*xlarr,*garr,*harr;
58: VecGetOwnershipRange(x,&xstart,NULL);
60: VecGetArrayRead(x,&xarr);
61: VecGetArrayRead(tao->XU,&xuarr);
62: VecGetArrayRead(tao->XL,&xlarr);
64: /* (1) Update ce vector */
65: VecGetArrayWrite(pdipm->ce,&carr);
67: if (pdipm->Ng) {
68: /* (1.a) Inserting updated g(x) */
69: VecGetArrayRead(tao->constraints_equality,&garr);
70: PetscMemcpy(carr,garr,pdipm->ng*sizeof(PetscScalar));
71: VecRestoreArrayRead(tao->constraints_equality,&garr);
72: }
74: /* (1.b) Update xfixed */
75: if (pdipm->Nxfixed) {
76: offset = pdipm->ng;
77: ISGetIndices(pdipm->isxfixed,&fxptr); /* global indices in x */
78: for (k=0;k < pdipm->nxfixed;k++) {
79: i = fxptr[k]-xstart;
80: carr[offset + k] = xarr[i] - xuarr[i];
81: }
82: }
83: VecRestoreArrayWrite(pdipm->ce,&carr);
85: /* (2) Update ci vector */
86: VecGetArrayWrite(pdipm->ci,&carr);
88: if (pdipm->Nh) {
89: /* (2.a) Inserting updated h(x) */
90: VecGetArrayRead(tao->constraints_inequality,&harr);
91: PetscMemcpy(carr,harr,pdipm->nh*sizeof(PetscScalar));
92: VecRestoreArrayRead(tao->constraints_inequality,&harr);
93: }
95: /* (2.b) Update xub */
96: offset = pdipm->nh;
97: if (pdipm->Nxub) {
98: ISGetIndices(pdipm->isxub,&ubptr);
99: for (k=0; k<pdipm->nxub; k++) {
100: i = ubptr[k]-xstart;
101: carr[offset + k] = xuarr[i] - xarr[i];
102: }
103: }
105: if (pdipm->Nxlb) {
106: /* (2.c) Update xlb */
107: offset += pdipm->nxub;
108: ISGetIndices(pdipm->isxlb,&lbptr); /* global indices in x */
109: for (k=0; k<pdipm->nxlb; k++) {
110: i = lbptr[k]-xstart;
111: carr[offset + k] = xarr[i] - xlarr[i];
112: }
113: }
115: if (pdipm->Nxbox) {
116: /* (2.d) Update xbox */
117: offset += pdipm->nxlb;
118: offset1 = offset + pdipm->nxbox;
119: ISGetIndices(pdipm->isxbox,&bxptr); /* global indices in x */
120: for (k=0; k<pdipm->nxbox; k++) {
121: i = bxptr[k]-xstart; /* local indices in x */
122: carr[offset+k] = xuarr[i] - xarr[i];
123: carr[offset1+k] = xarr[i] - xlarr[i];
124: }
125: }
126: VecRestoreArrayWrite(pdipm->ci,&carr);
128: /* Restoring Vectors */
129: VecRestoreArrayRead(x,&xarr);
130: VecRestoreArrayRead(tao->XU,&xuarr);
131: VecRestoreArrayRead(tao->XL,&xlarr);
132: return 0;
133: }
135: /*
136: TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x
138: Collective
140: Input Parameter:
141: . tao - holds pdipm and XL & XU
143: Level: beginner
145: .seealso: TaoPDIPMUpdateConstraints
146: */
147: static PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
148: {
149: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
150: const PetscScalar *xl,*xu;
151: PetscInt n,*ixlb,*ixub,*ixfixed,*ixfree,*ixbox,i,low,high,idx;
152: MPI_Comm comm;
153: PetscInt sendbuf[5],recvbuf[5];
155: /* Creates upper and lower bounds vectors on x, if not created already */
156: TaoComputeVariableBounds(tao);
158: VecGetLocalSize(tao->XL,&n);
159: PetscMalloc5(n,&ixlb,n,&ixub,n,&ixfree,n,&ixfixed,n,&ixbox);
161: VecGetOwnershipRange(tao->XL,&low,&high);
162: VecGetArrayRead(tao->XL,&xl);
163: VecGetArrayRead(tao->XU,&xu);
164: for (i=0; i<n; i++) {
165: idx = low + i;
166: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
167: if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
168: ixfixed[pdipm->nxfixed++] = idx;
169: } else ixbox[pdipm->nxbox++] = idx;
170: } else {
171: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
172: ixlb[pdipm->nxlb++] = idx;
173: } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
174: ixub[pdipm->nxlb++] = idx;
175: } else ixfree[pdipm->nxfree++] = idx;
176: }
177: }
178: VecRestoreArrayRead(tao->XL,&xl);
179: VecRestoreArrayRead(tao->XU,&xu);
181: PetscObjectGetComm((PetscObject)tao,&comm);
182: sendbuf[0] = pdipm->nxlb;
183: sendbuf[1] = pdipm->nxub;
184: sendbuf[2] = pdipm->nxfixed;
185: sendbuf[3] = pdipm->nxbox;
186: sendbuf[4] = pdipm->nxfree;
188: MPI_Allreduce(sendbuf,recvbuf,5,MPIU_INT,MPI_SUM,comm);
189: pdipm->Nxlb = recvbuf[0];
190: pdipm->Nxub = recvbuf[1];
191: pdipm->Nxfixed = recvbuf[2];
192: pdipm->Nxbox = recvbuf[3];
193: pdipm->Nxfree = recvbuf[4];
195: if (pdipm->Nxlb) {
196: ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);
197: }
198: if (pdipm->Nxub) {
199: ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);
200: }
201: if (pdipm->Nxfixed) {
202: ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);
203: }
204: if (pdipm->Nxbox) {
205: ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);
206: }
207: if (pdipm->Nxfree) {
208: ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);
209: }
210: PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);
211: return 0;
212: }
214: /*
215: TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
216: X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
217: four subvectors need to be initialized and its values copied over to X. Instead
218: of copying, we use VecPlace/ResetArray functions to share the memory locations for
219: X and the subvectors
221: Collective
223: Input Parameter:
224: . tao - Tao context
226: Level: beginner
227: */
228: static PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
229: {
230: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
231: PetscScalar *Xarr,*z,*lambdai;
232: PetscInt i;
233: const PetscScalar *xarr,*h;
235: VecGetArrayWrite(pdipm->X,&Xarr);
237: /* Set Initialize X.x = tao->solution */
238: VecGetArrayRead(tao->solution,&xarr);
239: PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));
240: VecRestoreArrayRead(tao->solution,&xarr);
242: /* Initialize X.lambdae = 0.0 */
243: if (pdipm->lambdae) {
244: VecSet(pdipm->lambdae,0.0);
245: }
247: /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
248: if (pdipm->Nci) {
249: VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
250: VecSet(pdipm->z,pdipm->push_init_slack);
252: /* Additional modification for X.lambdai and X.z */
253: VecGetArrayWrite(pdipm->lambdai,&lambdai);
254: VecGetArrayWrite(pdipm->z,&z);
255: if (pdipm->Nh) {
256: VecGetArrayRead(tao->constraints_inequality,&h);
257: for (i=0; i < pdipm->nh; i++) {
258: if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
259: if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
260: }
261: VecRestoreArrayRead(tao->constraints_inequality,&h);
262: }
263: VecRestoreArrayWrite(pdipm->lambdai,&lambdai);
264: VecRestoreArrayWrite(pdipm->z,&z);
265: }
267: VecRestoreArrayWrite(pdipm->X,&Xarr);
268: return 0;
269: }
271: /*
272: TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X
274: Input Parameter:
275: snes - SNES context
276: X - KKT Vector
277: *ctx - pdipm context
279: Output Parameter:
280: J - Hessian matrix
281: Jpre - Preconditioner
282: */
283: static PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
284: {
285: Tao tao=(Tao)ctx;
286: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
287: PetscInt i,row,cols[2],Jrstart,rjstart,nc,j;
288: const PetscInt *aj,*ranges,*Jranges,*rranges,*cranges;
289: const PetscScalar *Xarr,*aa;
290: PetscScalar vals[2];
291: PetscInt proc,nx_all,*nce_all=pdipm->nce_all;
292: MPI_Comm comm;
293: PetscMPIInt rank,size;
294: Mat jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;
296: PetscObjectGetComm((PetscObject)snes,&comm);
297: MPI_Comm_rank(comm,&rank);
298: MPI_Comm_rank(comm,&size);
300: MatGetOwnershipRanges(Jpre,&Jranges);
301: MatGetOwnershipRange(Jpre,&Jrstart,NULL);
302: MatGetOwnershipRangesColumn(tao->hessian,&rranges);
303: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
305: VecGetArrayRead(X,&Xarr);
307: /* (1) insert Z and Ci to the 4th block of Jpre -- overwrite existing values */
308: if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
309: vals[0] = 1.0;
310: for (i=0; i < pdipm->nci; i++) {
311: row = Jrstart + pdipm->off_z + i;
312: cols[0] = Jrstart + pdipm->off_lambdai + i;
313: cols[1] = row;
314: vals[1] = Xarr[pdipm->off_lambdai + i]/Xarr[pdipm->off_z + i];
315: MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
316: }
317: } else {
318: for (i=0; i < pdipm->nci; i++) {
319: row = Jrstart + pdipm->off_z + i;
320: cols[0] = Jrstart + pdipm->off_lambdai + i;
321: cols[1] = row;
322: vals[0] = Xarr[pdipm->off_z + i];
323: vals[1] = Xarr[pdipm->off_lambdai + i];
324: MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
325: }
326: }
328: /* (2) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
329: if (pdipm->Ng) {
330: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
331: for (i=0; i<pdipm->ng; i++) {
332: row = Jrstart + pdipm->off_lambdae + i;
334: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
335: proc = 0;
336: for (j=0; j < nc; j++) {
337: while (aj[j] >= cranges[proc+1]) proc++;
338: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
339: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
340: }
341: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
342: if (pdipm->kkt_pd) {
343: /* add shift \delta_c */
344: MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
345: }
346: }
347: }
349: /* (3) insert 3rd row block of Jpre: [ -grad h, 0, deltac, I] */
350: if (pdipm->Nh) {
351: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
352: for (i=0; i < pdipm->nh; i++) {
353: row = Jrstart + pdipm->off_lambdai + i;
354: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
355: proc = 0;
356: for (j=0; j < nc; j++) {
357: while (aj[j] >= cranges[proc+1]) proc++;
358: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
359: MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
360: }
361: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
362: if (pdipm->kkt_pd) {
363: /* add shift \delta_c */
364: MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
365: }
366: }
367: }
369: /* (4) insert 1st row block of Jpre: [Wxx, grad g', -grad h', 0] */
370: if (pdipm->Ng) { /* grad g' */
371: MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
372: }
373: if (pdipm->Nh) { /* grad h' */
374: MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
375: }
377: VecPlaceArray(pdipm->x,Xarr);
378: TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
379: VecResetArray(pdipm->x);
381: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
382: for (i=0; i<pdipm->nx; i++) {
383: row = Jrstart + i;
385: /* insert Wxx = fxx + ... -- provided by user */
386: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
387: proc = 0;
388: for (j=0; j < nc; j++) {
389: while (aj[j] >= cranges[proc+1]) proc++;
390: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
391: if (row == cols[0] && pdipm->kkt_pd) {
392: /* add shift deltaw to Wxx component */
393: MatSetValue(Jpre,row,cols[0],aa[j]+pdipm->deltaw,INSERT_VALUES);
394: } else {
395: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
396: }
397: }
398: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
400: /* insert grad g' */
401: if (pdipm->ng) {
402: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
403: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
404: proc = 0;
405: for (j=0; j < nc; j++) {
406: /* find row ownership of */
407: while (aj[j] >= ranges[proc+1]) proc++;
408: nx_all = rranges[proc+1] - rranges[proc];
409: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
410: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
411: }
412: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
413: }
415: /* insert -grad h' */
416: if (pdipm->nh) {
417: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
418: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
419: proc = 0;
420: for (j=0; j < nc; j++) {
421: /* find row ownership of */
422: while (aj[j] >= ranges[proc+1]) proc++;
423: nx_all = rranges[proc+1] - rranges[proc];
424: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
425: MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
426: }
427: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
428: }
429: }
430: VecRestoreArrayRead(X,&Xarr);
432: /* (6) assemble Jpre and J */
433: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
434: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
436: if (J != Jpre) {
437: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
438: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
439: }
440: return 0;
441: }
443: /*
444: TaoSnesFunction_PDIPM - Evaluate KKT function at X
446: Input Parameter:
447: snes - SNES context
448: X - KKT Vector
449: *ctx - pdipm
451: Output Parameter:
452: F - Updated Lagrangian vector
453: */
454: static PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
455: {
456: Tao tao=(Tao)ctx;
457: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
458: PetscScalar *Farr;
459: Vec x,L1;
460: PetscInt i;
461: const PetscScalar *Xarr,*carr,*zarr,*larr;
463: VecSet(F,0.0);
465: VecGetArrayRead(X,&Xarr);
466: VecGetArrayWrite(F,&Farr);
468: /* (0) Evaluate f, fx, gradG, gradH at X.x Note: pdipm->x is not changed below */
469: x = pdipm->x;
470: VecPlaceArray(x,Xarr);
471: TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);
473: /* Update ce, ci, and Jci at X.x */
474: TaoPDIPMUpdateConstraints(tao,x);
475: VecResetArray(x);
477: /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
478: L1 = pdipm->x;
479: VecPlaceArray(L1,Farr); /* L1 = 0.0 */
480: if (pdipm->Nci) {
481: if (pdipm->Nh) {
482: /* L1 += gradH'*DI. Note: tao->DI is not changed below */
483: VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
484: MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
485: VecResetArray(tao->DI);
486: }
488: /* L1 += Jci_xb'*lambdai_xb */
489: VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
490: MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
491: VecResetArray(pdipm->lambdai_xb);
493: /* L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
494: VecScale(L1,-1.0);
495: }
497: /* L1 += fx */
498: VecAXPY(L1,1.0,tao->gradient);
500: if (pdipm->Nce) {
501: if (pdipm->Ng) {
502: /* L1 += gradG'*DE. Note: tao->DE is not changed below */
503: VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
504: MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
505: VecResetArray(tao->DE);
506: }
507: if (pdipm->Nxfixed) {
508: /* L1 += Jce_xfixed'*lambdae_xfixed */
509: VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
510: MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
511: VecResetArray(pdipm->lambdae_xfixed);
512: }
513: }
514: VecResetArray(L1);
516: /* (2) L2 = ce(x) */
517: if (pdipm->Nce) {
518: VecGetArrayRead(pdipm->ce,&carr);
519: for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
520: VecRestoreArrayRead(pdipm->ce,&carr);
521: }
523: if (pdipm->Nci) {
524: if (pdipm->solve_symmetric_kkt) {
525: /* (3) L3 = z - ci(x);
526: (4) L4 = Lambdai * e - mu/z *e */
527: VecGetArrayRead(pdipm->ci,&carr);
528: larr = Xarr+pdipm->off_lambdai;
529: zarr = Xarr+pdipm->off_z;
530: for (i=0; i<pdipm->nci; i++) {
531: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
532: Farr[pdipm->off_z + i] = larr[i] - pdipm->mu/zarr[i];
533: }
534: VecRestoreArrayRead(pdipm->ci,&carr);
535: } else {
536: /* (3) L3 = z - ci(x);
537: (4) L4 = Z * Lambdai * e - mu * e */
538: VecGetArrayRead(pdipm->ci,&carr);
539: larr = Xarr+pdipm->off_lambdai;
540: zarr = Xarr+pdipm->off_z;
541: for (i=0; i<pdipm->nci; i++) {
542: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
543: Farr[pdipm->off_z + i] = zarr[i]*larr[i] - pdipm->mu;
544: }
545: VecRestoreArrayRead(pdipm->ci,&carr);
546: }
547: }
549: VecRestoreArrayRead(X,&Xarr);
550: VecRestoreArrayWrite(F,&Farr);
551: return 0;
552: }
554: /*
555: Evaluate F(X); then update update tao->gnorm0, tao->step = mu,
556: tao->residual = norm2(F_x,F_z) and tao->cnorm = norm2(F_ce,F_ci).
557: */
558: static PetscErrorCode TaoSNESFunction_PDIPM_residual(SNES snes,Vec X,Vec F,void *ctx)
559: {
560: Tao tao=(Tao)ctx;
561: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
562: PetscScalar *Farr,*tmparr;
563: Vec L1;
564: PetscInt i;
565: PetscReal res[2],cnorm[2];
566: const PetscScalar *Xarr=NULL;
568: TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);
569: VecGetArrayWrite(F,&Farr);
570: VecGetArrayRead(X,&Xarr);
572: /* compute res[0] = norm2(F_x) */
573: L1 = pdipm->x;
574: VecPlaceArray(L1,Farr);
575: VecNorm(L1,NORM_2,&res[0]);
576: VecResetArray(L1);
578: /* compute res[1] = norm2(F_z), cnorm[1] = norm2(F_ci) */
579: if (pdipm->z) {
580: if (pdipm->solve_symmetric_kkt) {
581: VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
582: if (pdipm->Nci) {
583: VecGetArrayWrite(pdipm->z,&tmparr);
584: for (i=0; i<pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
585: VecRestoreArrayWrite(pdipm->z,&tmparr);
586: }
588: VecNorm(pdipm->z,NORM_2,&res[1]);
590: if (pdipm->Nci) {
591: VecGetArrayWrite(pdipm->z,&tmparr);
592: for (i=0; i<pdipm->nci; i++) {
593: tmparr[i] /= Xarr[pdipm->off_z + i];
594: }
595: VecRestoreArrayWrite(pdipm->z,&tmparr);
596: }
597: VecResetArray(pdipm->z);
598: } else { /* !solve_symmetric_kkt */
599: VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
600: VecNorm(pdipm->z,NORM_2,&res[1]);
601: VecResetArray(pdipm->z);
602: }
604: VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
605: VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
606: VecResetArray(pdipm->ci);
607: } else {
608: res[1] = 0.0; cnorm[1] = 0.0;
609: }
611: /* compute cnorm[0] = norm2(F_ce) */
612: if (pdipm->Nce) {
613: VecPlaceArray(pdipm->ce,Farr+pdipm->off_lambdae);
614: VecNorm(pdipm->ce,NORM_2,&cnorm[0]);
615: VecResetArray(pdipm->ce);
616: } else cnorm[0] = 0.0;
618: VecRestoreArrayWrite(F,&Farr);
619: VecRestoreArrayRead(X,&Xarr);
621: tao->gnorm0 = tao->residual;
622: tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);
623: tao->cnorm = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);
624: tao->step = pdipm->mu;
625: return 0;
626: }
628: /*
629: KKTAddShifts - Check the inertia of Cholesky factor of KKT matrix.
630: If it does not match the numbers of prime and dual variables, add shifts to the KKT matrix.
631: */
632: static PetscErrorCode KKTAddShifts(Tao tao,SNES snes,Vec X)
633: {
634: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
635: KSP ksp;
636: PC pc;
637: Mat Factor;
638: PetscBool isCHOL;
639: PetscInt nneg,nzero,npos;
641: /* Get the inertia of Cholesky factor */
642: SNESGetKSP(snes,&ksp);
643: KSPGetPC(ksp,&pc);
644: PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);
645: if (!isCHOL) return 0;
647: PCFactorGetMatrix(pc,&Factor);
648: MatGetInertia(Factor,&nneg,&nzero,&npos);
650: if (npos < pdipm->Nx+pdipm->Nci) {
651: pdipm->deltaw = PetscMax(pdipm->lastdeltaw/3, 1.e-4*PETSC_MACHINE_EPSILON);
652: PetscInfo(tao,"Test reduced deltaw=%g; previous MatInertia: nneg %D, nzero %D, npos %D(<%D)\n",(double)pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
653: TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
654: PCSetUp(pc);
655: MatGetInertia(Factor,&nneg,&nzero,&npos);
657: if (npos < pdipm->Nx+pdipm->Nci) {
658: pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
659: while (npos < pdipm->Nx+pdipm->Nci && pdipm->deltaw <= 1./PETSC_SMALL) { /* increase deltaw */
660: PetscInfo(tao," deltaw=%g fails, MatInertia: nneg %D, nzero %D, npos %D(<%D)\n",(double)pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
661: pdipm->deltaw = PetscMin(8*pdipm->deltaw,PetscPowReal(10,20));
662: TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
663: PCSetUp(pc);
664: MatGetInertia(Factor,&nneg,&nzero,&npos);
665: }
669: PetscInfo(tao,"Updated deltaw %g\n",(double)pdipm->deltaw);
670: pdipm->lastdeltaw = pdipm->deltaw;
671: pdipm->deltaw = 0.0;
672: }
673: }
675: if (nzero) { /* Jacobian is singular */
676: if (pdipm->deltac == 0.0) {
677: pdipm->deltac = PETSC_SQRT_MACHINE_EPSILON;
678: } else {
679: pdipm->deltac = pdipm->deltac*PetscPowReal(pdipm->mu,.25);
680: }
681: PetscInfo(tao,"Updated deltac=%g, MatInertia: nneg %D, nzero %D(!=0), npos %D\n",(double)pdipm->deltac,nneg,nzero,npos);
682: TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
683: PCSetUp(pc);
684: MatGetInertia(Factor,&nneg,&nzero,&npos);
685: }
686: return 0;
687: }
689: /*
690: PCPreSolve_PDIPM -- called betwee MatFactorNumeric() and MatSolve()
691: */
692: PetscErrorCode PCPreSolve_PDIPM(PC pc,KSP ksp)
693: {
694: Tao tao;
695: TAO_PDIPM *pdipm;
697: KSPGetApplicationContext(ksp,&tao);
698: pdipm = (TAO_PDIPM*)tao->data;
699: KKTAddShifts(tao,pdipm->snes,pdipm->X);
700: return 0;
701: }
703: /*
704: SNESLineSearch_PDIPM - Custom line search used with PDIPM.
706: Collective on TAO
708: Notes:
709: This routine employs a simple backtracking line-search to keep
710: the slack variables (z) and inequality constraints Lagrange multipliers
711: (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
712: alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
713: slack variables are updated as X = X - alpha_d*dx. The constraint multipliers
714: are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
715: is also updated as mu = mu + z'lambdai/Nci
716: */
717: static PetscErrorCode SNESLineSearch_PDIPM(SNESLineSearch linesearch,void *ctx)
718: {
719: Tao tao=(Tao)ctx;
720: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
721: SNES snes;
722: Vec X,F,Y;
723: PetscInt i,iter;
724: PetscReal alpha_p=1.0,alpha_d=1.0,alpha[4];
725: PetscScalar *Xarr,*z,*lambdai,dot,*taosolarr;
726: const PetscScalar *dXarr,*dz,*dlambdai;
728: SNESLineSearchGetSNES(linesearch,&snes);
729: SNESGetIterationNumber(snes,&iter);
731: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
732: SNESLineSearchGetVecs(linesearch,&X,&F,&Y,NULL,NULL);
734: VecGetArrayWrite(X,&Xarr);
735: VecGetArrayRead(Y,&dXarr);
736: z = Xarr + pdipm->off_z;
737: dz = dXarr + pdipm->off_z;
738: for (i=0; i < pdipm->nci; i++) {
739: if (z[i] - dz[i] < 0.0) alpha_p = PetscMin(alpha_p, 0.9999*z[i]/dz[i]);
740: }
742: lambdai = Xarr + pdipm->off_lambdai;
743: dlambdai = dXarr + pdipm->off_lambdai;
745: for (i=0; i<pdipm->nci; i++) {
746: if (lambdai[i] - dlambdai[i] < 0.0) alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i], alpha_d);
747: }
749: alpha[0] = alpha_p;
750: alpha[1] = alpha_d;
751: VecRestoreArrayRead(Y,&dXarr);
752: VecRestoreArrayWrite(X,&Xarr);
754: /* alpha = min(alpha) over all processes */
755: MPI_Allreduce(alpha,alpha+2,2,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)tao));
757: alpha_p = alpha[2];
758: alpha_d = alpha[3];
760: /* X = X - alpha * Y */
761: VecGetArrayWrite(X,&Xarr);
762: VecGetArrayRead(Y,&dXarr);
763: for (i=0; i<pdipm->nx; i++) Xarr[i] -= alpha_p * dXarr[i];
764: for (i=0; i<pdipm->nce; i++) Xarr[i+pdipm->off_lambdae] -= alpha_d * dXarr[i+pdipm->off_lambdae];
766: for (i=0; i<pdipm->nci; i++) {
767: Xarr[i+pdipm->off_lambdai] -= alpha_d * dXarr[i+pdipm->off_lambdai];
768: Xarr[i+pdipm->off_z] -= alpha_p * dXarr[i+pdipm->off_z];
769: }
770: VecGetArrayWrite(tao->solution,&taosolarr);
771: PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
772: VecRestoreArrayWrite(tao->solution,&taosolarr);
774: VecRestoreArrayWrite(X,&Xarr);
775: VecRestoreArrayRead(Y,&dXarr);
777: /* Update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
778: if (pdipm->z) {
779: VecDot(pdipm->z,pdipm->lambdai,&dot);
780: } else dot = 0.0;
782: /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu) */
783: pdipm->mu = pdipm->mu_update_factor * dot/pdipm->Nci;
785: /* Update F; get tao->residual and tao->cnorm */
786: TaoSNESFunction_PDIPM_residual(snes,X,F,(void*)tao);
788: tao->niter++;
789: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
790: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
792: (*tao->ops->convergencetest)(tao,tao->cnvP);
793: if (tao->reason) {
794: SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
795: }
796: return 0;
797: }
799: /*
800: TaoSolve_PDIPM
802: Input Parameter:
803: tao - TAO context
805: Output Parameter:
806: tao - TAO context
807: */
808: PetscErrorCode TaoSolve_PDIPM(Tao tao)
809: {
810: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
811: SNESLineSearch linesearch; /* SNESLineSearch context */
812: Vec dummy;
816: /* Initialize all variables */
817: TaoPDIPMInitializeSolution(tao);
819: /* Set linesearch */
820: SNESGetLineSearch(pdipm->snes,&linesearch);
821: SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
822: SNESLineSearchShellSetUserFunc(linesearch,SNESLineSearch_PDIPM,tao);
823: SNESLineSearchSetFromOptions(linesearch);
825: tao->reason = TAO_CONTINUE_ITERATING;
827: /* -tao_monitor for iteration 0 and check convergence */
828: VecDuplicate(pdipm->X,&dummy);
829: TaoSNESFunction_PDIPM_residual(pdipm->snes,pdipm->X,dummy,(void*)tao);
831: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
832: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
833: VecDestroy(&dummy);
834: (*tao->ops->convergencetest)(tao,tao->cnvP);
835: if (tao->reason) {
836: SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
837: }
839: while (tao->reason == TAO_CONTINUE_ITERATING) {
840: SNESConvergedReason reason;
841: SNESSolve(pdipm->snes,NULL,pdipm->X);
843: /* Check SNES convergence */
844: SNESGetConvergedReason(pdipm->snes,&reason);
845: if (reason < 0) {
846: PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
847: }
849: /* Check TAO convergence */
851: }
852: return 0;
853: }
855: /*
856: TaoView_PDIPM - View PDIPM
858: Input Parameter:
859: tao - TAO object
860: viewer - PetscViewer
862: Output:
863: */
864: PetscErrorCode TaoView_PDIPM(Tao tao,PetscViewer viewer)
865: {
866: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
868: tao->constrained = PETSC_TRUE;
869: PetscViewerASCIIPushTab(viewer);
870: PetscViewerASCIIPrintf(viewer,"Number of prime=%D, Number of dual=%D\n",pdipm->Nx+pdipm->Nci,pdipm->Nce + pdipm->Nci);
871: if (pdipm->kkt_pd) {
872: PetscViewerASCIIPrintf(viewer,"KKT shifts deltaw=%g, deltac=%g\n",(double)pdipm->deltaw,(double)pdipm->deltac);
873: }
874: PetscViewerASCIIPopTab(viewer);
875: return 0;
876: }
878: /*
879: TaoSetup_PDIPM - Sets up tao and pdipm
881: Input Parameter:
882: tao - TAO object
884: Output: pdipm - initialized object
885: */
886: PetscErrorCode TaoSetup_PDIPM(Tao tao)
887: {
888: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
889: PetscErrorCode ierr;
890: MPI_Comm comm;
891: PetscMPIInt size;
892: PetscInt row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
893: PetscInt offset,*xa,*xb,i,j,rstart,rend;
894: PetscScalar one=1.0,neg_one=-1.0;
895: const PetscInt *cols,*rranges,*cranges,*aj,*ranges;
896: const PetscScalar *aa,*Xarr;
897: Mat J,jac_equality_trans,jac_inequality_trans;
898: Mat Jce_xfixed_trans,Jci_xb_trans;
899: PetscInt *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];
901: PetscObjectGetComm((PetscObject)tao,&comm);
902: MPI_Comm_size(comm,&size);
904: /* (1) Setup Bounds and create Tao vectors */
905: TaoPDIPMSetUpBounds(tao);
907: if (!tao->gradient) {
908: VecDuplicate(tao->solution,&tao->gradient);
909: VecDuplicate(tao->solution,&tao->stepdirection);
910: }
912: /* (2) Get sizes */
913: /* Size of vector x - This is set by TaoSetSolution */
914: VecGetSize(tao->solution,&pdipm->Nx);
915: VecGetLocalSize(tao->solution,&pdipm->nx);
917: /* Size of equality constraints and vectors */
918: if (tao->constraints_equality) {
919: VecGetSize(tao->constraints_equality,&pdipm->Ng);
920: VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
921: } else {
922: pdipm->ng = pdipm->Ng = 0;
923: }
925: pdipm->nce = pdipm->ng + pdipm->nxfixed;
926: pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;
928: /* Size of inequality constraints and vectors */
929: if (tao->constraints_inequality) {
930: VecGetSize(tao->constraints_inequality,&pdipm->Nh);
931: VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
932: } else {
933: pdipm->nh = pdipm->Nh = 0;
934: }
936: pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2*pdipm->nxbox;
937: pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2*pdipm->Nxbox;
939: /* Full size of the KKT system to be solved */
940: pdipm->n = pdipm->nx + pdipm->nce + 2*pdipm->nci;
941: pdipm->N = pdipm->Nx + pdipm->Nce + 2*pdipm->Nci;
943: /* (3) Offsets for subvectors */
944: pdipm->off_lambdae = pdipm->nx;
945: pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
946: pdipm->off_z = pdipm->off_lambdai + pdipm->nci;
948: /* (4) Create vectors and subvectors */
949: /* Ce and Ci vectors */
950: VecCreate(comm,&pdipm->ce);
951: VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
952: VecSetFromOptions(pdipm->ce);
954: VecCreate(comm,&pdipm->ci);
955: VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
956: VecSetFromOptions(pdipm->ci);
958: /* X=[x; lambdae; lambdai; z] for the big KKT system */
959: VecCreate(comm,&pdipm->X);
960: VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
961: VecSetFromOptions(pdipm->X);
963: /* Subvectors; they share local arrays with X */
964: VecGetArrayRead(pdipm->X,&Xarr);
965: /* x shares local array with X.x */
966: if (pdipm->Nx) {
967: VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
968: }
970: /* lambdae shares local array with X.lambdae */
971: if (pdipm->Nce) {
972: VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
973: }
975: /* tao->DE shares local array with X.lambdae_g */
976: if (pdipm->Ng) {
977: VecCreateMPIWithArray(comm,1,pdipm->ng,pdipm->Ng,Xarr+pdipm->off_lambdae,&tao->DE);
979: VecCreate(comm,&pdipm->lambdae_xfixed);
980: VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
981: VecSetFromOptions(pdipm->lambdae_xfixed);
982: }
984: if (pdipm->Nci) {
985: /* lambdai shares local array with X.lambdai */
986: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_lambdai,&pdipm->lambdai);
988: /* z for slack variables; it shares local array with X.z */
989: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_z,&pdipm->z);
990: }
992: /* tao->DI which shares local array with X.lambdai_h */
993: if (pdipm->Nh) {
994: VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
995: }
996: VecCreate(comm,&pdipm->lambdai_xb);
997: VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
998: VecSetFromOptions(pdipm->lambdai_xb);
1000: VecRestoreArrayRead(pdipm->X,&Xarr);
1002: /* (5) Create Jacobians Jce_xfixed and Jci */
1003: /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
1004: if (pdipm->Nxfixed) {
1005: /* Create Jce_xfixed */
1006: MatCreate(comm,&pdipm->Jce_xfixed);
1007: MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1008: MatSetFromOptions(pdipm->Jce_xfixed);
1009: MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
1010: MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);
1012: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
1013: ISGetIndices(pdipm->isxfixed,&cols);
1014: k = 0;
1015: for (row = Jcrstart; row < Jcrend; row++) {
1016: MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
1017: k++;
1018: }
1019: ISRestoreIndices(pdipm->isxfixed, &cols);
1020: MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1021: MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1022: }
1024: /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
1025: MatCreate(comm,&pdipm->Jci_xb);
1026: MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1027: MatSetFromOptions(pdipm->Jci_xb);
1028: MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
1029: MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);
1031: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
1032: offset = Jcrstart;
1033: if (pdipm->Nxub) {
1034: /* Add xub to Jci_xb */
1035: ISGetIndices(pdipm->isxub,&cols);
1036: k = 0;
1037: for (row = offset; row < offset + pdipm->nxub; row++) {
1038: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1039: k++;
1040: }
1041: ISRestoreIndices(pdipm->isxub, &cols);
1042: }
1044: if (pdipm->Nxlb) {
1045: /* Add xlb to Jci_xb */
1046: ISGetIndices(pdipm->isxlb,&cols);
1047: k = 0;
1048: offset += pdipm->nxub;
1049: for (row = offset; row < offset + pdipm->nxlb; row++) {
1050: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
1051: k++;
1052: }
1053: ISRestoreIndices(pdipm->isxlb, &cols);
1054: }
1056: /* Add xbox to Jci_xb */
1057: if (pdipm->Nxbox) {
1058: ISGetIndices(pdipm->isxbox,&cols);
1059: k = 0;
1060: offset += pdipm->nxlb;
1061: for (row = offset; row < offset + pdipm->nxbox; row++) {
1062: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1063: tmp = row + pdipm->nxbox;
1064: MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
1065: k++;
1066: }
1067: ISRestoreIndices(pdipm->isxbox, &cols);
1068: }
1070: MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1071: MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1072: /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */
1074: /* (6) Set up ISs for PC Fieldsplit */
1075: if (pdipm->solve_reduced_kkt) {
1076: PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
1077: for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1078: for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;
1080: ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
1081: ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
1082: }
1084: /* (7) Gather offsets from all processes */
1085: PetscMalloc1(size,&pdipm->nce_all);
1087: /* Get rstart of KKT matrix */
1088: MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
1089: rstart -= pdipm->n;
1091: MPI_Allgather(&pdipm->nce,1,MPIU_INT,pdipm->nce_all,1,MPIU_INT,comm);
1093: PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
1094: MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
1095: MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
1096: MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);
1098: MatGetOwnershipRanges(tao->hessian,&rranges);
1099: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
1101: if (pdipm->Ng) {
1102: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
1103: MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
1104: }
1105: if (pdipm->Nh) {
1106: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
1107: MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
1108: }
1110: /* Count dnz,onz for preallocation of KKT matrix */
1111: jac_equality_trans = pdipm->jac_equality_trans;
1112: jac_inequality_trans = pdipm->jac_inequality_trans;
1113: nce_all = pdipm->nce_all;
1115: if (pdipm->Nxfixed) {
1116: MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
1117: }
1118: MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);
1120: MatPreallocateInitialize(comm,pdipm->n,pdipm->n,dnz,onz);
1122: /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
1123: TaoPDIPMEvaluateFunctionsAndJacobians(tao,pdipm->x);
1124: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
1126: /* Insert tao->hessian */
1127: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
1128: for (i=0; i<pdipm->nx; i++) {
1129: row = rstart + i;
1131: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1132: proc = 0;
1133: for (j=0; j < nc; j++) {
1134: while (aj[j] >= cranges[proc+1]) proc++;
1135: col = aj[j] - cranges[proc] + Jranges[proc];
1136: MatPreallocateSet(row,1,&col,dnz,onz);
1137: }
1138: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1140: if (pdipm->ng) {
1141: /* Insert grad g' */
1142: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1143: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
1144: proc = 0;
1145: for (j=0; j < nc; j++) {
1146: /* find row ownership of */
1147: while (aj[j] >= ranges[proc+1]) proc++;
1148: nx_all = rranges[proc+1] - rranges[proc];
1149: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1150: MatPreallocateSet(row,1,&col,dnz,onz);
1151: }
1152: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1153: }
1155: /* Insert Jce_xfixed^T' */
1156: if (pdipm->nxfixed) {
1157: MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1158: MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1159: proc = 0;
1160: for (j=0; j < nc; j++) {
1161: /* find row ownership of */
1162: while (aj[j] >= ranges[proc+1]) proc++;
1163: nx_all = rranges[proc+1] - rranges[proc];
1164: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1165: MatPreallocateSet(row,1,&col,dnz,onz);
1166: }
1167: MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1168: }
1170: if (pdipm->nh) {
1171: /* Insert -grad h' */
1172: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1173: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1174: proc = 0;
1175: for (j=0; j < nc; j++) {
1176: /* find row ownership of */
1177: while (aj[j] >= ranges[proc+1]) proc++;
1178: nx_all = rranges[proc+1] - rranges[proc];
1179: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1180: MatPreallocateSet(row,1,&col,dnz,onz);
1181: }
1182: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1183: }
1185: /* Insert Jci_xb^T' */
1186: MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1187: MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1188: proc = 0;
1189: for (j=0; j < nc; j++) {
1190: /* find row ownership of */
1191: while (aj[j] >= ranges[proc+1]) proc++;
1192: nx_all = rranges[proc+1] - rranges[proc];
1193: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1194: MatPreallocateSet(row,1,&col,dnz,onz);
1195: }
1196: MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1197: }
1199: /* 2nd Row block of KKT matrix: [grad Ce, deltac*I, 0, 0] */
1200: if (pdipm->Ng) {
1201: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
1202: for (i=0; i < pdipm->ng; i++) {
1203: row = rstart + pdipm->off_lambdae + i;
1205: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1206: proc = 0;
1207: for (j=0; j < nc; j++) {
1208: while (aj[j] >= cranges[proc+1]) proc++;
1209: col = aj[j] - cranges[proc] + Jranges[proc];
1210: MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1211: }
1212: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1213: }
1214: }
1215: /* Jce_xfixed */
1216: if (pdipm->Nxfixed) {
1217: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1218: for (i=0; i < (pdipm->nce - pdipm->ng); i++) {
1219: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1221: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1224: proc = 0;
1225: j = 0;
1226: while (cols[j] >= cranges[proc+1]) proc++;
1227: col = cols[j] - cranges[proc] + Jranges[proc];
1228: MatPreallocateSet(row,1,&col,dnz,onz);
1229: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1230: }
1231: }
1233: /* 3rd Row block of KKT matrix: [ gradCi, 0, deltac*I, -I] */
1234: if (pdipm->Nh) {
1235: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
1236: for (i=0; i < pdipm->nh; i++) {
1237: row = rstart + pdipm->off_lambdai + i;
1239: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1240: proc = 0;
1241: for (j=0; j < nc; j++) {
1242: while (aj[j] >= cranges[proc+1]) proc++;
1243: col = aj[j] - cranges[proc] + Jranges[proc];
1244: MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1245: }
1246: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1247: }
1248: /* I */
1249: for (i=0; i < pdipm->nh; i++) {
1250: row = rstart + pdipm->off_lambdai + i;
1251: col = rstart + pdipm->off_z + i;
1252: MatPreallocateSet(row,1,&col,dnz,onz);
1253: }
1254: }
1256: /* Jci_xb */
1257: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1258: for (i=0; i < (pdipm->nci - pdipm->nh); i++) {
1259: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1261: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1263: proc = 0;
1264: for (j=0; j < nc; j++) {
1265: while (cols[j] >= cranges[proc+1]) proc++;
1266: col = cols[j] - cranges[proc] + Jranges[proc];
1267: MatPreallocateSet(row,1,&col,dnz,onz);
1268: }
1269: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1270: /* I */
1271: col = rstart + pdipm->off_z + pdipm->nh + i;
1272: MatPreallocateSet(row,1,&col,dnz,onz);
1273: }
1275: /* 4-th Row block of KKT matrix: Z and Ci */
1276: for (i=0; i < pdipm->nci; i++) {
1277: row = rstart + pdipm->off_z + i;
1278: cols1[0] = rstart + pdipm->off_lambdai + i;
1279: cols1[1] = row;
1280: MatPreallocateSet(row,2,cols1,dnz,onz);
1281: }
1283: /* diagonal entry */
1284: for (i=0; i<pdipm->n; i++) dnz[i]++; /* diagonal entry */
1286: /* Create KKT matrix */
1287: MatCreate(comm,&J);
1288: MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1289: MatSetFromOptions(J);
1290: MatSeqAIJSetPreallocation(J,0,dnz);
1291: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1292: MatPreallocateFinalize(dnz,onz);
1293: pdipm->K = J;
1295: /* (8) Insert constant entries to K */
1296: /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1297: MatGetOwnershipRange(J,&rstart,&rend);
1298: for (i=rstart; i<rend; i++) {
1299: MatSetValue(J,i,i,0.0,INSERT_VALUES);
1300: }
1301: /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1302: if (pdipm->kkt_pd) {
1303: for (i=0; i<pdipm->nh; i++) {
1304: row = rstart + i;
1305: MatSetValue(J,row,row,pdipm->deltaw,INSERT_VALUES);
1306: }
1307: }
1309: /* Row block of K: [ grad Ce, 0, 0, 0] */
1310: if (pdipm->Nxfixed) {
1311: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1312: for (i=0; i < (pdipm->nce - pdipm->ng); i++) {
1313: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1315: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1316: proc = 0;
1317: for (j=0; j < nc; j++) {
1318: while (cols[j] >= cranges[proc+1]) proc++;
1319: col = cols[j] - cranges[proc] + Jranges[proc];
1320: MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1321: MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1322: }
1323: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1324: }
1325: }
1327: /* Row block of K: [ -grad Ci, 0, 0, I] */
1328: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1329: for (i=0; i < pdipm->nci - pdipm->nh; i++) {
1330: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1332: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1333: proc = 0;
1334: for (j=0; j < nc; j++) {
1335: while (cols[j] >= cranges[proc+1]) proc++;
1336: col = cols[j] - cranges[proc] + Jranges[proc];
1337: MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1338: MatSetValue(J,row,col,-aa[j],INSERT_VALUES);
1339: }
1340: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1342: col = rstart + pdipm->off_z + pdipm->nh + i;
1343: MatSetValue(J,row,col,1,INSERT_VALUES);
1344: }
1346: for (i=0; i < pdipm->nh; i++) {
1347: row = rstart + pdipm->off_lambdai + i;
1348: col = rstart + pdipm->off_z + i;
1349: MatSetValue(J,row,col,1,INSERT_VALUES);
1350: }
1352: /* Row block of K: [ 0, 0, I, ...] */
1353: for (i=0; i < pdipm->nci; i++) {
1354: row = rstart + pdipm->off_z + i;
1355: col = rstart + pdipm->off_lambdai + i;
1356: MatSetValue(J,row,col,1,INSERT_VALUES);
1357: }
1359: if (pdipm->Nxfixed) {
1360: MatDestroy(&Jce_xfixed_trans);
1361: }
1362: MatDestroy(&Jci_xb_trans);
1363: PetscFree3(ng_all,nh_all,Jranges);
1365: /* (9) Set up nonlinear solver SNES */
1366: SNESSetFunction(pdipm->snes,NULL,TaoSNESFunction_PDIPM,(void*)tao);
1367: SNESSetJacobian(pdipm->snes,J,J,TaoSNESJacobian_PDIPM,(void*)tao);
1369: if (pdipm->solve_reduced_kkt) {
1370: PC pc;
1371: KSPGetPC(tao->ksp,&pc);
1372: PCSetType(pc,PCFIELDSPLIT);
1373: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1374: PCFieldSplitSetIS(pc,"2",pdipm->is2);
1375: PCFieldSplitSetIS(pc,"1",pdipm->is1);
1376: }
1377: SNESSetFromOptions(pdipm->snes);
1379: /* (10) Setup PCPreSolve() for pdipm->solve_symmetric_kkt */
1380: if (pdipm->solve_symmetric_kkt) {
1381: KSP ksp;
1382: PC pc;
1383: PetscBool isCHOL;
1384: SNESGetKSP(pdipm->snes,&ksp);
1385: KSPGetPC(ksp,&pc);
1386: PCSetPreSolve(pc,PCPreSolve_PDIPM);
1388: PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);
1389: if (isCHOL) {
1390: Mat Factor;
1391: PetscBool isMUMPS;
1392: PCFactorGetMatrix(pc,&Factor);
1393: PetscObjectTypeCompare((PetscObject)Factor,"mumps",&isMUMPS);
1394: if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
1395: #if defined(PETSC_HAVE_MUMPS)
1396: MatMumpsSetIcntl(Factor,24,1); /* detection of null pivot rows */
1397: if (size > 1) {
1398: MatMumpsSetIcntl(Factor,13,1); /* parallelism of the root node (enable ScaLAPACK) and its splitting */
1399: }
1400: #else
1401: SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires external package MUMPS");
1402: #endif
1403: }
1404: }
1405: }
1406: return 0;
1407: }
1409: /*
1410: TaoDestroy_PDIPM - Destroys the pdipm object
1412: Input:
1413: full pdipm
1415: Output:
1416: Destroyed pdipm
1417: */
1418: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1419: {
1420: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1422: /* Freeing Vectors assocaiated with KKT (X) */
1423: VecDestroy(&pdipm->x); /* Solution x */
1424: VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1425: VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1426: VecDestroy(&pdipm->z); /* Slack variables */
1427: VecDestroy(&pdipm->X); /* Big KKT system vector [x; lambdae; lambdai; z] */
1429: /* work vectors */
1430: VecDestroy(&pdipm->lambdae_xfixed);
1431: VecDestroy(&pdipm->lambdai_xb);
1433: /* Legrangian equality and inequality Vec */
1434: VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1435: VecDestroy(&pdipm->ci); /* Vec of inequality constraints */
1437: /* Matrices */
1438: MatDestroy(&pdipm->Jce_xfixed);
1439: MatDestroy(&pdipm->Jci_xb)); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb; J(nxbx)] */
1440: MatDestroy(&pdipm->K);
1442: /* Index Sets */
1443: if (pdipm->Nxub) {
1444: ISDestroy(&pdipm->isxub); /* Finite upper bound only -inf < x < ub */
1445: }
1447: if (pdipm->Nxlb) {
1448: ISDestroy(&pdipm->isxlb); /* Finite lower bound only lb <= x < inf */
1449: }
1451: if (pdipm->Nxfixed) {
1452: ISDestroy(&pdipm->isxfixed); /* Fixed variables lb = x = ub */
1453: }
1455: if (pdipm->Nxbox) {
1456: ISDestroy(&pdipm->isxbox); /* Boxed variables lb <= x <= ub */
1457: }
1459: if (pdipm->Nxfree) {
1460: ISDestroy(&pdipm->isxfree); /* Free variables -inf <= x <= inf */
1461: }
1463: if (pdipm->solve_reduced_kkt) {
1464: ISDestroy(&pdipm->is1);
1465: ISDestroy(&pdipm->is2);
1466: }
1468: /* SNES */
1469: SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1470: PetscFree(pdipm->nce_all);
1471: MatDestroy(&pdipm->jac_equality_trans);
1472: MatDestroy(&pdipm->jac_inequality_trans);
1474: /* Destroy pdipm */
1475: PetscFree(tao->data); /* Holding locations of pdipm */
1477: /* Destroy Dual */
1478: VecDestroy(&tao->DE); /* equality dual */
1479: VecDestroy(&tao->DI); /* dinequality dual */
1480: return 0;
1481: }
1483: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1484: {
1485: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1487: PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1488: PetscOptionsReal("-tao_pdipm_push_init_slack","parameter to push initial slack variables away from bounds",NULL,pdipm->push_init_slack,&pdipm->push_init_slack,NULL);
1489: PetscOptionsReal("-tao_pdipm_push_init_lambdai","parameter to push initial (inequality) dual variables away from bounds",NULL,pdipm->push_init_lambdai,&pdipm->push_init_lambdai,NULL);
1490: PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1491: PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1492: PetscOptionsBool("-tao_pdipm_symmetric_kkt","Solve non reduced symmetric KKT system",NULL,pdipm->solve_symmetric_kkt,&pdipm->solve_symmetric_kkt,NULL);
1493: PetscOptionsBool("-tao_pdipm_kkt_shift_pd","Add shifts to make KKT matrix positive definite",NULL,pdipm->kkt_pd,&pdipm->kkt_pd,NULL);
1494: PetscOptionsTail();
1495: return 0;
1496: }
1498: /*MC
1499: TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.
1501: Option Database Keys:
1502: + -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1503: . -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1504: . -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1505: . -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1506: - -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite
1508: Level: beginner
1509: M*/
1510: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1511: {
1512: TAO_PDIPM *pdipm;
1514: tao->ops->setup = TaoSetup_PDIPM;
1515: tao->ops->solve = TaoSolve_PDIPM;
1516: tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1517: tao->ops->view = TaoView_PDIPM;
1518: tao->ops->destroy = TaoDestroy_PDIPM;
1520: PetscNewLog(tao,&pdipm);
1521: tao->data = (void*)pdipm;
1523: pdipm->nx = pdipm->Nx = 0;
1524: pdipm->nxfixed = pdipm->Nxfixed = 0;
1525: pdipm->nxlb = pdipm->Nxlb = 0;
1526: pdipm->nxub = pdipm->Nxub = 0;
1527: pdipm->nxbox = pdipm->Nxbox = 0;
1528: pdipm->nxfree = pdipm->Nxfree = 0;
1530: pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1531: pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1532: pdipm->n = pdipm->N = 0;
1533: pdipm->mu = 1.0;
1534: pdipm->mu_update_factor = 0.1;
1536: pdipm->deltaw = 0.0;
1537: pdipm->lastdeltaw = 3*1.e-4;
1538: pdipm->deltac = 0.0;
1539: pdipm->kkt_pd = PETSC_FALSE;
1541: pdipm->push_init_slack = 1.0;
1542: pdipm->push_init_lambdai = 1.0;
1543: pdipm->solve_reduced_kkt = PETSC_FALSE;
1544: pdipm->solve_symmetric_kkt = PETSC_TRUE;
1546: /* Override default settings (unless already changed) */
1547: if (!tao->max_it_changed) tao->max_it = 200;
1548: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1550: SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1551: SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1552: SNESGetKSP(pdipm->snes,&tao->ksp);
1553: PetscObjectReference((PetscObject)tao->ksp);
1554: KSPSetApplicationContext(tao->ksp,(void *)tao);
1555: return 0;
1556: }