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