Actual source code: pdipm.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/constrained/impls/ipm/pdipm.h>
3: #include <petscsnes.h>
5: /*
6: TaoPDIPMEvaluateFunctionsAndJacobians - Evaluate the objective function f, gradient fx, constraints, and all the Jacobians at current vector
8: Collective on tao
10: Input Parameter:
11: + tao - solver context
12: - x - vector at which all objects to be evaluated
14: Level: beginner
16: .seealso: TaoPDIPMUpdateConstraints(), TaoPDIPMSetUpBounds()
17: */
18: PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)
19: {
21: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
24: /* Compute user objective function and gradient */
25: TaoComputeObjectiveAndGradient(tao,x,&pdipm->obj,tao->gradient);
27: /* Equality constraints and Jacobian */
28: if (pdipm->Ng) {
29: TaoComputeEqualityConstraints(tao,x,tao->constraints_equality);
30: TaoComputeJacobianEquality(tao,x,tao->jacobian_equality,tao->jacobian_equality_pre);
31: }
33: /* Inequality constraints and Jacobian */
34: if (pdipm->Nh) {
35: TaoComputeInequalityConstraints(tao,x,tao->constraints_inequality);
36: TaoComputeJacobianInequality(tao,x,tao->jacobian_inequality,tao->jacobian_inequality_pre);
37: }
38: return(0);
39: }
41: /*
42: TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x
44: Collective
46: Input Parameter:
47: + tao - Tao context
48: - x - vector at which constraints to be evaluted
50: Level: beginner
52: .seealso: TaoPDIPMEvaluateFunctionsAndJacobians()
53: */
54: PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao,Vec x)
55: {
56: PetscErrorCode ierr;
57: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
58: PetscInt i,offset,offset1,k,xstart;
59: PetscScalar *carr;
60: const PetscInt *ubptr,*lbptr,*bxptr,*fxptr;
61: const PetscScalar *xarr,*xuarr,*xlarr,*garr,*harr;
64: VecGetOwnershipRange(x,&xstart,NULL);
66: VecGetArrayRead(x,&xarr);
67: VecGetArrayRead(tao->XU,&xuarr);
68: VecGetArrayRead(tao->XL,&xlarr);
70: /* (1) Update ce vector */
71: VecGetArray(pdipm->ce,&carr);
73: if (pdipm->Ng) {
74: /* (1.a) Inserting updated g(x) */
75: VecGetArrayRead(tao->constraints_equality,&garr);
76: PetscMemcpy(carr,garr,pdipm->ng*sizeof(PetscScalar));
77: VecRestoreArrayRead(tao->constraints_equality,&garr);
78: }
80: /* (1.b) Update xfixed */
81: if (pdipm->Nxfixed) {
82: offset = pdipm->ng;
83: ISGetIndices(pdipm->isxfixed,&fxptr); /* global indices in x */
84: for (k=0;k < pdipm->nxfixed;k++){
85: i = fxptr[k]-xstart;
86: carr[offset + k] = xarr[i] - xuarr[i];
87: }
88: }
89: VecRestoreArray(pdipm->ce,&carr);
91: /* (2) Update ci vector */
92: VecGetArray(pdipm->ci,&carr);
94: if (pdipm->Nh) {
95: /* (2.a) Inserting updated h(x) */
96: VecGetArrayRead(tao->constraints_inequality,&harr);
97: PetscMemcpy(carr,harr,pdipm->nh*sizeof(PetscScalar));
98: VecRestoreArrayRead(tao->constraints_inequality,&harr);
99: }
101: /* (2.b) Update xub */
102: offset = pdipm->nh;
103: if (pdipm->Nxub) {
104: ISGetIndices(pdipm->isxub,&ubptr);
105: for (k=0; k<pdipm->nxub; k++){
106: i = ubptr[k]-xstart;
107: carr[offset + k] = xuarr[i] - xarr[i];
108: }
109: }
111: if (pdipm->Nxlb) {
112: /* (2.c) Update xlb */
113: offset += pdipm->nxub;
114: ISGetIndices(pdipm->isxlb,&lbptr); /* global indices in x */
115: for (k=0; k<pdipm->nxlb; k++){
116: i = lbptr[k]-xstart;
117: carr[offset + k] = xarr[i] - xlarr[i];
118: }
119: }
121: if (pdipm->Nxbox) {
122: /* (2.d) Update xbox */
123: offset += pdipm->nxlb;
124: offset1 = offset + pdipm->nxbox;
125: ISGetIndices(pdipm->isxbox,&bxptr); /* global indices in x */
126: for (k=0; k<pdipm->nxbox; k++){
127: i = bxptr[k]-xstart; /* local indices in x */
128: carr[offset+k] = xuarr[i] - xarr[i];
129: carr[offset1+k] = xarr[i] - xlarr[i];
130: }
131: }
132: VecRestoreArray(pdipm->ci,&carr);
134: /* Restoring Vectors */
135: VecRestoreArrayRead(x,&xarr);
136: VecRestoreArrayRead(tao->XU,&xuarr);
137: VecRestoreArrayRead(tao->XL,&xlarr);
138: return(0);
139: }
141: /*
142: TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x
144: Collective
146: Input Parameter:
147: . tao - holds pdipm and XL & XU
149: Level: beginner
151: .seealso: TaoPDIPMUpdateConstraints
152: */
153: PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
154: {
155: PetscErrorCode ierr;
156: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
157: const PetscScalar *xl,*xu;
158: PetscInt n,*ixlb,*ixub,*ixfixed,*ixfree,*ixbox,i,low,high,idx;
159: MPI_Comm comm;
160: PetscInt sendbuf[5],recvbuf[5];
163: /* Creates upper and lower bounds vectors on x, if not created already */
164: TaoComputeVariableBounds(tao);
166: VecGetLocalSize(tao->XL,&n);
167: PetscMalloc5(n,&ixlb,n,&ixub,n,&ixfree,n,&ixfixed,n,&ixbox);
169: VecGetOwnershipRange(tao->XL,&low,&high);
170: VecGetArrayRead(tao->XL,&xl);
171: VecGetArrayRead(tao->XU,&xu);
172: for (i=0; i<n; i++) {
173: idx = low + i;
174: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
175: if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
176: ixfixed[pdipm->nxfixed++] = idx;
177: } else ixbox[pdipm->nxbox++] = idx;
178: } else {
179: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
180: ixlb[pdipm->nxlb++] = idx;
181: } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
182: ixub[pdipm->nxlb++] = idx;
183: } else ixfree[pdipm->nxfree++] = idx;
184: }
185: }
186: VecRestoreArrayRead(tao->XL,&xl);
187: VecRestoreArrayRead(tao->XU,&xu);
189: PetscObjectGetComm((PetscObject)tao,&comm);
190: sendbuf[0] = pdipm->nxlb;
191: sendbuf[1] = pdipm->nxub;
192: sendbuf[2] = pdipm->nxfixed;
193: sendbuf[3] = pdipm->nxbox;
194: sendbuf[4] = pdipm->nxfree;
196: MPI_Allreduce(sendbuf,recvbuf,5,MPIU_INT,MPI_SUM,comm);
197: pdipm->Nxlb = recvbuf[0];
198: pdipm->Nxub = recvbuf[1];
199: pdipm->Nxfixed = recvbuf[2];
200: pdipm->Nxbox = recvbuf[3];
201: pdipm->Nxfree = recvbuf[4];
203: if (pdipm->Nxlb) {
204: ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);
205: }
206: if (pdipm->Nxub) {
207: ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);
208: }
209: if (pdipm->Nxfixed) {
210: ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);
211: }
212: if (pdipm->Nxbox) {
213: ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);
214: }
215: if (pdipm->Nxfree) {
216: ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);
217: }
218: PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);
219: return(0);
220: }
222: /*
223: TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
224: X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
225: four subvectors need to be initialized and its values copied over to X. Instead
226: of copying, we use VecPlace/ResetArray functions to share the memory locations for
227: X and the subvectors
229: Collective
231: Input Parameter:
232: . tao - Tao context
234: Level: beginner
235: */
236: PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
237: {
238: PetscErrorCode ierr;
239: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
240: PetscScalar *Xarr,*z,*lambdai;
241: PetscInt i;
242: const PetscScalar *xarr,*h;
245: VecGetArray(pdipm->X,&Xarr);
247: /* Set Initialize X.x = tao->solution */
248: VecGetArrayRead(tao->solution,&xarr);
249: PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));
250: VecRestoreArrayRead(tao->solution,&xarr);
252: /* Initialize X.lambdae = 0.0 */
253: if (pdipm->lambdae) {
254: VecSet(pdipm->lambdae,0.0);
255: }
256: /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
257: if (pdipm->lambdai) {
258: VecSet(pdipm->lambdai,pdipm->push_init_lambdai);
259: }
260: if (pdipm->z) {
261: VecSet(pdipm->z,pdipm->push_init_slack);
262: }
264: /* Additional modification for X.lambdai and X.z */
265: if (pdipm->lambdai) {
266: VecGetArray(pdipm->lambdai,&lambdai);
267: }
268: if (pdipm->z) {
269: VecGetArray(pdipm->z,&z);
270: }
271: if (pdipm->Nh) {
272: VecGetArrayRead(tao->constraints_inequality,&h);
273: for (i=0; i < pdipm->nh; i++) {
274: if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
275: if (pdipm->mu/z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu/z[i];
276: }
277: VecRestoreArrayRead(tao->constraints_inequality,&h);
278: }
279: if (pdipm->lambdai) {
280: VecRestoreArray(pdipm->lambdai,&lambdai);
281: }
282: if (pdipm->z) {
283: VecRestoreArray(pdipm->z,&z);
284: }
286: VecRestoreArray(pdipm->X,&Xarr);
287: return(0);
288: }
290: /*
291: TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X
293: Input Parameter:
294: snes - SNES context
295: X - KKT Vector
296: *ctx - pdipm context
298: Output Parameter:
299: J - Hessian matrix
300: Jpre - Preconditioner
301: */
302: PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
303: {
304: PetscErrorCode ierr;
305: Tao tao=(Tao)ctx;
306: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
307: PetscInt i,row,cols[2],Jrstart,rjstart,nc,j;
308: const PetscInt *aj,*ranges,*Jranges,*rranges,*cranges;
309: const PetscScalar *Xarr,*aa;
310: PetscScalar vals[2];
311: PetscInt proc,nx_all,*nce_all=pdipm->nce_all;
312: MPI_Comm comm;
313: PetscMPIInt rank,size;
314: Mat jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;
317: PetscObjectGetComm((PetscObject)snes,&comm);
318: MPI_Comm_rank(comm,&rank);
319: MPI_Comm_rank(comm,&size);
321: MatGetOwnershipRanges(Jpre,&Jranges);
322: MatGetOwnershipRange(Jpre,&Jrstart,NULL);
323: MatGetOwnershipRangesColumn(tao->hessian,&rranges);
324: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
326: VecGetArrayRead(X,&Xarr);
328: if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
329: vals[0] = 1.0;
330: for (i=0; i < pdipm->nci; i++) {
331: row = Jrstart + pdipm->off_z + i;
332: cols[0] = Jrstart + pdipm->off_lambdai + i;
333: cols[1] = row;
334: vals[1] = Xarr[pdipm->off_lambdai + i]/Xarr[pdipm->off_z + i];
335: MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
336: }
337: } else {
338: /* (2) insert Z and Ci to Jpre -- overwrite existing values */
339: for (i=0; i < pdipm->nci; i++) {
340: row = Jrstart + pdipm->off_z + i;
341: cols[0] = Jrstart + pdipm->off_lambdai + i;
342: cols[1] = row;
343: vals[0] = Xarr[pdipm->off_z + i];
344: vals[1] = Xarr[pdipm->off_lambdai + i];
345: MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
346: }
347: }
349: /* (3) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
350: if (pdipm->Ng) {
351: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
352: for (i=0; i<pdipm->ng; i++){
353: row = Jrstart + pdipm->off_lambdae + i;
355: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
356: proc = 0;
357: for (j=0; j < nc; j++) {
358: while (aj[j] >= cranges[proc+1]) proc++;
359: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
360: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
361: }
362: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
363: if (pdipm->kkt_pd) {
364: /* (3) insert 2nd row block of Jpre: [ grad g, \delta_c*I, 0, 0] */
365: MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
366: }
367: }
368: }
370: if (pdipm->Nh) {
371: /* (4) insert 3nd row block of Jpre: [ -grad h, 0, deltac, I] */
372: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
373: for (i=0; i < pdipm->nh; i++){
374: row = Jrstart + pdipm->off_lambdai + i;
375: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
376: proc = 0;
377: for (j=0; j < nc; j++) {
378: while (aj[j] >= cranges[proc+1]) proc++;
379: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
380: MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
381: }
382: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
383: if (pdipm->kkt_pd) {
384: MatSetValue(Jpre,row,row,-pdipm->deltac,INSERT_VALUES);
385: }
386: }
387: }
389: /* (5) insert Wxx, grad g' and -grad h' to Jpre */
390: if (pdipm->Ng) {
391: MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
392: }
393: if (pdipm->Nh) {
394: MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
395: }
397: VecPlaceArray(pdipm->x,Xarr);
398: TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
399: VecResetArray(pdipm->x);
401: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
402: for (i=0; i<pdipm->nx; i++) {
403: row = Jrstart + i;
405: /* insert Wxx */
406: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
407: proc = 0;
408: for (j=0; j < nc; j++) {
409: while (aj[j] >= cranges[proc+1]) proc++;
410: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
411: if (row == cols[0] && pdipm->kkt_pd) {
412: /* add diagonal shift to Wxx component */
413: MatSetValue(Jpre,row,cols[0],aa[j]+pdipm->deltaw,INSERT_VALUES);
414: } else {
415: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
416: }
417: }
418: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
420: if (pdipm->ng) {
421: /* insert grad g' */
422: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
423: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
424: proc = 0;
425: for (j=0; j < nc; j++) {
426: /* find row ownership of */
427: while (aj[j] >= ranges[proc+1]) proc++;
428: nx_all = rranges[proc+1] - rranges[proc];
429: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
430: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
431: }
432: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
433: }
435: if (pdipm->nh) {
436: /* insert -grad h' */
437: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
438: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
439: proc = 0;
440: for (j=0; j < nc; j++) {
441: /* find row ownership of */
442: while (aj[j] >= ranges[proc+1]) proc++;
443: nx_all = rranges[proc+1] - rranges[proc];
444: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
445: MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
446: }
447: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
448: }
449: }
450: VecRestoreArrayRead(X,&Xarr);
452: /* (6) assemble Jpre and J */
453: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
454: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
456: if (J != Jpre) {
457: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
458: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
459: }
460: return(0);
461: }
463: /*
464: TaoSnesFunction_PDIPM - Evaluate KKT function at X
466: Input Parameter:
467: snes - SNES context
468: X - KKT Vector
469: *ctx - pdipm
471: Output Parameter:
472: F - Updated Lagrangian vector
473: */
474: PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
475: {
476: PetscErrorCode ierr;
477: Tao tao=(Tao)ctx;
478: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
479: PetscScalar *Farr,*tmparr;
480: Vec x,L1;
481: PetscInt i;
482: PetscReal res[2],cnorm[2];
483: const PetscScalar *Xarr,*carr,*zarr,*larr;
486: VecSet(F,0.0);
488: VecGetArrayRead(X,&Xarr);
489: VecGetArray(F,&Farr);
491: /* (0) Evaluate f, fx, Gx, Hx at X.x Note: pdipm->x is not changed below */
492: x = pdipm->x;
493: VecPlaceArray(x,Xarr);
494: TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);
496: /* Update ce, ci, and Jci at X.x */
497: TaoPDIPMUpdateConstraints(tao,x);
498: VecResetArray(x);
500: /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
501: L1 = pdipm->x;
502: VecPlaceArray(L1,Farr);
503: if (pdipm->Nci) {
504: if (pdipm->Nh) {
505: /* L1 += gradH'*DI. Note: tao->DI is not changed below */
506: VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
507: MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
508: VecResetArray(tao->DI);
509: }
511: /* L1 += Jci_xb'*lambdai_xb */
512: VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
513: MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
514: VecResetArray(pdipm->lambdai_xb);
516: /* (1.4) L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
517: VecScale(L1,-1.0);
518: }
520: /* L1 += fx */
521: VecAXPY(L1,1.0,tao->gradient);
523: if (pdipm->Nce) {
524: if (pdipm->Ng) {
525: /* L1 += gradG'*DE. Note: tao->DE is not changed below */
526: VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
527: MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
528: VecResetArray(tao->DE);
529: }
530: if (pdipm->Nxfixed) {
531: /* L1 += Jce_xfixed'*lambdae_xfixed */
532: VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
533: MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
534: VecResetArray(pdipm->lambdae_xfixed);
535: }
536: }
537: VecNorm(L1,NORM_2,&res[0]);
538: VecResetArray(L1);
540: /* (2) L2 = ce(x) */
541: if (pdipm->Nce) {
542: VecGetArrayRead(pdipm->ce,&carr);
543: for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
544: VecRestoreArrayRead(pdipm->ce,&carr);
545: }
546: VecNorm(pdipm->ce,NORM_2,&cnorm[0]);
548: if (pdipm->Nci) {
549: if (pdipm->solve_symmetric_kkt) {
550: /* (3) L3 = ci(x) - z;
551: (4) L4 = Lambdai * e - mu/z *e
552: */
553: VecGetArrayRead(pdipm->ci,&carr);
554: larr = Xarr+pdipm->off_lambdai;
555: zarr = Xarr+pdipm->off_z;
556: for (i=0; i<pdipm->nci; i++) {
557: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
558: Farr[pdipm->off_z + i] = larr[i] - pdipm->mu/zarr[i];
559: }
560: VecRestoreArrayRead(pdipm->ci,&carr);
561: } else {
562: /* (3) L3 = ci(x) - z;
563: (4) L4 = Z * Lambdai * e - mu * e
564: */
565: VecGetArrayRead(pdipm->ci,&carr);
566: larr = Xarr+pdipm->off_lambdai;
567: zarr = Xarr+pdipm->off_z;
568: for (i=0; i<pdipm->nci; i++) {
569: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
570: Farr[pdipm->off_z + i] = zarr[i]*larr[i] - pdipm->mu;
571: }
572: VecRestoreArrayRead(pdipm->ci,&carr);
573: }
574: }
576: VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
577: VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
578: VecResetArray(pdipm->ci);
580: /* note: pdipm->z is not changed below */
581: if (pdipm->z) {
582: if (pdipm->solve_symmetric_kkt) {
583: VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
584: if (pdipm->Nci) {
585: zarr = Xarr+pdipm->off_z;
586: VecGetArrayWrite(pdipm->z,&tmparr);
587: for (i=0; i<pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
588: VecRestoreArrayWrite(pdipm->z,&tmparr);
589: }
591: VecNorm(pdipm->z,NORM_2,&res[1]);
592: if (pdipm->Nci) {
593: zarr = Xarr+pdipm->off_z;
594: VecGetArray(pdipm->z,&tmparr);
595: for (i=0; i<pdipm->nci; i++) {
596: tmparr[i] /= Xarr[pdipm->off_z + i];
597: }
598: VecRestoreArray(pdipm->z,&tmparr);
599: }
600: VecResetArray(pdipm->z);
601: } else {
602: VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
603: VecNorm(pdipm->z,NORM_2,&res[1]);
604: VecResetArray(pdipm->z);
605: }
606: } else res[1] = 0.0;
608: tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);
609: tao->cnorm = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);
611: VecRestoreArrayRead(X,&Xarr);
612: VecRestoreArray(F,&Farr);
613: return(0);
614: }
616: #include <petsc/private/matimpl.h>
617: /*
618: PDIPMLineSearch - Custom line search used with PDIPM.
620: Collective on TAO
622: Notes:
623: PDIPMLineSearch employs a simple backtracking line-search to keep
624: the slack variables (z) and inequality constraints lagrange multipliers
625: (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
626: alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
627: slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
628: are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
629: is also updated as mu = mu + z'lambdai/Nci
630: */
631: PetscErrorCode PDIPMLineSearch(SNESLineSearch linesearch,void *ctx)
632: {
633: PetscErrorCode ierr;
634: Tao tao=(Tao)ctx;
635: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
636: SNES snes;
637: KSP ksp;
638: PC pc;
639: PCType ptype;
640: Mat Factor;
641: Vec X,F,Y,W,G;
642: PetscInt i,iter,nneg,nzero,npos;
643: PetscReal alpha_p=1.0,alpha_d=1.0,alpha[4];
644: PetscScalar *Xarr,*z,*lambdai,dot,*taosolarr;
645: const PetscScalar *dXarr,*dz,*dlambdai;
646: PetscBool isCHOL;
649: SNESLineSearchGetSNES(linesearch,&snes);
650: SNESGetIterationNumber(snes,&iter);
652: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
653: SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);
655: VecGetArray(X,&Xarr);
656: VecGetArrayRead(Y,&dXarr);
657: z = Xarr + pdipm->off_z;
658: dz = dXarr + pdipm->off_z;
659: for (i=0; i < pdipm->nci; i++) {
660: if (z[i] - dz[i] < 0.0) {
661: alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
662: }
663: }
665: lambdai = Xarr + pdipm->off_lambdai;
666: dlambdai = dXarr + pdipm->off_lambdai;
668: for (i=0; i<pdipm->nci; i++) {
669: if (lambdai[i] - dlambdai[i] < 0.0) {
670: alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
671: }
672: }
674: alpha[0] = alpha_p;
675: alpha[1] = alpha_d;
676: VecRestoreArrayRead(Y,&dXarr);
677: VecRestoreArray(X,&Xarr);
679: /* alpha = min(alpha) over all processes */
680: MPI_Allreduce(alpha,alpha+2,2,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)tao));
682: alpha_p = alpha[2];
683: alpha_d = alpha[3];
685: VecGetArray(X,&Xarr);
686: VecGetArrayRead(Y,&dXarr);
687: for (i=0; i<pdipm->nx; i++) {
688: Xarr[i] = Xarr[i] - alpha_p * dXarr[i];
689: }
691: for (i=0; i<pdipm->nce; i++) {
692: Xarr[i+pdipm->off_lambdae] = Xarr[i+pdipm->off_lambdae] - alpha_d * dXarr[i+pdipm->off_lambdae];
693: }
695: for (i=0; i<pdipm->nci; i++) {
696: Xarr[i+pdipm->off_lambdai] = Xarr[i+pdipm->off_lambdai] - alpha_d * dXarr[i+pdipm->off_lambdai];
697: }
699: for (i=0; i<pdipm->nci; i++) {
700: Xarr[i+pdipm->off_z] = Xarr[i+pdipm->off_z] - alpha_p * dXarr[i+pdipm->off_z];
701: }
703: VecGetArray(tao->solution,&taosolarr);
704: PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
705: VecRestoreArray(tao->solution,&taosolarr);
708: VecRestoreArray(X,&Xarr);
709: VecRestoreArrayRead(Y,&dXarr);
711: /* Evaluate F at X */
712: SNESComputeFunction(snes,X,F);
713: SNESLineSearchComputeNorms(linesearch); /* must call this func, do not know why */
715: /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
716: if (pdipm->z) {
717: VecDot(pdipm->z,pdipm->lambdai,&dot);
718: } else dot = 0.0;
720: /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu) */
721: pdipm->mu = pdipm->mu_update_factor * dot/pdipm->Nci;
723: /* Update F; get tao->residual and tao->cnorm */
724: TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);
726: tao->niter++;
727: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
728: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
730: (*tao->ops->convergencetest)(tao,tao->cnvP);
731: if (tao->reason) {
732: SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
733: }
735: if (!pdipm->kkt_pd) return(0);
737: /* Get the inertia of Cholesky factor to set shifts for next SNES interation */
738: SNESGetKSP(snes,&ksp);
739: KSPGetPC(ksp,&pc);
740: PCGetType(pc,&ptype);
741: PetscObjectTypeCompare((PetscObject)pc,PCCHOLESKY,&isCHOL);
743: if (isCHOL) {
744: PetscMPIInt size;
745: PCFactorGetMatrix(pc,&Factor);
746: MPI_Comm_size(PetscObjectComm((PetscObject)Factor),&size);
747: if (Factor->ops->getinertia) {
748: #if defined(PETSC_HAVE_MUMPS)
749: MatSolverType stype;
750: PetscBool isMUMPS;
751: PCFactorGetMatSolverType(pc,&stype);
752: PetscStrcmp(stype, MATSOLVERMUMPS, &isMUMPS);
753: if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
754: MatMumpsSetIcntl(Factor,24,1);
755: if (size > 1) {
756: MatMumpsSetIcntl(Factor,13,1);
757: }
758: }
759: #else
760: if (size > 1) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires external package MUMPS");
761: #endif
762: MatGetInertia(Factor,&nneg,&nzero,&npos);
764: if (npos < pdipm->Nx+pdipm->Nci) {
765: pdipm->deltaw = PetscMax(pdipm->lastdeltaw/3, 1.e-4*PETSC_MACHINE_EPSILON);
766: PetscInfo5(tao,"Test reduced deltaw=%g; previous MatInertia: nneg %d, nzero %d, npos %d(<%d)\n",pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
767: TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
768: PCSetUp(pc);
769: MatGetInertia(Factor,&nneg,&nzero,&npos);
771: if (npos < pdipm->Nx+pdipm->Nci) {
772: pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
773: while (npos < pdipm->Nx+pdipm->Nci && pdipm->deltaw <= 1.e10) { /* increase deltaw */
774: PetscInfo5(tao," deltaw=%g fails, MatInertia: nneg %d, nzero %d, npos %d(<%d)\n",pdipm->deltaw,nneg,nzero,npos,pdipm->Nx+pdipm->Nci);
775: pdipm->deltaw = PetscMin(8*pdipm->deltaw,PetscPowReal(10,20));
776: TaoSNESJacobian_PDIPM(snes,X, pdipm->K, pdipm->K, tao);
777: PCSetUp(pc);
778: MatGetInertia(Factor,&nneg,&nzero,&npos);
779: }
781: if (pdipm->deltaw >= 1.e10) {
782: SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_CONV_FAILED,"Reached maximum delta w will not converge, try different inital x0");
783: }
784: PetscInfo1(tao,"Updated deltaw %g\n",pdipm->deltaw);
785: pdipm->lastdeltaw = pdipm->deltaw;
786: pdipm->deltaw = 0.0;
787: }
788: }
790: if (nzero) { /* Jacobian is singular */
791: if (pdipm->deltac == 0.0) {
792: pdipm->deltac = 1.e8*PETSC_MACHINE_EPSILON;
793: } else {
794: pdipm->deltac = pdipm->deltac*PetscPowReal(pdipm->mu,.25);
795: }
796: PetscInfo4(tao,"Updated deltac=%g, MatInertia: nneg %D, nzero %D(!=0), npos %D\n",pdipm->deltac,nneg,nzero,npos);
797: }
798: } else
799: SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Requires an external package that supports MatGetInertia()");
800: }
801: return(0);
802: }
804: /*
805: TaoSolve_PDIPM
807: Input Parameter:
808: tao - TAO context
810: Output Parameter:
811: tao - TAO context
812: */
813: PetscErrorCode TaoSolve_PDIPM(Tao tao)
814: {
815: PetscErrorCode ierr;
816: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
817: SNESLineSearch linesearch; /* SNESLineSearch context */
818: Vec dummy;
821: if (!tao->constraints_equality && !tao->constraints_inequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_NULL,"Equality and inequality constraints are not set. Either set them or switch to a different algorithm");
823: /* Initialize all variables */
824: TaoPDIPMInitializeSolution(tao);
826: /* Set linesearch */
827: SNESGetLineSearch(pdipm->snes,&linesearch);
828: SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
829: SNESLineSearchShellSetUserFunc(linesearch,PDIPMLineSearch,tao);
830: SNESLineSearchSetFromOptions(linesearch);
832: tao->reason = TAO_CONTINUE_ITERATING;
834: /* -tao_monitor for iteration 0 and check convergence */
835: VecDuplicate(pdipm->X,&dummy);
836: TaoSNESFunction_PDIPM(pdipm->snes,pdipm->X,dummy,(void*)tao);
838: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
839: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
840: VecDestroy(&dummy);
841: (*tao->ops->convergencetest)(tao,tao->cnvP);
842: if (tao->reason) {
843: SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
844: }
846: while (tao->reason == TAO_CONTINUE_ITERATING) {
847: SNESConvergedReason reason;
848: SNESSolve(pdipm->snes,NULL,pdipm->X);
850: /* Check SNES convergence */
851: SNESGetConvergedReason(pdipm->snes,&reason);
852: if (reason < 0) {
853: PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes),"SNES solve did not converged due to reason %D\n",reason);
854: }
856: /* Check TAO convergence */
857: if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
858: }
859: return(0);
860: }
862: /*
863: TaoView_PDIPM - View PDIPM
865: Input Parameter:
866: tao - TAO object
867: viewer - PetscViewer
869: Output:
870: */
871: PetscErrorCode TaoView_PDIPM(Tao tao,PetscViewer viewer)
872: {
873: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
877: tao->constrained = PETSC_TRUE;
878: PetscViewerASCIIPushTab(viewer);
879: PetscViewerASCIIPrintf(viewer,"Number of prime=%D, Number of dual=%D\n",pdipm->Nx+pdipm->Nci,pdipm->Nce + pdipm->Nci);
880: if (pdipm->kkt_pd) {
881: PetscViewerASCIIPrintf(viewer,"KKT shifts deltaw=%g, deltac=%g\n",pdipm->deltaw,pdipm->deltac);
882: }
883: PetscViewerASCIIPopTab(viewer);
884: return(0);
885: }
887: /*
888: TaoSetup_PDIPM - Sets up tao and pdipm
890: Input Parameter:
891: tao - TAO object
893: Output: pdipm - initialized object
894: */
895: PetscErrorCode TaoSetup_PDIPM(Tao tao)
896: {
897: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
898: PetscErrorCode ierr;
899: MPI_Comm comm;
900: PetscMPIInt rank,size;
901: PetscInt row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
902: PetscInt offset,*xa,*xb,i,j,rstart,rend;
903: PetscScalar one=1.0,neg_one=-1.0,*Xarr;
904: const PetscInt *cols,*rranges,*cranges,*aj,*ranges;
905: const PetscScalar *aa;
906: Mat J,jac_equality_trans,jac_inequality_trans;
907: Mat Jce_xfixed_trans,Jci_xb_trans;
908: PetscInt *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];
911: PetscObjectGetComm((PetscObject)tao,&comm);
912: MPI_Comm_rank(comm,&rank);
913: MPI_Comm_size(comm,&size);
915: /* (1) Setup Bounds and create Tao vectors */
916: TaoPDIPMSetUpBounds(tao);
918: if (!tao->gradient) {
919: VecDuplicate(tao->solution,&tao->gradient);
920: VecDuplicate(tao->solution,&tao->stepdirection);
921: }
923: /* (2) Get sizes */
924: /* Size of vector x - This is set by TaoSetInitialVector */
925: VecGetSize(tao->solution,&pdipm->Nx);
926: VecGetLocalSize(tao->solution,&pdipm->nx);
928: /* Size of equality constraints and vectors */
929: if (tao->constraints_equality) {
930: VecGetSize(tao->constraints_equality,&pdipm->Ng);
931: VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
932: } else {
933: pdipm->ng = pdipm->Ng = 0;
934: }
936: pdipm->nce = pdipm->ng + pdipm->nxfixed;
937: pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;
939: /* Size of inequality constraints and vectors */
940: if (tao->constraints_inequality) {
941: VecGetSize(tao->constraints_inequality,&pdipm->Nh);
942: VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
943: } else {
944: pdipm->nh = pdipm->Nh = 0;
945: }
947: pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2*pdipm->nxbox;
948: pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2*pdipm->Nxbox;
950: /* Full size of the KKT system to be solved */
951: pdipm->n = pdipm->nx + pdipm->nce + 2*pdipm->nci;
952: pdipm->N = pdipm->Nx + pdipm->Nce + 2*pdipm->Nci;
954: /* (3) Offsets for subvectors */
955: pdipm->off_lambdae = pdipm->nx;
956: pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
957: pdipm->off_z = pdipm->off_lambdai + pdipm->nci;
959: /* (4) Create vectors and subvectors */
960: /* Ce and Ci vectors */
961: VecCreate(comm,&pdipm->ce);
962: VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
963: VecSetFromOptions(pdipm->ce);
965: VecCreate(comm,&pdipm->ci);
966: VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
967: VecSetFromOptions(pdipm->ci);
969: /* X=[x; lambdae; lambdai; z] for the big KKT system */
970: VecCreate(comm,&pdipm->X);
971: VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
972: VecSetFromOptions(pdipm->X);
974: /* Subvectors; they share local arrays with X */
975: VecGetArray(pdipm->X,&Xarr);
976: /* x shares local array with X.x */
977: if (pdipm->Nx) {
978: VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
979: }
981: /* lambdae shares local array with X.lambdae */
982: if (pdipm->Nce) {
983: VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
984: }
986: /* tao->DE shares local array with X.lambdae_g */
987: if (pdipm->Ng) {
988: VecCreateMPIWithArray(comm,1,pdipm->ng,pdipm->Ng,Xarr+pdipm->off_lambdae,&tao->DE);
990: VecCreate(comm,&pdipm->lambdae_xfixed);
991: VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
992: VecSetFromOptions(pdipm->lambdae_xfixed);
993: }
995: if (pdipm->Nci) {
996: /* lambdai shares local array with X.lambdai */
997: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_lambdai,&pdipm->lambdai);
999: /* z for slack variables; it shares local array with X.z */
1000: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_z,&pdipm->z);
1001: }
1003: /* tao->DI which shares local array with X.lambdai_h */
1004: if (pdipm->Nh) {
1005: VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
1006: }
1008: VecCreate(comm,&pdipm->lambdai_xb);
1009: VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
1010: VecSetFromOptions(pdipm->lambdai_xb);
1012: VecRestoreArray(pdipm->X,&Xarr);
1014: /* (5) Create Jacobians Jce_xfixed and Jci */
1015: /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
1016: if (pdipm->Nxfixed) {
1017: /* Create Jce_xfixed */
1018: MatCreate(comm,&pdipm->Jce_xfixed);
1019: MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1020: MatSetFromOptions(pdipm->Jce_xfixed);
1021: MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
1022: MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);
1024: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
1025: ISGetIndices(pdipm->isxfixed,&cols);
1026: k = 0;
1027: for (row = Jcrstart; row < Jcrend; row++) {
1028: MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
1029: k++;
1030: }
1031: ISRestoreIndices(pdipm->isxfixed, &cols);
1032: MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1033: MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
1034: }
1036: /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
1037: MatCreate(comm,&pdipm->Jci_xb);
1038: MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
1039: MatSetFromOptions(pdipm->Jci_xb);
1040: MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
1041: MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);
1043: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
1044: offset = Jcrstart;
1045: if (pdipm->Nxub) {
1046: /* Add xub to Jci_xb */
1047: ISGetIndices(pdipm->isxub,&cols);
1048: k = 0;
1049: for (row = offset; row < offset + pdipm->nxub; row++) {
1050: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1051: k++;
1052: }
1053: ISRestoreIndices(pdipm->isxub, &cols);
1054: }
1056: if (pdipm->Nxlb) {
1057: /* Add xlb to Jci_xb */
1058: ISGetIndices(pdipm->isxlb,&cols);
1059: k = 0;
1060: offset += pdipm->nxub;
1061: for (row = offset; row < offset + pdipm->nxlb; row++) {
1062: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
1063: k++;
1064: }
1065: ISRestoreIndices(pdipm->isxlb, &cols);
1066: }
1068: /* Add xbox to Jci_xb */
1069: if (pdipm->Nxbox) {
1070: ISGetIndices(pdipm->isxbox,&cols);
1071: k = 0;
1072: offset += pdipm->nxlb;
1073: for (row = offset; row < offset + pdipm->nxbox; row++) {
1074: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
1075: tmp = row + pdipm->nxbox;
1076: MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
1077: k++;
1078: }
1079: ISRestoreIndices(pdipm->isxbox, &cols);
1080: }
1082: MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1083: MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
1084: /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */
1086: /* (6) Set up ISs for PC Fieldsplit */
1087: if (pdipm->solve_reduced_kkt) {
1088: PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
1089: for (i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1090: for (i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;
1092: ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
1093: ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
1094: }
1096: /* (7) Gather offsets from all processes */
1097: PetscMalloc1(size,&pdipm->nce_all);
1099: /* Get rstart of KKT matrix */
1100: MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
1101: rstart -= pdipm->n;
1103: MPI_Allgather(&pdipm->nce,1,MPIU_INT,pdipm->nce_all,1,MPIU_INT,comm);
1105: PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
1106: MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
1107: MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
1108: MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);
1110: MatGetOwnershipRanges(tao->hessian,&rranges);
1111: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
1113: if (pdipm->Ng) {
1114: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
1115: MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
1116: }
1117: if (pdipm->Nh) {
1118: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
1119: MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
1120: }
1122: /* Count dnz,onz for preallocation of KKT matrix */
1123: jac_equality_trans = pdipm->jac_equality_trans;
1124: jac_inequality_trans = pdipm->jac_inequality_trans;
1125: nce_all = pdipm->nce_all;
1127: if (pdipm->Nxfixed) {
1128: MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
1129: }
1130: MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);
1132: MatPreallocateInitialize(comm,pdipm->n,pdipm->n,dnz,onz);
1134: /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
1135: TaoPDIPMEvaluateFunctionsAndJacobians(tao,pdipm->x);
1136: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
1138: /* Insert tao->hessian */
1139: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
1140: for (i=0; i<pdipm->nx; i++){
1141: row = rstart + i;
1143: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1144: proc = 0;
1145: for (j=0; j < nc; j++) {
1146: while (aj[j] >= cranges[proc+1]) proc++;
1147: col = aj[j] - cranges[proc] + Jranges[proc];
1148: MatPreallocateSet(row,1,&col,dnz,onz);
1149: }
1150: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
1152: if (pdipm->ng) {
1153: /* Insert grad g' */
1154: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1155: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
1156: proc = 0;
1157: for (j=0; j < nc; j++) {
1158: /* find row ownership of */
1159: while (aj[j] >= ranges[proc+1]) proc++;
1160: nx_all = rranges[proc+1] - rranges[proc];
1161: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1162: MatPreallocateSet(row,1,&col,dnz,onz);
1163: }
1164: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1165: }
1167: /* Insert Jce_xfixed^T' */
1168: if (pdipm->nxfixed) {
1169: MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1170: MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1171: proc = 0;
1172: for (j=0; j < nc; j++) {
1173: /* find row ownership of */
1174: while (aj[j] >= ranges[proc+1]) proc++;
1175: nx_all = rranges[proc+1] - rranges[proc];
1176: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1177: MatPreallocateSet(row,1,&col,dnz,onz);
1178: }
1179: MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1180: }
1182: if (pdipm->nh) {
1183: /* Insert -grad h' */
1184: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1185: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1186: proc = 0;
1187: for (j=0; j < nc; j++) {
1188: /* find row ownership of */
1189: while (aj[j] >= ranges[proc+1]) proc++;
1190: nx_all = rranges[proc+1] - rranges[proc];
1191: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1192: MatPreallocateSet(row,1,&col,dnz,onz);
1193: }
1194: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1195: }
1197: /* Insert Jci_xb^T' */
1198: MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1199: MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1200: proc = 0;
1201: for (j=0; j < nc; j++) {
1202: /* find row ownership of */
1203: while (aj[j] >= ranges[proc+1]) proc++;
1204: nx_all = rranges[proc+1] - rranges[proc];
1205: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1206: MatPreallocateSet(row,1,&col,dnz,onz);
1207: }
1208: MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1209: }
1211: /* 2nd Row block of KKT matrix: [grad Ce, deltac*I, 0, 0] */
1212: if (pdipm->Ng) {
1213: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
1214: for (i=0; i < pdipm->ng; i++){
1215: row = rstart + pdipm->off_lambdae + i;
1217: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1218: proc = 0;
1219: for (j=0; j < nc; j++) {
1220: while (aj[j] >= cranges[proc+1]) proc++;
1221: col = aj[j] - cranges[proc] + Jranges[proc];
1222: MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1223: }
1224: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1225: }
1226: }
1227: /* Jce_xfixed */
1228: if (pdipm->Nxfixed) {
1229: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1230: for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1231: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1233: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1234: if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1236: proc = 0;
1237: j = 0;
1238: while (cols[j] >= cranges[proc+1]) proc++;
1239: col = cols[j] - cranges[proc] + Jranges[proc];
1240: MatPreallocateSet(row,1,&col,dnz,onz);
1241: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1242: }
1243: }
1245: /* 3rd Row block of KKT matrix: [ gradCi, 0, deltac*I, -I] */
1246: if (pdipm->Nh) {
1247: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
1248: for (i=0; i < pdipm->nh; i++){
1249: row = rstart + pdipm->off_lambdai + i;
1251: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1252: proc = 0;
1253: for (j=0; j < nc; j++) {
1254: while (aj[j] >= cranges[proc+1]) proc++;
1255: col = aj[j] - cranges[proc] + Jranges[proc];
1256: MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1257: }
1258: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1259: }
1260: /* I */
1261: for (i=0; i < pdipm->nh; i++){
1262: row = rstart + pdipm->off_lambdai + i;
1263: col = rstart + pdipm->off_z + i;
1264: MatPreallocateSet(row,1,&col,dnz,onz);
1265: }
1266: }
1268: /* Jci_xb */
1269: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1270: for (i=0; i < (pdipm->nci - pdipm->nh); i++){
1271: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1273: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1274: if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1275: proc = 0;
1276: for (j=0; j < nc; j++) {
1277: while (cols[j] >= cranges[proc+1]) proc++;
1278: col = cols[j] - cranges[proc] + Jranges[proc];
1279: MatPreallocateSet(row,1,&col,dnz,onz);
1280: }
1281: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1282: /* I */
1283: col = rstart + pdipm->off_z + pdipm->nh + i;
1284: MatPreallocateSet(row,1,&col,dnz,onz);
1285: }
1287: /* 4-th Row block of KKT matrix: Z and Ci */
1288: for (i=0; i < pdipm->nci; i++) {
1289: row = rstart + pdipm->off_z + i;
1290: cols1[0] = rstart + pdipm->off_lambdai + i;
1291: cols1[1] = row;
1292: MatPreallocateSet(row,2,cols1,dnz,onz);
1293: }
1295: /* diagonal entry */
1296: for (i=0; i<pdipm->n; i++) dnz[i]++; /* diagonal entry */
1298: /* Create KKT matrix */
1299: MatCreate(comm,&J);
1300: MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1301: MatSetFromOptions(J);
1302: MatSeqAIJSetPreallocation(J,0,dnz);
1303: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1304: /* MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); */
1305: MatPreallocateFinalize(dnz,onz);
1306: pdipm->K = J;
1308: /* (8) Set up nonlinear solver SNES */
1309: SNESSetFunction(pdipm->snes,NULL,TaoSNESFunction_PDIPM,(void*)tao);
1310: SNESSetJacobian(pdipm->snes,J,J,TaoSNESJacobian_PDIPM,(void*)tao);
1312: if (pdipm->solve_reduced_kkt) {
1313: PC pc;
1314: KSPGetPC(tao->ksp,&pc);
1315: PCSetType(pc,PCFIELDSPLIT);
1316: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1317: PCFieldSplitSetIS(pc,"2",pdipm->is2);
1318: PCFieldSplitSetIS(pc,"1",pdipm->is1);
1319: }
1320: SNESSetFromOptions(pdipm->snes);
1322: /* (9) Insert constant entries to K */
1323: /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1324: MatGetOwnershipRange(J,&rstart,&rend);
1325: for (i=rstart; i<rend; i++){
1326: MatSetValue(J,i,i,0.0,INSERT_VALUES);
1327: }
1328: /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1329: if (pdipm->kkt_pd){
1330: for (i=0; i<pdipm->nh; i++){
1331: row = rstart + i;
1332: MatSetValue(J,row,row,pdipm->deltaw,INSERT_VALUES);
1333: }
1334: }
1336: /* Row block of K: [ grad Ce, 0, 0, 0] */
1337: if (pdipm->Nxfixed) {
1338: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1339: for (i=0; i < (pdipm->nce - pdipm->ng); i++){
1340: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1342: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1343: proc = 0;
1344: for (j=0; j < nc; j++) {
1345: while (cols[j] >= cranges[proc+1]) proc++;
1346: col = cols[j] - cranges[proc] + Jranges[proc];
1347: MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1348: MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1349: }
1350: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1351: }
1352: }
1354: /* Row block of K: [ -grad Ci, 0, 0, I] */
1355: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1356: for (i=0; i < pdipm->nci - pdipm->nh; i++){
1357: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1359: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1360: proc = 0;
1361: for (j=0; j < nc; j++) {
1362: while (cols[j] >= cranges[proc+1]) proc++;
1363: col = cols[j] - cranges[proc] + Jranges[proc];
1364: MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1365: MatSetValue(J,row,col,-aa[j],INSERT_VALUES);
1366: }
1367: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1369: col = rstart + pdipm->off_z + pdipm->nh + i;
1370: MatSetValue(J,row,col,1,INSERT_VALUES);
1371: }
1373: for (i=0; i < pdipm->nh; i++){
1374: row = rstart + pdipm->off_lambdai + i;
1375: col = rstart + pdipm->off_z + i;
1376: MatSetValue(J,row,col,1,INSERT_VALUES);
1377: }
1379: /* Row block of K: [ 0, 0, I, ...] */
1380: for (i=0; i < pdipm->nci; i++){
1381: row = rstart + pdipm->off_z + i;
1382: col = rstart + pdipm->off_lambdai + i;
1383: MatSetValue(J,row,col,1,INSERT_VALUES);
1384: }
1386: if (pdipm->Nxfixed) {
1387: MatDestroy(&Jce_xfixed_trans);
1388: }
1389: MatDestroy(&Jci_xb_trans);
1390: PetscFree3(ng_all,nh_all,Jranges);
1391: return(0);
1392: }
1394: /*
1395: TaoDestroy_PDIPM - Destroys the pdipm object
1397: Input:
1398: full pdipm
1400: Output:
1401: Destroyed pdipm
1402: */
1403: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1404: {
1405: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1409: /* Freeing Vectors assocaiated with KKT (X) */
1410: VecDestroy(&pdipm->x); /* Solution x */
1411: VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1412: VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1413: VecDestroy(&pdipm->z); /* Slack variables */
1414: VecDestroy(&pdipm->X); /* Big KKT system vector [x; lambdae; lambdai; z] */
1416: /* work vectors */
1417: VecDestroy(&pdipm->lambdae_xfixed);
1418: VecDestroy(&pdipm->lambdai_xb);
1420: /* Legrangian equality and inequality Vec */
1421: VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1422: VecDestroy(&pdipm->ci); /* Vec of inequality constraints */
1424: /* Matrices */
1425: MatDestroy(&pdipm->Jce_xfixed);
1426: MatDestroy(&pdipm->Jci_xb); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1427: MatDestroy(&pdipm->K);
1429: /* Index Sets */
1430: if (pdipm->Nxub) {
1431: ISDestroy(&pdipm->isxub); /* Finite upper bound only -inf < x < ub */
1432: }
1434: if (pdipm->Nxlb) {
1435: ISDestroy(&pdipm->isxlb); /* Finite lower bound only lb <= x < inf */
1436: }
1438: if (pdipm->Nxfixed) {
1439: ISDestroy(&pdipm->isxfixed); /* Fixed variables lb = x = ub */
1440: }
1442: if (pdipm->Nxbox) {
1443: ISDestroy(&pdipm->isxbox); /* Boxed variables lb <= x <= ub */
1444: }
1446: if (pdipm->Nxfree) {
1447: ISDestroy(&pdipm->isxfree); /* Free variables -inf <= x <= inf */
1448: }
1450: if (pdipm->solve_reduced_kkt) {
1451: ISDestroy(&pdipm->is1);
1452: ISDestroy(&pdipm->is2);
1453: }
1455: /* SNES */
1456: SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1457: PetscFree(pdipm->nce_all);
1458: MatDestroy(&pdipm->jac_equality_trans);
1459: MatDestroy(&pdipm->jac_inequality_trans);
1461: /* Destroy pdipm */
1462: PetscFree(tao->data); /* Holding locations of pdipm */
1464: /* Destroy Dual */
1465: VecDestroy(&tao->DE); /* equality dual */
1466: VecDestroy(&tao->DI); /* dinequality dual */
1467: return(0);
1468: }
1470: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1471: {
1472: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1476: PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1477: PetscOptionsReal("-tao_pdipm_push_init_slack","parameter to push initial slack variables away from bounds",NULL,pdipm->push_init_slack,&pdipm->push_init_slack,NULL);
1478: PetscOptionsReal("-tao_pdipm_push_init_lambdai","parameter to push initial (inequality) dual variables away from bounds",NULL,pdipm->push_init_lambdai,&pdipm->push_init_lambdai,NULL);
1479: PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1480: PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1481: PetscOptionsBool("-tao_pdipm_symmetric_kkt","Solve non reduced symmetric KKT system",NULL,pdipm->solve_symmetric_kkt,&pdipm->solve_symmetric_kkt,NULL);
1482: PetscOptionsBool("-tao_pdipm_kkt_shift_pd","Add shifts to make KKT matrix positive definite",NULL,pdipm->kkt_pd,&pdipm->kkt_pd,NULL);
1483: PetscOptionsTail();
1484: return(0);
1485: }
1487: /*MC
1488: TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.
1490: Option Database Keys:
1491: + -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1492: . -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1493: . -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1494: . -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1495: - -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite
1497: Level: beginner
1498: M*/
1499: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1500: {
1501: TAO_PDIPM *pdipm;
1505: tao->ops->setup = TaoSetup_PDIPM;
1506: tao->ops->solve = TaoSolve_PDIPM;
1507: tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1508: tao->ops->view = TaoView_PDIPM;
1509: tao->ops->destroy = TaoDestroy_PDIPM;
1511: PetscNewLog(tao,&pdipm);
1512: tao->data = (void*)pdipm;
1514: pdipm->nx = pdipm->Nx = 0;
1515: pdipm->nxfixed = pdipm->Nxfixed = 0;
1516: pdipm->nxlb = pdipm->Nxlb = 0;
1517: pdipm->nxub = pdipm->Nxub = 0;
1518: pdipm->nxbox = pdipm->Nxbox = 0;
1519: pdipm->nxfree = pdipm->Nxfree = 0;
1521: pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1522: pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1523: pdipm->n = pdipm->N = 0;
1524: pdipm->mu = 1.0;
1525: pdipm->mu_update_factor = 0.1;
1527: pdipm->deltaw = 0.0;
1528: pdipm->lastdeltaw = 3*1.e-4;
1529: pdipm->deltac = 0.0;
1530: pdipm->kkt_pd = PETSC_FALSE;
1532: pdipm->push_init_slack = 1.0;
1533: pdipm->push_init_lambdai = 1.0;
1534: pdipm->solve_reduced_kkt = PETSC_FALSE;
1535: pdipm->solve_symmetric_kkt = PETSC_TRUE;
1537: /* Override default settings (unless already changed) */
1538: if (!tao->max_it_changed) tao->max_it = 200;
1539: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1541: SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1542: SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1543: SNESGetKSP(pdipm->snes,&tao->ksp);
1544: PetscObjectReference((PetscObject)tao->ksp);
1545: return(0);
1546: }