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
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) ISCreateGeneral(comm, pdipm->nxlb, ixlb, PETSC_COPY_VALUES, &pdipm->isxlb);
196: if (pdipm->Nxub) ISCreateGeneral(comm, pdipm->nxub, ixub, PETSC_COPY_VALUES, &pdipm->isxub);
197: if (pdipm->Nxfixed) ISCreateGeneral(comm, pdipm->nxfixed, ixfixed, PETSC_COPY_VALUES, &pdipm->isxfixed);
198: if (pdipm->Nxbox) ISCreateGeneral(comm, pdipm->nxbox, ixbox, PETSC_COPY_VALUES, &pdipm->isxbox);
199: if (pdipm->Nxfree) ISCreateGeneral(comm, pdipm->nxfree, ixfree, PETSC_COPY_VALUES, &pdipm->isxfree);
200: PetscFree5(ixlb, ixub, ixfixed, ixbox, ixfree);
201: return 0;
202: }
204: /*
205: TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
206: X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
207: four subvectors need to be initialized and its values copied over to X. Instead
208: of copying, we use VecPlace/ResetArray functions to share the memory locations for
209: X and the subvectors
211: Collective
213: Input Parameter:
214: . tao - Tao context
216: Level: beginner
217: */
218: static PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
219: {
220: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
221: PetscScalar *Xarr, *z, *lambdai;
222: PetscInt i;
223: const PetscScalar *xarr, *h;
225: VecGetArrayWrite(pdipm->X, &Xarr);
227: /* Set Initialize X.x = tao->solution */
228: VecGetArrayRead(tao->solution, &xarr);
229: PetscMemcpy(Xarr, xarr, pdipm->nx * sizeof(PetscScalar));
230: VecRestoreArrayRead(tao->solution, &xarr);
232: /* Initialize X.lambdae = 0.0 */
233: if (pdipm->lambdae) VecSet(pdipm->lambdae, 0.0);
235: /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
236: if (pdipm->Nci) {
237: VecSet(pdipm->lambdai, pdipm->push_init_lambdai);
238: VecSet(pdipm->z, pdipm->push_init_slack);
240: /* Additional modification for X.lambdai and X.z */
241: VecGetArrayWrite(pdipm->lambdai, &lambdai);
242: VecGetArrayWrite(pdipm->z, &z);
243: if (pdipm->Nh) {
244: VecGetArrayRead(tao->constraints_inequality, &h);
245: for (i = 0; i < pdipm->nh; i++) {
246: if (h[i] < -pdipm->push_init_slack) z[i] = -h[i];
247: if (pdipm->mu / z[i] > pdipm->push_init_lambdai) lambdai[i] = pdipm->mu / z[i];
248: }
249: VecRestoreArrayRead(tao->constraints_inequality, &h);
250: }
251: VecRestoreArrayWrite(pdipm->lambdai, &lambdai);
252: VecRestoreArrayWrite(pdipm->z, &z);
253: }
255: VecRestoreArrayWrite(pdipm->X, &Xarr);
256: return 0;
257: }
259: /*
260: TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X
262: Input Parameter:
263: snes - SNES context
264: X - KKT Vector
265: *ctx - pdipm context
267: Output Parameter:
268: J - Hessian matrix
269: Jpre - Preconditioner
270: */
271: static PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes, Vec X, Mat J, Mat Jpre, void *ctx)
272: {
273: Tao tao = (Tao)ctx;
274: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
275: PetscInt i, row, cols[2], Jrstart, rjstart, nc, j;
276: const PetscInt *aj, *ranges, *Jranges, *rranges, *cranges;
277: const PetscScalar *Xarr, *aa;
278: PetscScalar vals[2];
279: PetscInt proc, nx_all, *nce_all = pdipm->nce_all;
280: MPI_Comm comm;
281: PetscMPIInt rank, size;
283: PetscObjectGetComm((PetscObject)snes, &comm);
284: MPI_Comm_rank(comm, &rank);
285: MPI_Comm_rank(comm, &size);
287: MatGetOwnershipRanges(Jpre, &Jranges);
288: MatGetOwnershipRange(Jpre, &Jrstart, NULL);
289: MatGetOwnershipRangesColumn(tao->hessian, &rranges);
290: MatGetOwnershipRangesColumn(tao->hessian, &cranges);
292: VecGetArrayRead(X, &Xarr);
294: /* (1) insert Z and Ci to the 4th block of Jpre -- overwrite existing values */
295: if (pdipm->solve_symmetric_kkt) { /* 1 for eq 17 revised pdipm doc 0 for eq 18 (symmetric KKT) */
296: vals[0] = 1.0;
297: for (i = 0; i < pdipm->nci; i++) {
298: row = Jrstart + pdipm->off_z + i;
299: cols[0] = Jrstart + pdipm->off_lambdai + i;
300: cols[1] = row;
301: vals[1] = Xarr[pdipm->off_lambdai + i] / Xarr[pdipm->off_z + i];
302: MatSetValues(Jpre, 1, &row, 2, cols, vals, INSERT_VALUES);
303: }
304: } else {
305: for (i = 0; i < pdipm->nci; i++) {
306: row = Jrstart + pdipm->off_z + i;
307: cols[0] = Jrstart + pdipm->off_lambdai + i;
308: cols[1] = row;
309: vals[0] = Xarr[pdipm->off_z + i];
310: vals[1] = Xarr[pdipm->off_lambdai + i];
311: MatSetValues(Jpre, 1, &row, 2, cols, vals, INSERT_VALUES);
312: }
313: }
315: /* (2) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
316: if (pdipm->Ng) {
317: MatGetOwnershipRange(tao->jacobian_equality, &rjstart, NULL);
318: for (i = 0; i < pdipm->ng; i++) {
319: row = Jrstart + pdipm->off_lambdae + i;
321: MatGetRow(tao->jacobian_equality, i + rjstart, &nc, &aj, &aa);
322: proc = 0;
323: for (j = 0; j < nc; j++) {
324: while (aj[j] >= cranges[proc + 1]) proc++;
325: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
326: MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
327: }
328: MatRestoreRow(tao->jacobian_equality, i + rjstart, &nc, &aj, &aa);
329: if (pdipm->kkt_pd) {
330: /* add shift \delta_c */
331: MatSetValue(Jpre, row, row, -pdipm->deltac, INSERT_VALUES);
332: }
333: }
334: }
336: /* (3) insert 3rd row block of Jpre: [ -grad h, 0, deltac, I] */
337: if (pdipm->Nh) {
338: MatGetOwnershipRange(tao->jacobian_inequality, &rjstart, NULL);
339: for (i = 0; i < pdipm->nh; i++) {
340: row = Jrstart + pdipm->off_lambdai + i;
341: MatGetRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, &aa);
342: proc = 0;
343: for (j = 0; j < nc; j++) {
344: while (aj[j] >= cranges[proc + 1]) proc++;
345: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
346: MatSetValue(Jpre, row, cols[0], -aa[j], INSERT_VALUES);
347: }
348: MatRestoreRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, &aa);
349: if (pdipm->kkt_pd) {
350: /* add shift \delta_c */
351: MatSetValue(Jpre, row, row, -pdipm->deltac, INSERT_VALUES);
352: }
353: }
354: }
356: /* (4) insert 1st row block of Jpre: [Wxx, grad g', -grad h', 0] */
357: if (pdipm->Ng) { /* grad g' */
358: MatTranspose(tao->jacobian_equality, MAT_REUSE_MATRIX, &pdipm->jac_equality_trans);
359: }
360: if (pdipm->Nh) { /* grad h' */
361: MatTranspose(tao->jacobian_inequality, MAT_REUSE_MATRIX, &pdipm->jac_inequality_trans);
362: }
364: VecPlaceArray(pdipm->x, Xarr);
365: TaoComputeHessian(tao, pdipm->x, tao->hessian, tao->hessian_pre);
366: VecResetArray(pdipm->x);
368: MatGetOwnershipRange(tao->hessian, &rjstart, NULL);
369: for (i = 0; i < pdipm->nx; i++) {
370: row = Jrstart + i;
372: /* insert Wxx = fxx + ... -- provided by user */
373: MatGetRow(tao->hessian, i + rjstart, &nc, &aj, &aa);
374: proc = 0;
375: for (j = 0; j < nc; j++) {
376: while (aj[j] >= cranges[proc + 1]) proc++;
377: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
378: if (row == cols[0] && pdipm->kkt_pd) {
379: /* add shift deltaw to Wxx component */
380: MatSetValue(Jpre, row, cols[0], aa[j] + pdipm->deltaw, INSERT_VALUES);
381: } else {
382: MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
383: }
384: }
385: MatRestoreRow(tao->hessian, i + rjstart, &nc, &aj, &aa);
387: /* insert grad g' */
388: if (pdipm->ng) {
389: MatGetRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, &aa);
390: MatGetOwnershipRanges(tao->jacobian_equality, &ranges);
391: proc = 0;
392: for (j = 0; j < nc; j++) {
393: /* find row ownership of */
394: while (aj[j] >= ranges[proc + 1]) proc++;
395: nx_all = rranges[proc + 1] - rranges[proc];
396: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
397: MatSetValue(Jpre, row, cols[0], aa[j], INSERT_VALUES);
398: }
399: MatRestoreRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, &aa);
400: }
402: /* insert -grad h' */
403: if (pdipm->nh) {
404: MatGetRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, &aa);
405: MatGetOwnershipRanges(tao->jacobian_inequality, &ranges);
406: proc = 0;
407: for (j = 0; j < nc; j++) {
408: /* find row ownership of */
409: while (aj[j] >= ranges[proc + 1]) proc++;
410: nx_all = rranges[proc + 1] - rranges[proc];
411: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
412: MatSetValue(Jpre, row, cols[0], -aa[j], INSERT_VALUES);
413: }
414: MatRestoreRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, &aa);
415: }
416: }
417: VecRestoreArrayRead(X, &Xarr);
419: /* (6) assemble Jpre and J */
420: MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
421: MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
423: if (J != Jpre) {
424: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
425: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
426: }
427: return 0;
428: }
430: /*
431: TaoSnesFunction_PDIPM - Evaluate KKT function at X
433: Input Parameter:
434: snes - SNES context
435: X - KKT Vector
436: *ctx - pdipm
438: Output Parameter:
439: F - Updated Lagrangian vector
440: */
441: static PetscErrorCode TaoSNESFunction_PDIPM(SNES snes, Vec X, Vec F, void *ctx)
442: {
443: Tao tao = (Tao)ctx;
444: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
445: PetscScalar *Farr;
446: Vec x, L1;
447: PetscInt i;
448: const PetscScalar *Xarr, *carr, *zarr, *larr;
450: VecSet(F, 0.0);
452: VecGetArrayRead(X, &Xarr);
453: VecGetArrayWrite(F, &Farr);
455: /* (0) Evaluate f, fx, gradG, gradH at X.x Note: pdipm->x is not changed below */
456: x = pdipm->x;
457: VecPlaceArray(x, Xarr);
458: TaoPDIPMEvaluateFunctionsAndJacobians(tao, x);
460: /* Update ce, ci, and Jci at X.x */
461: TaoPDIPMUpdateConstraints(tao, x);
462: VecResetArray(x);
464: /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
465: L1 = pdipm->x;
466: VecPlaceArray(L1, Farr); /* L1 = 0.0 */
467: if (pdipm->Nci) {
468: if (pdipm->Nh) {
469: /* L1 += gradH'*DI. Note: tao->DI is not changed below */
470: VecPlaceArray(tao->DI, Xarr + pdipm->off_lambdai);
471: MatMultTransposeAdd(tao->jacobian_inequality, tao->DI, L1, L1);
472: VecResetArray(tao->DI);
473: }
475: /* L1 += Jci_xb'*lambdai_xb */
476: VecPlaceArray(pdipm->lambdai_xb, Xarr + pdipm->off_lambdai + pdipm->nh);
477: MatMultTransposeAdd(pdipm->Jci_xb, pdipm->lambdai_xb, L1, L1);
478: VecResetArray(pdipm->lambdai_xb);
480: /* L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
481: VecScale(L1, -1.0);
482: }
484: /* L1 += fx */
485: VecAXPY(L1, 1.0, tao->gradient);
487: if (pdipm->Nce) {
488: if (pdipm->Ng) {
489: /* L1 += gradG'*DE. Note: tao->DE is not changed below */
490: VecPlaceArray(tao->DE, Xarr + pdipm->off_lambdae);
491: MatMultTransposeAdd(tao->jacobian_equality, tao->DE, L1, L1);
492: VecResetArray(tao->DE);
493: }
494: if (pdipm->Nxfixed) {
495: /* L1 += Jce_xfixed'*lambdae_xfixed */
496: VecPlaceArray(pdipm->lambdae_xfixed, Xarr + pdipm->off_lambdae + pdipm->ng);
497: MatMultTransposeAdd(pdipm->Jce_xfixed, pdipm->lambdae_xfixed, L1, L1);
498: VecResetArray(pdipm->lambdae_xfixed);
499: }
500: }
501: VecResetArray(L1);
503: /* (2) L2 = ce(x) */
504: if (pdipm->Nce) {
505: VecGetArrayRead(pdipm->ce, &carr);
506: for (i = 0; i < pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
507: VecRestoreArrayRead(pdipm->ce, &carr);
508: }
510: if (pdipm->Nci) {
511: if (pdipm->solve_symmetric_kkt) {
512: /* (3) L3 = z - ci(x);
513: (4) L4 = Lambdai * e - mu/z *e */
514: VecGetArrayRead(pdipm->ci, &carr);
515: larr = Xarr + pdipm->off_lambdai;
516: zarr = Xarr + pdipm->off_z;
517: for (i = 0; i < pdipm->nci; i++) {
518: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
519: Farr[pdipm->off_z + i] = larr[i] - pdipm->mu / zarr[i];
520: }
521: VecRestoreArrayRead(pdipm->ci, &carr);
522: } else {
523: /* (3) L3 = z - ci(x);
524: (4) L4 = Z * Lambdai * e - mu * e */
525: VecGetArrayRead(pdipm->ci, &carr);
526: larr = Xarr + pdipm->off_lambdai;
527: zarr = Xarr + pdipm->off_z;
528: for (i = 0; i < pdipm->nci; i++) {
529: Farr[pdipm->off_lambdai + i] = zarr[i] - carr[i];
530: Farr[pdipm->off_z + i] = zarr[i] * larr[i] - pdipm->mu;
531: }
532: VecRestoreArrayRead(pdipm->ci, &carr);
533: }
534: }
536: VecRestoreArrayRead(X, &Xarr);
537: VecRestoreArrayWrite(F, &Farr);
538: return 0;
539: }
541: /*
542: Evaluate F(X); then update update tao->gnorm0, tao->step = mu,
543: tao->residual = norm2(F_x,F_z) and tao->cnorm = norm2(F_ce,F_ci).
544: */
545: static PetscErrorCode TaoSNESFunction_PDIPM_residual(SNES snes, Vec X, Vec F, void *ctx)
546: {
547: Tao tao = (Tao)ctx;
548: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
549: PetscScalar *Farr, *tmparr;
550: Vec L1;
551: PetscInt i;
552: PetscReal res[2], cnorm[2];
553: const PetscScalar *Xarr = NULL;
555: TaoSNESFunction_PDIPM(snes, X, F, (void *)tao);
556: VecGetArrayWrite(F, &Farr);
557: VecGetArrayRead(X, &Xarr);
559: /* compute res[0] = norm2(F_x) */
560: L1 = pdipm->x;
561: VecPlaceArray(L1, Farr);
562: VecNorm(L1, NORM_2, &res[0]);
563: VecResetArray(L1);
565: /* compute res[1] = norm2(F_z), cnorm[1] = norm2(F_ci) */
566: if (pdipm->z) {
567: if (pdipm->solve_symmetric_kkt) {
568: VecPlaceArray(pdipm->z, Farr + pdipm->off_z);
569: if (pdipm->Nci) {
570: VecGetArrayWrite(pdipm->z, &tmparr);
571: for (i = 0; i < pdipm->nci; i++) tmparr[i] *= Xarr[pdipm->off_z + i];
572: VecRestoreArrayWrite(pdipm->z, &tmparr);
573: }
575: VecNorm(pdipm->z, NORM_2, &res[1]);
577: if (pdipm->Nci) {
578: VecGetArrayWrite(pdipm->z, &tmparr);
579: for (i = 0; i < pdipm->nci; i++) tmparr[i] /= Xarr[pdipm->off_z + i];
580: VecRestoreArrayWrite(pdipm->z, &tmparr);
581: }
582: VecResetArray(pdipm->z);
583: } else { /* !solve_symmetric_kkt */
584: VecPlaceArray(pdipm->z, Farr + pdipm->off_z);
585: VecNorm(pdipm->z, NORM_2, &res[1]);
586: VecResetArray(pdipm->z);
587: }
589: VecPlaceArray(pdipm->ci, Farr + pdipm->off_lambdai);
590: VecNorm(pdipm->ci, NORM_2, &cnorm[1]);
591: VecResetArray(pdipm->ci);
592: } else {
593: res[1] = 0.0;
594: cnorm[1] = 0.0;
595: }
597: /* compute cnorm[0] = norm2(F_ce) */
598: if (pdipm->Nce) {
599: VecPlaceArray(pdipm->ce, Farr + pdipm->off_lambdae);
600: VecNorm(pdipm->ce, NORM_2, &cnorm[0]);
601: VecResetArray(pdipm->ce);
602: } else cnorm[0] = 0.0;
604: VecRestoreArrayWrite(F, &Farr);
605: VecRestoreArrayRead(X, &Xarr);
607: tao->gnorm0 = tao->residual;
608: tao->residual = PetscSqrtReal(res[0] * res[0] + res[1] * res[1]);
609: tao->cnorm = PetscSqrtReal(cnorm[0] * cnorm[0] + cnorm[1] * cnorm[1]);
610: tao->step = pdipm->mu;
611: return 0;
612: }
614: /*
615: KKTAddShifts - Check the inertia of Cholesky factor of KKT matrix.
616: If it does not match the numbers of prime and dual variables, add shifts to the KKT matrix.
617: */
618: static PetscErrorCode KKTAddShifts(Tao tao, SNES snes, Vec X)
619: {
620: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
621: KSP ksp;
622: PC pc;
623: Mat Factor;
624: PetscBool isCHOL;
625: PetscInt nneg, nzero, npos;
627: /* Get the inertia of Cholesky factor */
628: SNESGetKSP(snes, &ksp);
629: KSPGetPC(ksp, &pc);
630: PetscObjectTypeCompare((PetscObject)pc, PCCHOLESKY, &isCHOL);
631: if (!isCHOL) return 0;
633: PCFactorGetMatrix(pc, &Factor);
634: MatGetInertia(Factor, &nneg, &nzero, &npos);
636: if (npos < pdipm->Nx + pdipm->Nci) {
637: pdipm->deltaw = PetscMax(pdipm->lastdeltaw / 3, 1.e-4 * PETSC_MACHINE_EPSILON);
638: PetscInfo(tao, "Test reduced deltaw=%g; previous MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT ", npos %" PetscInt_FMT "(<%" PetscInt_FMT ")\n", (double)pdipm->deltaw, nneg, nzero, npos, pdipm->Nx + pdipm->Nci);
639: TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
640: PCSetUp(pc);
641: MatGetInertia(Factor, &nneg, &nzero, &npos);
643: if (npos < pdipm->Nx + pdipm->Nci) {
644: pdipm->deltaw = pdipm->lastdeltaw; /* in case reduction update does not help, this prevents that step from impacting increasing update */
645: while (npos < pdipm->Nx + pdipm->Nci && pdipm->deltaw <= 1. / PETSC_SMALL) { /* increase deltaw */
646: PetscInfo(tao, " deltaw=%g fails, MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT ", npos %" PetscInt_FMT "(<%" PetscInt_FMT ")\n", (double)pdipm->deltaw, nneg, nzero, npos, pdipm->Nx + pdipm->Nci);
647: pdipm->deltaw = PetscMin(8 * pdipm->deltaw, PetscPowReal(10, 20));
648: TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
649: PCSetUp(pc);
650: MatGetInertia(Factor, &nneg, &nzero, &npos);
651: }
655: PetscInfo(tao, "Updated deltaw %g\n", (double)pdipm->deltaw);
656: pdipm->lastdeltaw = pdipm->deltaw;
657: pdipm->deltaw = 0.0;
658: }
659: }
661: if (nzero) { /* Jacobian is singular */
662: if (pdipm->deltac == 0.0) {
663: pdipm->deltac = PETSC_SQRT_MACHINE_EPSILON;
664: } else {
665: pdipm->deltac = pdipm->deltac * PetscPowReal(pdipm->mu, .25);
666: }
667: PetscInfo(tao, "Updated deltac=%g, MatInertia: nneg %" PetscInt_FMT ", nzero %" PetscInt_FMT "(!=0), npos %" PetscInt_FMT "\n", (double)pdipm->deltac, nneg, nzero, npos);
668: TaoSNESJacobian_PDIPM(snes, X, pdipm->K, pdipm->K, tao);
669: PCSetUp(pc);
670: MatGetInertia(Factor, &nneg, &nzero, &npos);
671: }
672: return 0;
673: }
675: /*
676: PCPreSolve_PDIPM -- called between MatFactorNumeric() and MatSolve()
677: */
678: PetscErrorCode PCPreSolve_PDIPM(PC pc, KSP ksp)
679: {
680: Tao tao;
681: TAO_PDIPM *pdipm;
683: KSPGetApplicationContext(ksp, &tao);
684: pdipm = (TAO_PDIPM *)tao->data;
685: KKTAddShifts(tao, pdipm->snes, pdipm->X);
686: return 0;
687: }
689: /*
690: SNESLineSearch_PDIPM - Custom line search used with PDIPM.
692: Collective
694: Notes:
695: This routine employs a simple backtracking line-search to keep
696: the slack variables (z) and inequality constraints Lagrange multipliers
697: (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
698: alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
699: slack variables are updated as X = X - alpha_d*dx. The constraint multipliers
700: are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
701: is also updated as mu = mu + z'lambdai/Nci
702: */
703: static PetscErrorCode SNESLineSearch_PDIPM(SNESLineSearch linesearch, void *ctx)
704: {
705: Tao tao = (Tao)ctx;
706: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
707: SNES snes;
708: Vec X, F, Y;
709: PetscInt i, iter;
710: PetscReal alpha_p = 1.0, alpha_d = 1.0, alpha[4];
711: PetscScalar *Xarr, *z, *lambdai, dot, *taosolarr;
712: const PetscScalar *dXarr, *dz, *dlambdai;
714: SNESLineSearchGetSNES(linesearch, &snes);
715: SNESGetIterationNumber(snes, &iter);
717: SNESLineSearchSetReason(linesearch, SNES_LINESEARCH_SUCCEEDED);
718: SNESLineSearchGetVecs(linesearch, &X, &F, &Y, NULL, NULL);
720: VecGetArrayWrite(X, &Xarr);
721: VecGetArrayRead(Y, &dXarr);
722: z = Xarr + pdipm->off_z;
723: dz = dXarr + pdipm->off_z;
724: for (i = 0; i < pdipm->nci; i++) {
725: if (z[i] - dz[i] < 0.0) alpha_p = PetscMin(alpha_p, 0.9999 * z[i] / dz[i]);
726: }
728: lambdai = Xarr + pdipm->off_lambdai;
729: dlambdai = dXarr + pdipm->off_lambdai;
731: for (i = 0; i < pdipm->nci; i++) {
732: if (lambdai[i] - dlambdai[i] < 0.0) alpha_d = PetscMin(0.9999 * lambdai[i] / dlambdai[i], alpha_d);
733: }
735: alpha[0] = alpha_p;
736: alpha[1] = alpha_d;
737: VecRestoreArrayRead(Y, &dXarr);
738: VecRestoreArrayWrite(X, &Xarr);
740: /* alpha = min(alpha) over all processes */
741: MPI_Allreduce(alpha, alpha + 2, 2, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)tao));
743: alpha_p = alpha[2];
744: alpha_d = alpha[3];
746: /* X = X - alpha * Y */
747: VecGetArrayWrite(X, &Xarr);
748: VecGetArrayRead(Y, &dXarr);
749: for (i = 0; i < pdipm->nx; i++) Xarr[i] -= alpha_p * dXarr[i];
750: for (i = 0; i < pdipm->nce; i++) Xarr[i + pdipm->off_lambdae] -= alpha_d * dXarr[i + pdipm->off_lambdae];
752: for (i = 0; i < pdipm->nci; i++) {
753: Xarr[i + pdipm->off_lambdai] -= alpha_d * dXarr[i + pdipm->off_lambdai];
754: Xarr[i + pdipm->off_z] -= alpha_p * dXarr[i + pdipm->off_z];
755: }
756: VecGetArrayWrite(tao->solution, &taosolarr);
757: PetscMemcpy(taosolarr, Xarr, pdipm->nx * sizeof(PetscScalar));
758: VecRestoreArrayWrite(tao->solution, &taosolarr);
760: VecRestoreArrayWrite(X, &Xarr);
761: VecRestoreArrayRead(Y, &dXarr);
763: /* Update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
764: if (pdipm->z) {
765: VecDot(pdipm->z, pdipm->lambdai, &dot);
766: } else dot = 0.0;
768: /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu) */
769: pdipm->mu = pdipm->mu_update_factor * dot / pdipm->Nci;
771: /* Update F; get tao->residual and tao->cnorm */
772: TaoSNESFunction_PDIPM_residual(snes, X, F, (void *)tao);
774: tao->niter++;
775: TaoLogConvergenceHistory(tao, pdipm->obj, tao->residual, tao->cnorm, tao->niter);
776: TaoMonitor(tao, tao->niter, pdipm->obj, tao->residual, tao->cnorm, pdipm->mu);
778: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
779: if (tao->reason) SNESSetConvergedReason(snes, SNES_CONVERGED_FNORM_ABS);
780: return 0;
781: }
783: /*
784: TaoSolve_PDIPM
786: Input Parameter:
787: tao - TAO context
789: Output Parameter:
790: tao - TAO context
791: */
792: PetscErrorCode TaoSolve_PDIPM(Tao tao)
793: {
794: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
795: SNESLineSearch linesearch; /* SNESLineSearch context */
796: Vec dummy;
800: /* Initialize all variables */
801: TaoPDIPMInitializeSolution(tao);
803: /* Set linesearch */
804: SNESGetLineSearch(pdipm->snes, &linesearch);
805: SNESLineSearchSetType(linesearch, SNESLINESEARCHSHELL);
806: SNESLineSearchShellSetUserFunc(linesearch, SNESLineSearch_PDIPM, tao);
807: SNESLineSearchSetFromOptions(linesearch);
809: tao->reason = TAO_CONTINUE_ITERATING;
811: /* -tao_monitor for iteration 0 and check convergence */
812: VecDuplicate(pdipm->X, &dummy);
813: TaoSNESFunction_PDIPM_residual(pdipm->snes, pdipm->X, dummy, (void *)tao);
815: TaoLogConvergenceHistory(tao, pdipm->obj, tao->residual, tao->cnorm, tao->niter);
816: TaoMonitor(tao, tao->niter, pdipm->obj, tao->residual, tao->cnorm, pdipm->mu);
817: VecDestroy(&dummy);
818: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
819: if (tao->reason) SNESSetConvergedReason(pdipm->snes, SNES_CONVERGED_FNORM_ABS);
821: while (tao->reason == TAO_CONTINUE_ITERATING) {
822: SNESConvergedReason reason;
823: SNESSolve(pdipm->snes, NULL, pdipm->X);
825: /* Check SNES convergence */
826: SNESGetConvergedReason(pdipm->snes, &reason);
827: if (reason < 0) PetscPrintf(PetscObjectComm((PetscObject)pdipm->snes), "SNES solve did not converged due to reason %s\n", SNESConvergedReasons[reason]);
829: /* Check TAO convergence */
831: }
832: return 0;
833: }
835: /*
836: TaoView_PDIPM - View PDIPM
838: Input Parameter:
839: tao - TAO object
840: viewer - PetscViewer
842: Output:
843: */
844: PetscErrorCode TaoView_PDIPM(Tao tao, PetscViewer viewer)
845: {
846: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
848: tao->constrained = PETSC_TRUE;
849: PetscViewerASCIIPushTab(viewer);
850: PetscViewerASCIIPrintf(viewer, "Number of prime=%" PetscInt_FMT ", Number of dual=%" PetscInt_FMT "\n", pdipm->Nx + pdipm->Nci, pdipm->Nce + pdipm->Nci);
851: if (pdipm->kkt_pd) PetscViewerASCIIPrintf(viewer, "KKT shifts deltaw=%g, deltac=%g\n", (double)pdipm->deltaw, (double)pdipm->deltac);
852: PetscViewerASCIIPopTab(viewer);
853: return 0;
854: }
856: /*
857: TaoSetup_PDIPM - Sets up tao and pdipm
859: Input Parameter:
860: tao - TAO object
862: Output: pdipm - initialized object
863: */
864: PetscErrorCode TaoSetup_PDIPM(Tao tao)
865: {
866: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
867: MPI_Comm comm;
868: PetscMPIInt size;
869: PetscInt row, col, Jcrstart, Jcrend, k, tmp, nc, proc, *nh_all, *ng_all;
870: PetscInt offset, *xa, *xb, i, j, rstart, rend;
871: PetscScalar one = 1.0, neg_one = -1.0;
872: const PetscInt *cols, *rranges, *cranges, *aj, *ranges;
873: const PetscScalar *aa, *Xarr;
874: Mat J;
875: Mat Jce_xfixed_trans, Jci_xb_trans;
876: PetscInt *dnz, *onz, rjstart, nx_all, *nce_all, *Jranges, cols1[2];
878: PetscObjectGetComm((PetscObject)tao, &comm);
879: MPI_Comm_size(comm, &size);
881: /* (1) Setup Bounds and create Tao vectors */
882: TaoPDIPMSetUpBounds(tao);
884: if (!tao->gradient) {
885: VecDuplicate(tao->solution, &tao->gradient);
886: VecDuplicate(tao->solution, &tao->stepdirection);
887: }
889: /* (2) Get sizes */
890: /* Size of vector x - This is set by TaoSetSolution */
891: VecGetSize(tao->solution, &pdipm->Nx);
892: VecGetLocalSize(tao->solution, &pdipm->nx);
894: /* Size of equality constraints and vectors */
895: if (tao->constraints_equality) {
896: VecGetSize(tao->constraints_equality, &pdipm->Ng);
897: VecGetLocalSize(tao->constraints_equality, &pdipm->ng);
898: } else {
899: pdipm->ng = pdipm->Ng = 0;
900: }
902: pdipm->nce = pdipm->ng + pdipm->nxfixed;
903: pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;
905: /* Size of inequality constraints and vectors */
906: if (tao->constraints_inequality) {
907: VecGetSize(tao->constraints_inequality, &pdipm->Nh);
908: VecGetLocalSize(tao->constraints_inequality, &pdipm->nh);
909: } else {
910: pdipm->nh = pdipm->Nh = 0;
911: }
913: pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2 * pdipm->nxbox;
914: pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2 * pdipm->Nxbox;
916: /* Full size of the KKT system to be solved */
917: pdipm->n = pdipm->nx + pdipm->nce + 2 * pdipm->nci;
918: pdipm->N = pdipm->Nx + pdipm->Nce + 2 * pdipm->Nci;
920: /* (3) Offsets for subvectors */
921: pdipm->off_lambdae = pdipm->nx;
922: pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
923: pdipm->off_z = pdipm->off_lambdai + pdipm->nci;
925: /* (4) Create vectors and subvectors */
926: /* Ce and Ci vectors */
927: VecCreate(comm, &pdipm->ce);
928: VecSetSizes(pdipm->ce, pdipm->nce, pdipm->Nce);
929: VecSetFromOptions(pdipm->ce);
931: VecCreate(comm, &pdipm->ci);
932: VecSetSizes(pdipm->ci, pdipm->nci, pdipm->Nci);
933: VecSetFromOptions(pdipm->ci);
935: /* X=[x; lambdae; lambdai; z] for the big KKT system */
936: VecCreate(comm, &pdipm->X);
937: VecSetSizes(pdipm->X, pdipm->n, pdipm->N);
938: VecSetFromOptions(pdipm->X);
940: /* Subvectors; they share local arrays with X */
941: VecGetArrayRead(pdipm->X, &Xarr);
942: /* x shares local array with X.x */
943: if (pdipm->Nx) VecCreateMPIWithArray(comm, 1, pdipm->nx, pdipm->Nx, Xarr, &pdipm->x);
945: /* lambdae shares local array with X.lambdae */
946: if (pdipm->Nce) VecCreateMPIWithArray(comm, 1, pdipm->nce, pdipm->Nce, Xarr + pdipm->off_lambdae, &pdipm->lambdae);
948: /* tao->DE shares local array with X.lambdae_g */
949: if (pdipm->Ng) {
950: VecCreateMPIWithArray(comm, 1, pdipm->ng, pdipm->Ng, Xarr + pdipm->off_lambdae, &tao->DE);
952: VecCreate(comm, &pdipm->lambdae_xfixed);
953: VecSetSizes(pdipm->lambdae_xfixed, pdipm->nxfixed, PETSC_DECIDE);
954: VecSetFromOptions(pdipm->lambdae_xfixed);
955: }
957: if (pdipm->Nci) {
958: /* lambdai shares local array with X.lambdai */
959: VecCreateMPIWithArray(comm, 1, pdipm->nci, pdipm->Nci, Xarr + pdipm->off_lambdai, &pdipm->lambdai);
961: /* z for slack variables; it shares local array with X.z */
962: VecCreateMPIWithArray(comm, 1, pdipm->nci, pdipm->Nci, Xarr + pdipm->off_z, &pdipm->z);
963: }
965: /* tao->DI which shares local array with X.lambdai_h */
966: if (pdipm->Nh) VecCreateMPIWithArray(comm, 1, pdipm->nh, pdipm->Nh, Xarr + pdipm->off_lambdai, &tao->DI);
967: VecCreate(comm, &pdipm->lambdai_xb);
968: VecSetSizes(pdipm->lambdai_xb, (pdipm->nci - pdipm->nh), PETSC_DECIDE);
969: VecSetFromOptions(pdipm->lambdai_xb);
971: VecRestoreArrayRead(pdipm->X, &Xarr);
973: /* (5) Create Jacobians Jce_xfixed and Jci */
974: /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
975: if (pdipm->Nxfixed) {
976: /* Create Jce_xfixed */
977: MatCreate(comm, &pdipm->Jce_xfixed);
978: MatSetSizes(pdipm->Jce_xfixed, pdipm->nxfixed, pdipm->nx, PETSC_DECIDE, pdipm->Nx);
979: MatSetFromOptions(pdipm->Jce_xfixed);
980: MatSeqAIJSetPreallocation(pdipm->Jce_xfixed, 1, NULL);
981: MatMPIAIJSetPreallocation(pdipm->Jce_xfixed, 1, NULL, 1, NULL);
983: MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, &Jcrend);
984: ISGetIndices(pdipm->isxfixed, &cols);
985: k = 0;
986: for (row = Jcrstart; row < Jcrend; row++) {
987: MatSetValues(pdipm->Jce_xfixed, 1, &row, 1, cols + k, &one, INSERT_VALUES);
988: k++;
989: }
990: ISRestoreIndices(pdipm->isxfixed, &cols);
991: MatAssemblyBegin(pdipm->Jce_xfixed, MAT_FINAL_ASSEMBLY);
992: MatAssemblyEnd(pdipm->Jce_xfixed, MAT_FINAL_ASSEMBLY);
993: }
995: /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
996: MatCreate(comm, &pdipm->Jci_xb);
997: MatSetSizes(pdipm->Jci_xb, pdipm->nci - pdipm->nh, pdipm->nx, PETSC_DECIDE, pdipm->Nx);
998: MatSetFromOptions(pdipm->Jci_xb);
999: MatSeqAIJSetPreallocation(pdipm->Jci_xb, 1, NULL);
1000: MatMPIAIJSetPreallocation(pdipm->Jci_xb, 1, NULL, 1, NULL);
1002: MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, &Jcrend);
1003: offset = Jcrstart;
1004: if (pdipm->Nxub) {
1005: /* Add xub to Jci_xb */
1006: ISGetIndices(pdipm->isxub, &cols);
1007: k = 0;
1008: for (row = offset; row < offset + pdipm->nxub; row++) {
1009: MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &neg_one, INSERT_VALUES);
1010: k++;
1011: }
1012: ISRestoreIndices(pdipm->isxub, &cols);
1013: }
1015: if (pdipm->Nxlb) {
1016: /* Add xlb to Jci_xb */
1017: ISGetIndices(pdipm->isxlb, &cols);
1018: k = 0;
1019: offset += pdipm->nxub;
1020: for (row = offset; row < offset + pdipm->nxlb; row++) {
1021: MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &one, INSERT_VALUES);
1022: k++;
1023: }
1024: ISRestoreIndices(pdipm->isxlb, &cols);
1025: }
1027: /* Add xbox to Jci_xb */
1028: if (pdipm->Nxbox) {
1029: ISGetIndices(pdipm->isxbox, &cols);
1030: k = 0;
1031: offset += pdipm->nxlb;
1032: for (row = offset; row < offset + pdipm->nxbox; row++) {
1033: MatSetValues(pdipm->Jci_xb, 1, &row, 1, cols + k, &neg_one, INSERT_VALUES);
1034: tmp = row + pdipm->nxbox;
1035: MatSetValues(pdipm->Jci_xb, 1, &tmp, 1, cols + k, &one, INSERT_VALUES);
1036: k++;
1037: }
1038: ISRestoreIndices(pdipm->isxbox, &cols);
1039: }
1041: MatAssemblyBegin(pdipm->Jci_xb, MAT_FINAL_ASSEMBLY);
1042: MatAssemblyEnd(pdipm->Jci_xb, MAT_FINAL_ASSEMBLY);
1043: /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */
1045: /* (6) Set up ISs for PC Fieldsplit */
1046: if (pdipm->solve_reduced_kkt) {
1047: PetscMalloc2(pdipm->nx + pdipm->nce, &xa, 2 * pdipm->nci, &xb);
1048: for (i = 0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
1049: for (i = 0; i < 2 * pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;
1051: ISCreateGeneral(comm, pdipm->nx + pdipm->nce, xa, PETSC_OWN_POINTER, &pdipm->is1);
1052: ISCreateGeneral(comm, 2 * pdipm->nci, xb, PETSC_OWN_POINTER, &pdipm->is2);
1053: }
1055: /* (7) Gather offsets from all processes */
1056: PetscMalloc1(size, &pdipm->nce_all);
1058: /* Get rstart of KKT matrix */
1059: MPI_Scan(&pdipm->n, &rstart, 1, MPIU_INT, MPI_SUM, comm);
1060: rstart -= pdipm->n;
1062: MPI_Allgather(&pdipm->nce, 1, MPIU_INT, pdipm->nce_all, 1, MPIU_INT, comm);
1064: PetscMalloc3(size, &ng_all, size, &nh_all, size, &Jranges);
1065: MPI_Allgather(&rstart, 1, MPIU_INT, Jranges, 1, MPIU_INT, comm);
1066: MPI_Allgather(&pdipm->nh, 1, MPIU_INT, nh_all, 1, MPIU_INT, comm);
1067: MPI_Allgather(&pdipm->ng, 1, MPIU_INT, ng_all, 1, MPIU_INT, comm);
1069: MatGetOwnershipRanges(tao->hessian, &rranges);
1070: MatGetOwnershipRangesColumn(tao->hessian, &cranges);
1072: if (pdipm->Ng) {
1073: TaoComputeJacobianEquality(tao, tao->solution, tao->jacobian_equality, tao->jacobian_equality_pre);
1074: MatTranspose(tao->jacobian_equality, MAT_INITIAL_MATRIX, &pdipm->jac_equality_trans);
1075: }
1076: if (pdipm->Nh) {
1077: TaoComputeJacobianInequality(tao, tao->solution, tao->jacobian_inequality, tao->jacobian_inequality_pre);
1078: MatTranspose(tao->jacobian_inequality, MAT_INITIAL_MATRIX, &pdipm->jac_inequality_trans);
1079: }
1081: /* Count dnz,onz for preallocation of KKT matrix */
1082: nce_all = pdipm->nce_all;
1084: if (pdipm->Nxfixed) MatTranspose(pdipm->Jce_xfixed, MAT_INITIAL_MATRIX, &Jce_xfixed_trans);
1085: MatTranspose(pdipm->Jci_xb, MAT_INITIAL_MATRIX, &Jci_xb_trans);
1087: MatPreallocateBegin(comm, pdipm->n, pdipm->n, dnz, onz);
1089: /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
1090: TaoPDIPMEvaluateFunctionsAndJacobians(tao, pdipm->x);
1091: TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);
1093: /* Insert tao->hessian */
1094: MatGetOwnershipRange(tao->hessian, &rjstart, NULL);
1095: for (i = 0; i < pdipm->nx; i++) {
1096: row = rstart + i;
1098: MatGetRow(tao->hessian, i + rjstart, &nc, &aj, NULL);
1099: proc = 0;
1100: for (j = 0; j < nc; j++) {
1101: while (aj[j] >= cranges[proc + 1]) proc++;
1102: col = aj[j] - cranges[proc] + Jranges[proc];
1103: MatPreallocateSet(row, 1, &col, dnz, onz);
1104: }
1105: MatRestoreRow(tao->hessian, i + rjstart, &nc, &aj, NULL);
1107: if (pdipm->ng) {
1108: /* Insert grad g' */
1109: MatGetRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, NULL);
1110: MatGetOwnershipRanges(tao->jacobian_equality, &ranges);
1111: proc = 0;
1112: for (j = 0; j < nc; j++) {
1113: /* find row ownership of */
1114: while (aj[j] >= ranges[proc + 1]) proc++;
1115: nx_all = rranges[proc + 1] - rranges[proc];
1116: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
1117: MatPreallocateSet(row, 1, &col, dnz, onz);
1118: }
1119: MatRestoreRow(pdipm->jac_equality_trans, i + rjstart, &nc, &aj, NULL);
1120: }
1122: /* Insert Jce_xfixed^T' */
1123: if (pdipm->nxfixed) {
1124: MatGetRow(Jce_xfixed_trans, i + rjstart, &nc, &aj, NULL);
1125: MatGetOwnershipRanges(pdipm->Jce_xfixed, &ranges);
1126: proc = 0;
1127: for (j = 0; j < nc; j++) {
1128: /* find row ownership of */
1129: while (aj[j] >= ranges[proc + 1]) proc++;
1130: nx_all = rranges[proc + 1] - rranges[proc];
1131: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1132: MatPreallocateSet(row, 1, &col, dnz, onz);
1133: }
1134: MatRestoreRow(Jce_xfixed_trans, i + rjstart, &nc, &aj, NULL);
1135: }
1137: if (pdipm->nh) {
1138: /* Insert -grad h' */
1139: MatGetRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, NULL);
1140: MatGetOwnershipRanges(tao->jacobian_inequality, &ranges);
1141: proc = 0;
1142: for (j = 0; j < nc; j++) {
1143: /* find row ownership of */
1144: while (aj[j] >= ranges[proc + 1]) proc++;
1145: nx_all = rranges[proc + 1] - rranges[proc];
1146: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1147: MatPreallocateSet(row, 1, &col, dnz, onz);
1148: }
1149: MatRestoreRow(pdipm->jac_inequality_trans, i + rjstart, &nc, &aj, NULL);
1150: }
1152: /* Insert Jci_xb^T' */
1153: MatGetRow(Jci_xb_trans, i + rjstart, &nc, &aj, NULL);
1154: MatGetOwnershipRanges(pdipm->Jci_xb, &ranges);
1155: proc = 0;
1156: for (j = 0; j < nc; j++) {
1157: /* find row ownership of */
1158: while (aj[j] >= ranges[proc + 1]) proc++;
1159: nx_all = rranges[proc + 1] - rranges[proc];
1160: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1161: MatPreallocateSet(row, 1, &col, dnz, onz);
1162: }
1163: MatRestoreRow(Jci_xb_trans, i + rjstart, &nc, &aj, NULL);
1164: }
1166: /* 2nd Row block of KKT matrix: [grad Ce, deltac*I, 0, 0] */
1167: if (pdipm->Ng) {
1168: MatGetOwnershipRange(tao->jacobian_equality, &rjstart, NULL);
1169: for (i = 0; i < pdipm->ng; i++) {
1170: row = rstart + pdipm->off_lambdae + i;
1172: MatGetRow(tao->jacobian_equality, i + rjstart, &nc, &aj, NULL);
1173: proc = 0;
1174: for (j = 0; j < nc; j++) {
1175: while (aj[j] >= cranges[proc + 1]) proc++;
1176: col = aj[j] - cranges[proc] + Jranges[proc];
1177: MatPreallocateSet(row, 1, &col, dnz, onz); /* grad g */
1178: }
1179: MatRestoreRow(tao->jacobian_equality, i + rjstart, &nc, &aj, NULL);
1180: }
1181: }
1182: /* Jce_xfixed */
1183: if (pdipm->Nxfixed) {
1184: MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, NULL);
1185: for (i = 0; i < (pdipm->nce - pdipm->ng); i++) {
1186: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1188: MatGetRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, NULL);
1191: proc = 0;
1192: j = 0;
1193: while (cols[j] >= cranges[proc + 1]) proc++;
1194: col = cols[j] - cranges[proc] + Jranges[proc];
1195: MatPreallocateSet(row, 1, &col, dnz, onz);
1196: MatRestoreRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, NULL);
1197: }
1198: }
1200: /* 3rd Row block of KKT matrix: [ gradCi, 0, deltac*I, -I] */
1201: if (pdipm->Nh) {
1202: MatGetOwnershipRange(tao->jacobian_inequality, &rjstart, NULL);
1203: for (i = 0; i < pdipm->nh; i++) {
1204: row = rstart + pdipm->off_lambdai + i;
1206: MatGetRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, NULL);
1207: proc = 0;
1208: for (j = 0; j < nc; j++) {
1209: while (aj[j] >= cranges[proc + 1]) proc++;
1210: col = aj[j] - cranges[proc] + Jranges[proc];
1211: MatPreallocateSet(row, 1, &col, dnz, onz); /* grad h */
1212: }
1213: MatRestoreRow(tao->jacobian_inequality, i + rjstart, &nc, &aj, NULL);
1214: }
1215: /* I */
1216: for (i = 0; i < pdipm->nh; i++) {
1217: row = rstart + pdipm->off_lambdai + i;
1218: col = rstart + pdipm->off_z + i;
1219: MatPreallocateSet(row, 1, &col, dnz, onz);
1220: }
1221: }
1223: /* Jci_xb */
1224: MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, NULL);
1225: for (i = 0; i < (pdipm->nci - pdipm->nh); i++) {
1226: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1228: MatGetRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, NULL);
1230: proc = 0;
1231: for (j = 0; j < nc; j++) {
1232: while (cols[j] >= cranges[proc + 1]) proc++;
1233: col = cols[j] - cranges[proc] + Jranges[proc];
1234: MatPreallocateSet(row, 1, &col, dnz, onz);
1235: }
1236: MatRestoreRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, NULL);
1237: /* I */
1238: col = rstart + pdipm->off_z + pdipm->nh + i;
1239: MatPreallocateSet(row, 1, &col, dnz, onz);
1240: }
1242: /* 4-th Row block of KKT matrix: Z and Ci */
1243: for (i = 0; i < pdipm->nci; i++) {
1244: row = rstart + pdipm->off_z + i;
1245: cols1[0] = rstart + pdipm->off_lambdai + i;
1246: cols1[1] = row;
1247: MatPreallocateSet(row, 2, cols1, dnz, onz);
1248: }
1250: /* diagonal entry */
1251: for (i = 0; i < pdipm->n; i++) dnz[i]++; /* diagonal entry */
1253: /* Create KKT matrix */
1254: MatCreate(comm, &J);
1255: MatSetSizes(J, pdipm->n, pdipm->n, PETSC_DECIDE, PETSC_DECIDE);
1256: MatSetFromOptions(J);
1257: MatSeqAIJSetPreallocation(J, 0, dnz);
1258: MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz);
1259: MatPreallocateEnd(dnz, onz);
1260: pdipm->K = J;
1262: /* (8) Insert constant entries to K */
1263: /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1264: MatGetOwnershipRange(J, &rstart, &rend);
1265: for (i = rstart; i < rend; i++) MatSetValue(J, i, i, 0.0, INSERT_VALUES);
1266: /* In case Wxx has no diagonal entries preset set diagonal to deltaw given */
1267: if (pdipm->kkt_pd) {
1268: for (i = 0; i < pdipm->nh; i++) {
1269: row = rstart + i;
1270: MatSetValue(J, row, row, pdipm->deltaw, INSERT_VALUES);
1271: }
1272: }
1274: /* Row block of K: [ grad Ce, 0, 0, 0] */
1275: if (pdipm->Nxfixed) {
1276: MatGetOwnershipRange(pdipm->Jce_xfixed, &Jcrstart, NULL);
1277: for (i = 0; i < (pdipm->nce - pdipm->ng); i++) {
1278: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1280: MatGetRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, &aa);
1281: proc = 0;
1282: for (j = 0; j < nc; j++) {
1283: while (cols[j] >= cranges[proc + 1]) proc++;
1284: col = cols[j] - cranges[proc] + Jranges[proc];
1285: MatSetValue(J, row, col, aa[j], INSERT_VALUES); /* grad Ce */
1286: MatSetValue(J, col, row, aa[j], INSERT_VALUES); /* grad Ce' */
1287: }
1288: MatRestoreRow(pdipm->Jce_xfixed, i + Jcrstart, &nc, &cols, &aa);
1289: }
1290: }
1292: /* Row block of K: [ -grad Ci, 0, 0, I] */
1293: MatGetOwnershipRange(pdipm->Jci_xb, &Jcrstart, NULL);
1294: for (i = 0; i < pdipm->nci - pdipm->nh; i++) {
1295: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1297: MatGetRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, &aa);
1298: proc = 0;
1299: for (j = 0; j < nc; j++) {
1300: while (cols[j] >= cranges[proc + 1]) proc++;
1301: col = cols[j] - cranges[proc] + Jranges[proc];
1302: MatSetValue(J, col, row, -aa[j], INSERT_VALUES);
1303: MatSetValue(J, row, col, -aa[j], INSERT_VALUES);
1304: }
1305: MatRestoreRow(pdipm->Jci_xb, i + Jcrstart, &nc, &cols, &aa);
1307: col = rstart + pdipm->off_z + pdipm->nh + i;
1308: MatSetValue(J, row, col, 1, INSERT_VALUES);
1309: }
1311: for (i = 0; i < pdipm->nh; i++) {
1312: row = rstart + pdipm->off_lambdai + i;
1313: col = rstart + pdipm->off_z + i;
1314: MatSetValue(J, row, col, 1, INSERT_VALUES);
1315: }
1317: /* Row block of K: [ 0, 0, I, ...] */
1318: for (i = 0; i < pdipm->nci; i++) {
1319: row = rstart + pdipm->off_z + i;
1320: col = rstart + pdipm->off_lambdai + i;
1321: MatSetValue(J, row, col, 1, INSERT_VALUES);
1322: }
1324: if (pdipm->Nxfixed) MatDestroy(&Jce_xfixed_trans);
1325: MatDestroy(&Jci_xb_trans);
1326: PetscFree3(ng_all, nh_all, Jranges);
1328: /* (9) Set up nonlinear solver SNES */
1329: SNESSetFunction(pdipm->snes, NULL, TaoSNESFunction_PDIPM, (void *)tao);
1330: SNESSetJacobian(pdipm->snes, J, J, TaoSNESJacobian_PDIPM, (void *)tao);
1332: if (pdipm->solve_reduced_kkt) {
1333: PC pc;
1334: KSPGetPC(tao->ksp, &pc);
1335: PCSetType(pc, PCFIELDSPLIT);
1336: PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1337: PCFieldSplitSetIS(pc, "2", pdipm->is2);
1338: PCFieldSplitSetIS(pc, "1", pdipm->is1);
1339: }
1340: SNESSetFromOptions(pdipm->snes);
1342: /* (10) Setup PCPreSolve() for pdipm->solve_symmetric_kkt */
1343: if (pdipm->solve_symmetric_kkt) {
1344: KSP ksp;
1345: PC pc;
1346: PetscBool isCHOL;
1347: SNESGetKSP(pdipm->snes, &ksp);
1348: KSPGetPC(ksp, &pc);
1349: PCSetPreSolve(pc, PCPreSolve_PDIPM);
1351: PetscObjectTypeCompare((PetscObject)pc, PCCHOLESKY, &isCHOL);
1352: if (isCHOL) {
1353: Mat Factor;
1354: PetscBool isMUMPS;
1355: PCFactorGetMatrix(pc, &Factor);
1356: PetscObjectTypeCompare((PetscObject)Factor, "mumps", &isMUMPS);
1357: if (isMUMPS) { /* must set mumps ICNTL(13)=1 and ICNTL(24)=1 to call MatGetInertia() */
1358: #if defined(PETSC_HAVE_MUMPS)
1359: MatMumpsSetIcntl(Factor, 24, 1); /* detection of null pivot rows */
1360: if (size > 1) { MatMumpsSetIcntl(Factor, 13, 1); /* parallelism of the root node (enable ScaLAPACK) and its splitting */ }
1361: #else
1362: SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Requires external package MUMPS");
1363: #endif
1364: }
1365: }
1366: }
1367: return 0;
1368: }
1370: /*
1371: TaoDestroy_PDIPM - Destroys the pdipm object
1373: Input:
1374: full pdipm
1376: Output:
1377: Destroyed pdipm
1378: */
1379: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1380: {
1381: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
1383: /* Freeing Vectors assocaiated with KKT (X) */
1384: VecDestroy(&pdipm->x); /* Solution x */
1385: VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1386: VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1387: VecDestroy(&pdipm->z); /* Slack variables */
1388: VecDestroy(&pdipm->X); /* Big KKT system vector [x; lambdae; lambdai; z] */
1390: /* work vectors */
1391: VecDestroy(&pdipm->lambdae_xfixed);
1392: VecDestroy(&pdipm->lambdai_xb);
1394: /* Legrangian equality and inequality Vec */
1395: VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1396: VecDestroy(&pdipm->ci); /* Vec of inequality constraints */
1398: /* Matrices */
1399: MatDestroy(&pdipm->Jce_xfixed);
1400: MatDestroy(&pdipm->Jci_xb)); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb; J(nxbx)] */
1401: MatDestroy(&pdipm->K);
1403: /* Index Sets */
1404: if (pdipm->Nxub) { ISDestroy(&pdipm->isxub); /* Finite upper bound only -inf < x < ub */ }
1406: if (pdipm->Nxlb) { ISDestroy(&pdipm->isxlb); /* Finite lower bound only lb <= x < inf */ }
1408: if (pdipm->Nxfixed) { ISDestroy(&pdipm->isxfixed); /* Fixed variables lb = x = ub */ }
1410: if (pdipm->Nxbox) { ISDestroy(&pdipm->isxbox); /* Boxed variables lb <= x <= ub */ }
1412: if (pdipm->Nxfree) { ISDestroy(&pdipm->isxfree); /* Free variables -inf <= x <= inf */ }
1414: if (pdipm->solve_reduced_kkt) {
1415: ISDestroy(&pdipm->is1);
1416: ISDestroy(&pdipm->is2);
1417: }
1419: /* SNES */
1420: SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1421: PetscFree(pdipm->nce_all);
1422: MatDestroy(&pdipm->jac_equality_trans);
1423: MatDestroy(&pdipm->jac_inequality_trans);
1425: /* Destroy pdipm */
1426: PetscFree(tao->data); /* Holding locations of pdipm */
1428: /* Destroy Dual */
1429: VecDestroy(&tao->DE); /* equality dual */
1430: VecDestroy(&tao->DI); /* dinequality dual */
1431: return 0;
1432: }
1434: PetscErrorCode TaoSetFromOptions_PDIPM(Tao tao, PetscOptionItems *PetscOptionsObject)
1435: {
1436: TAO_PDIPM *pdipm = (TAO_PDIPM *)tao->data;
1438: PetscOptionsHeadBegin(PetscOptionsObject, "PDIPM method for constrained optimization");
1439: 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);
1440: 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);
1441: PetscOptionsBool("-tao_pdipm_solve_reduced_kkt", "Solve reduced KKT system using Schur-complement", NULL, pdipm->solve_reduced_kkt, &pdipm->solve_reduced_kkt, NULL);
1442: PetscOptionsReal("-tao_pdipm_mu_update_factor", "Update scalar for barrier parameter (mu) update", NULL, pdipm->mu_update_factor, &pdipm->mu_update_factor, NULL);
1443: PetscOptionsBool("-tao_pdipm_symmetric_kkt", "Solve non reduced symmetric KKT system", NULL, pdipm->solve_symmetric_kkt, &pdipm->solve_symmetric_kkt, NULL);
1444: PetscOptionsBool("-tao_pdipm_kkt_shift_pd", "Add shifts to make KKT matrix positive definite", NULL, pdipm->kkt_pd, &pdipm->kkt_pd, NULL);
1445: PetscOptionsHeadEnd();
1446: return 0;
1447: }
1449: /*MC
1450: TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.
1452: Option Database Keys:
1453: + -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1454: . -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1455: . -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1456: . -tao_pdipm_symmetric_kkt - Solve non-reduced symmetric KKT system
1457: - -tao_pdipm_kkt_shift_pd - Add shifts to make KKT matrix positive definite
1459: Level: beginner
1460: M*/
1461: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1462: {
1463: TAO_PDIPM *pdipm;
1465: tao->ops->setup = TaoSetup_PDIPM;
1466: tao->ops->solve = TaoSolve_PDIPM;
1467: tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1468: tao->ops->view = TaoView_PDIPM;
1469: tao->ops->destroy = TaoDestroy_PDIPM;
1471: PetscNew(&pdipm);
1472: tao->data = (void *)pdipm;
1474: pdipm->nx = pdipm->Nx = 0;
1475: pdipm->nxfixed = pdipm->Nxfixed = 0;
1476: pdipm->nxlb = pdipm->Nxlb = 0;
1477: pdipm->nxub = pdipm->Nxub = 0;
1478: pdipm->nxbox = pdipm->Nxbox = 0;
1479: pdipm->nxfree = pdipm->Nxfree = 0;
1481: pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1482: pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1483: pdipm->n = pdipm->N = 0;
1484: pdipm->mu = 1.0;
1485: pdipm->mu_update_factor = 0.1;
1487: pdipm->deltaw = 0.0;
1488: pdipm->lastdeltaw = 3 * 1.e-4;
1489: pdipm->deltac = 0.0;
1490: pdipm->kkt_pd = PETSC_FALSE;
1492: pdipm->push_init_slack = 1.0;
1493: pdipm->push_init_lambdai = 1.0;
1494: pdipm->solve_reduced_kkt = PETSC_FALSE;
1495: pdipm->solve_symmetric_kkt = PETSC_TRUE;
1497: /* Override default settings (unless already changed) */
1498: if (!tao->max_it_changed) tao->max_it = 200;
1499: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1501: SNESCreate(((PetscObject)tao)->comm, &pdipm->snes);
1502: SNESSetOptionsPrefix(pdipm->snes, tao->hdr.prefix);
1503: SNESGetKSP(pdipm->snes, &tao->ksp);
1504: PetscObjectReference((PetscObject)tao->ksp);
1505: KSPSetApplicationContext(tao->ksp, (void *)tao);
1506: return 0;
1507: }