Actual source code: pdipm.c
petsc-3.13.6 2020-09-29
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/constrained/impls/ipm/pdipm.h>
3: #include <petscsnes.h>
5: /*
6: TaoPDIPMEvaluateFunctionsAndJacobians - Evaluate the objective function f, gradient fx, constraints, and all the Jacobians at current vector
8: Collective on tao
10: Input Parameter:
11: + tao - solver context
12: - x - vector at which all objects to be evaluated
14: Level: beginner
16: .seealso: TaoPDIPMUpdateConstraints(), TaoPDIPMSetUpBounds()
17: */
18: PetscErrorCode TaoPDIPMEvaluateFunctionsAndJacobians(Tao tao,Vec x)
19: {
21: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
24: /* Compute user objective function and gradient */
25: TaoComputeObjectiveAndGradient(tao,x,&pdipm->obj,tao->gradient);
27: /* Equality constraints and Jacobian */
28: if (pdipm->Ng) {
29: TaoComputeEqualityConstraints(tao,x,tao->constraints_equality);
30: TaoComputeJacobianEquality(tao,x,tao->jacobian_equality,tao->jacobian_equality_pre);
31: }
33: /* Inequality constraints and Jacobian */
34: if (pdipm->Nh) {
35: TaoComputeInequalityConstraints(tao,x,tao->constraints_inequality);
36: TaoComputeJacobianInequality(tao,x,tao->jacobian_inequality,tao->jacobian_inequality_pre);
37: }
38: return(0);
39: }
41: /*
42: TaoPDIPMUpdateConstraints - Update the vectors ce and ci at x
44: Collective
46: Input Parameter:
47: + tao - Tao context
48: - x - vector at which constraints to be evaluted
50: Level: beginner
52: .seealso: TaoPDIPMEvaluateFunctionsAndJacobians()
53: */
54: PetscErrorCode TaoPDIPMUpdateConstraints(Tao tao,Vec x)
55: {
56: PetscErrorCode ierr;
57: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
58: PetscInt i,offset,offset1,k,xstart;
59: PetscScalar *carr;
60: const PetscInt *ubptr,*lbptr,*bxptr,*fxptr;
61: const PetscScalar *xarr,*xuarr,*xlarr,*garr,*harr;
64: VecGetOwnershipRange(x,&xstart,NULL);
66: VecGetArrayRead(x,&xarr);
67: VecGetArrayRead(tao->XU,&xuarr);
68: VecGetArrayRead(tao->XL,&xlarr);
70: /* (1) Update ce vector */
71: VecGetArray(pdipm->ce,&carr);
73: /* (1.a) Inserting updated g(x) */
74: if (pdipm->Ng) {
75: VecGetArrayRead(tao->constraints_equality,&garr);
76: PetscMemcpy(carr,garr,pdipm->ng*sizeof(PetscScalar));
77: VecRestoreArrayRead(tao->constraints_equality,&garr);
78: }
80: /* (1.b) Update xfixed */
81: if (pdipm->Nxfixed) {
82: offset = pdipm->ng;
83: ISGetIndices(pdipm->isxfixed,&fxptr); /* global indices in x */
84: for (k=0;k < pdipm->nxfixed;k++){
85: i = fxptr[k]-xstart;
86: carr[offset + k] = xarr[i] - xuarr[i];
87: }
88: }
89: VecRestoreArray(pdipm->ce,&carr);
91: /* (2) Update ci vector */
92: VecGetArray(pdipm->ci,&carr);
94: if (pdipm->Nh) {
95: /* (2.a) Inserting updated h(x) */
96: VecGetArrayRead(tao->constraints_inequality,&harr);
97: PetscMemcpy(carr,harr,pdipm->nh*sizeof(PetscScalar));
98: VecRestoreArrayRead(tao->constraints_inequality,&harr);
99: }
101: /* (2.b) Update xub */
102: offset = pdipm->nh;
103: if (pdipm->Nxub) {
104: ISGetIndices(pdipm->isxub,&ubptr);
105: for (k=0; k<pdipm->nxub; k++){
106: i = ubptr[k]-xstart;
107: carr[offset + k] = xuarr[i] - xarr[i];
108: }
109: }
111: if (pdipm->Nxlb) {
112: /* (2.c) Update xlb */
113: offset += pdipm->nxub;
114: ISGetIndices(pdipm->isxlb,&lbptr); /* global indices in x */
115: for (k=0; k<pdipm->nxlb; k++){
116: i = lbptr[k]-xstart;
117: carr[offset + k] = xarr[i] - xlarr[i];
118: }
119: }
121: if (pdipm->Nxbox) {
122: /* (2.d) Update xbox */
123: offset += pdipm->nxlb;
124: offset1 = offset + pdipm->nxbox;
125: ISGetIndices(pdipm->isxbox,&bxptr); /* global indices in x */
126: for (k=0; k<pdipm->nxbox; k++){
127: i = bxptr[k]-xstart; /* local indices in x */
128: carr[offset+k] = xuarr[i] - xarr[i];
129: carr[offset1+k] = xarr[i] - xlarr[i];
130: }
131: }
132: VecRestoreArray(pdipm->ci,&carr);
134: /* Restoring Vectors */
135: VecRestoreArrayRead(x,&xarr);
136: VecRestoreArrayRead(tao->XU,&xuarr);
137: VecRestoreArrayRead(tao->XL,&xlarr);
138: return(0);
139: }
141: /*
142: TaoPDIPMSetUpBounds - Create upper and lower bound vectors of x
144: Collective
146: Input Parameter:
147: . tao - holds pdipm and XL & XU
149: Level: beginner
151: .seealso: TaoPDIPMUpdateConstraints
152: */
153: PetscErrorCode TaoPDIPMSetUpBounds(Tao tao)
154: {
155: PetscErrorCode ierr;
156: TAO_PDIPM *pdipm=(TAO_PDIPM*)tao->data;
157: const PetscScalar *xl,*xu;
158: PetscInt n,*ixlb,*ixub,*ixfixed,*ixfree,*ixbox,i,low,high,idx;
159: MPI_Comm comm;
160: PetscInt sendbuf[5],recvbuf[5];
163: /* Creates upper and lower bounds vectors on x, if not created already */
164: TaoComputeVariableBounds(tao);
166: VecGetLocalSize(tao->XL,&n);
167: PetscMalloc5(n,&ixlb,n,&ixub,n,&ixfree,n,&ixfixed,n,&ixbox);
169: VecGetOwnershipRange(tao->XL,&low,&high);
170: VecGetArrayRead(tao->XL,&xl);
171: VecGetArrayRead(tao->XU,&xu);
172: for (i=0; i<n; i++) {
173: idx = low + i;
174: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
175: if (PetscRealPart(xl[i]) == PetscRealPart(xu[i])) {
176: ixfixed[pdipm->nxfixed++] = idx;
177: } else ixbox[pdipm->nxbox++] = idx;
178: } else {
179: if ((PetscRealPart(xl[i]) > PETSC_NINFINITY) && (PetscRealPart(xu[i]) >= PETSC_INFINITY)) {
180: ixlb[pdipm->nxlb++] = idx;
181: } else if ((PetscRealPart(xl[i]) <= PETSC_NINFINITY) && (PetscRealPart(xu[i]) < PETSC_INFINITY)) {
182: ixub[pdipm->nxlb++] = idx;
183: } else ixfree[pdipm->nxfree++] = idx;
184: }
185: }
186: VecRestoreArrayRead(tao->XL,&xl);
187: VecRestoreArrayRead(tao->XU,&xu);
189: PetscObjectGetComm((PetscObject)tao,&comm);
190: sendbuf[0] = pdipm->nxlb;
191: sendbuf[1] = pdipm->nxub;
192: sendbuf[2] = pdipm->nxfixed;
193: sendbuf[3] = pdipm->nxbox;
194: sendbuf[4] = pdipm->nxfree;
196: MPI_Allreduce(sendbuf,recvbuf,5,MPIU_INT,MPI_SUM,comm);
197: pdipm->Nxlb = recvbuf[0];
198: pdipm->Nxub = recvbuf[1];
199: pdipm->Nxfixed = recvbuf[2];
200: pdipm->Nxbox = recvbuf[3];
201: pdipm->Nxfree = recvbuf[4];
203: if (pdipm->Nxlb) {
204: ISCreateGeneral(comm,pdipm->nxlb,ixlb,PETSC_COPY_VALUES,&pdipm->isxlb);
205: }
206: if (pdipm->Nxub) {
207: ISCreateGeneral(comm,pdipm->nxub,ixub,PETSC_COPY_VALUES,&pdipm->isxub);
208: }
209: if (pdipm->Nxfixed) {
210: ISCreateGeneral(comm,pdipm->nxfixed,ixfixed,PETSC_COPY_VALUES,&pdipm->isxfixed);
211: }
212: if (pdipm->Nxbox) {
213: ISCreateGeneral(comm,pdipm->nxbox,ixbox,PETSC_COPY_VALUES,&pdipm->isxbox);
214: }
215: if (pdipm->Nxfree) {
216: ISCreateGeneral(comm,pdipm->nxfree,ixfree,PETSC_COPY_VALUES,&pdipm->isxfree);
217: }
218: PetscFree5(ixlb,ixub,ixfixed,ixbox,ixfree);
219: return(0);
220: }
222: /*
223: TaoPDIPMInitializeSolution - Initialize PDIPM solution X = [x; lambdae; lambdai; z].
224: X consists of four subvectors in the order [x; lambdae; lambdai; z]. These
225: four subvectors need to be initialized and its values copied over to X. Instead
226: of copying, we use VecPlace/ResetArray functions to share the memory locations for
227: X and the subvectors
229: Collective
231: Input Parameter:
232: . tao - Tao context
234: Level: beginner
235: */
236: PetscErrorCode TaoPDIPMInitializeSolution(Tao tao)
237: {
239: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
240: PetscScalar *Xarr,*z,*lambdai;
241: PetscInt i;
242: const PetscScalar *xarr,*h;
245: VecGetArray(pdipm->X,&Xarr);
247: /* Set Initialize X.x = tao->solution */
248: VecGetArrayRead(tao->solution,&xarr);
249: PetscMemcpy(Xarr,xarr,pdipm->nx*sizeof(PetscScalar));
250: VecRestoreArrayRead(tao->solution,&xarr);
252: /* Initialize X.lambdae = 0.0 */
253: if (pdipm->lambdae) {
254: VecSet(pdipm->lambdae,0.0);
255: }
256: /* Initialize X.lambdai = push_init_lambdai, X.z = push_init_slack */
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: VecGetArray(pdipm->lambdai,&lambdai);
262: VecGetArray(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: VecRestoreArray(pdipm->lambdai,&lambdai);
272: VecRestoreArray(pdipm->z,&z);
274: VecRestoreArray(pdipm->X,&Xarr);
275: return(0);
276: }
278: /*
279: TaoSNESJacobian_PDIPM - Evaluate the Hessian matrix at X
281: Input Parameter:
282: snes - SNES context
283: X - KKT Vector
284: *ctx - pdipm context
286: Output Parameter:
287: J - Hessian matrix
288: Jpre - Preconditioner
289: */
290: PetscErrorCode TaoSNESJacobian_PDIPM(SNES snes,Vec X, Mat J, Mat Jpre, void *ctx)
291: {
292: PetscErrorCode ierr;
293: Tao tao=(Tao)ctx;
294: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
295: PetscInt i,row,cols[2],Jrstart,rjstart,nc,j;
296: const PetscInt *aj,*ranges,*Jranges,*rranges,*cranges;
297: const PetscScalar *Xarr,*aa;
298: PetscScalar vals[2];
299: PetscInt proc,nx_all,*nce_all=pdipm->nce_all;
300: MPI_Comm comm;
301: PetscMPIInt rank,size;
302: Mat jac_equality_trans=pdipm->jac_equality_trans,jac_inequality_trans=pdipm->jac_inequality_trans;
305: PetscObjectGetComm((PetscObject)snes,&comm);
306: MPI_Comm_rank(comm,&rank);
307: MPI_Comm_rank(comm,&size);
309: MatGetOwnershipRanges(Jpre,&Jranges);
310: MatGetOwnershipRange(Jpre,&Jrstart,NULL);
311: MatGetOwnershipRangesColumn(tao->hessian,&rranges);
312: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
314: VecGetArrayRead(X,&Xarr);
316: /* (2) insert Z and Ci to Jpre -- overwrite existing values */
317: for (i=0; i < pdipm->nci; i++) {
318: row = Jrstart + pdipm->off_z + i;
319: cols[0] = Jrstart + pdipm->off_lambdai + i;
320: cols[1] = row;
321: vals[0] = Xarr[pdipm->off_z + i];
322: vals[1] = Xarr[pdipm->off_lambdai + i];
323: MatSetValues(Jpre,1,&row,2,cols,vals,INSERT_VALUES);
324: }
326: /* (3) insert 2nd row block of Jpre: [ grad g, 0, 0, 0] */
327: if (pdipm->Ng) {
328: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
329: for (i=0; i<pdipm->ng; i++){
330: row = Jrstart + pdipm->off_lambdae + i;
331:
332: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
333: proc = 0;
334: for (j=0; j < nc; j++) {
335: while (aj[j] >= cranges[proc+1]) proc++;
336: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
337: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
338: }
339: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,&aa);
340: }
341: }
343: if (pdipm->Nh) {
344: /* (4) insert 3nd row block of Jpre: [ grad h, 0, 0, 0] */
345: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
346: for (i=0; i < pdipm->nh; i++){
347: row = Jrstart + pdipm->off_lambdai + i;
348:
349: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
350: proc = 0;
351: for (j=0; j < nc; j++) {
352: while (aj[j] >= cranges[proc+1]) proc++;
353: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
354: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
355: }
356: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,&aa);
357: }
358: }
360: /* (5) insert Wxx, grad g' and -grad h' to Jpre */
361: if (pdipm->Ng) {
362: MatTranspose(tao->jacobian_equality,MAT_REUSE_MATRIX,&jac_equality_trans);
363: }
364: if (pdipm->Nh) {
365: MatTranspose(tao->jacobian_inequality,MAT_REUSE_MATRIX,&jac_inequality_trans);
366: }
368: VecPlaceArray(pdipm->x,Xarr);
369: TaoComputeHessian(tao,pdipm->x,tao->hessian,tao->hessian_pre);
370: VecResetArray(pdipm->x);
372: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
373: for (i=0; i<pdipm->nx; i++){
374: row = Jrstart + i;
376: /* insert Wxx */
377: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
378: proc = 0;
379: for (j=0; j < nc; j++) {
380: while (aj[j] >= cranges[proc+1]) proc++;
381: cols[0] = aj[j] - cranges[proc] + Jranges[proc];
382: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
383: }
384: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,&aa);
386: if (pdipm->ng) {
387: /* insert grad g' */
388: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
389: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
390: proc = 0;
391: for (j=0; j < nc; j++) {
392: /* find row ownership of */
393: while (aj[j] >= ranges[proc+1]) proc++;
394: nx_all = rranges[proc+1] - rranges[proc];
395: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
396: MatSetValue(Jpre,row,cols[0],aa[j],INSERT_VALUES);
397: }
398: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,&aa);
399: }
401: if (pdipm->nh) {
402: /* insert -grad h' */
403: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
404: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
405: proc = 0;
406: for (j=0; j < nc; j++) {
407: /* find row ownership of */
408: while (aj[j] >= ranges[proc+1]) proc++;
409: nx_all = rranges[proc+1] - rranges[proc];
410: cols[0] = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
411: MatSetValue(Jpre,row,cols[0],-aa[j],INSERT_VALUES);
412: }
413: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,&aa);
414: }
415: }
416: VecRestoreArrayRead(X,&Xarr);
418: /* (6) assemble Jpre and J */
419: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
420: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
422: if (J != Jpre) {
423: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
424: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
425: }
426: return(0);
427: }
429: /*
430: TaoSnesFunction_PDIPM - Evaluate KKT function at X
432: Input Parameter:
433: snes - SNES context
434: X - KKT Vector
435: *ctx - pdipm
437: Output Parameter:
438: F - Updated Lagrangian vector
439: */
440: PetscErrorCode TaoSNESFunction_PDIPM(SNES snes,Vec X,Vec F,void *ctx)
441: {
443: Tao tao=(Tao)ctx;
444: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
445: PetscScalar *Farr;
446: Vec x,L1;
447: PetscInt i;
448: PetscReal res[2],cnorm[2];
449: const PetscScalar *Xarr,*carr,*zarr,*larr;
452: VecSet(F,0.0);
454: VecGetArrayRead(X,&Xarr);
455: VecGetArray(F,&Farr);
457: /* (0) Evaluate f, fx, Gx, Hx at X.x Note: pdipm->x is not changed below */
458: x = pdipm->x;
459: VecPlaceArray(x,Xarr);
460: TaoPDIPMEvaluateFunctionsAndJacobians(tao,x);
462: /* Update ce, ci, and Jci at X.x */
463: TaoPDIPMUpdateConstraints(tao,x);
464: VecResetArray(x);
466: /* (1) L1 = fx + (gradG'*DE + Jce_xfixed'*lambdae_xfixed) - (gradH'*DI + Jci_xb'*lambdai_xb) */
467: L1 = pdipm->x;
468: VecPlaceArray(L1,Farr);
469: if (pdipm->Nci) {
470: if (pdipm->Nh) {
471: /* L1 += gradH'*DI. Note: tao->DI is not changed below */
472: VecPlaceArray(tao->DI,Xarr+pdipm->off_lambdai);
473: MatMultTransposeAdd(tao->jacobian_inequality,tao->DI,L1,L1);
474: VecResetArray(tao->DI);
475: }
477: /* L1 += Jci_xb'*lambdai_xb */
478: VecPlaceArray(pdipm->lambdai_xb,Xarr+pdipm->off_lambdai+pdipm->nh);
479: MatMultTransposeAdd(pdipm->Jci_xb,pdipm->lambdai_xb,L1,L1);
480: VecResetArray(pdipm->lambdai_xb);
482: /* (1.4) L1 = - (gradH'*DI + Jci_xb'*lambdai_xb) */
483: VecScale(L1,-1.0);
484: }
486: /* L1 += fx */
487: VecAXPY(L1,1.0,tao->gradient);
489: if (pdipm->Nce) {
490: if (pdipm->Ng) {
491: /* L1 += gradG'*DE. Note: tao->DE is not changed below */
492: VecPlaceArray(tao->DE,Xarr+pdipm->off_lambdae);
493: MatMultTransposeAdd(tao->jacobian_equality,tao->DE,L1,L1);
494: VecResetArray(tao->DE);
495: }
496: if (pdipm->Nxfixed) {
497: /* L1 += Jce_xfixed'*lambdae_xfixed */
498: VecPlaceArray(pdipm->lambdae_xfixed,Xarr+pdipm->off_lambdae+pdipm->ng);
499: MatMultTransposeAdd(pdipm->Jce_xfixed,pdipm->lambdae_xfixed,L1,L1);
500: VecResetArray(pdipm->lambdae_xfixed);
501: }
502: }
503: VecNorm(L1,NORM_2,&res[0]);
504: VecResetArray(L1);
506: /* (2) L2 = ce(x) */
507: if (pdipm->Nce) {
508: VecGetArrayRead(pdipm->ce,&carr);
509: for (i=0; i<pdipm->nce; i++) Farr[pdipm->off_lambdae + i] = carr[i];
510: VecRestoreArrayRead(pdipm->ce,&carr);
511: }
512: VecNorm(pdipm->ce,NORM_2,&cnorm[0]);
514: if (pdipm->Nci) {
515: /* (3) L3 = ci(x) - z;
516: (4) L4 = Z * Lambdai * e - mu * e
517: */
518: VecGetArrayRead(pdipm->ci,&carr);
519: larr = Xarr+pdipm->off_lambdai;
520: zarr = Xarr+pdipm->off_z;
521: for (i=0; i<pdipm->nci; i++) {
522: Farr[pdipm->off_lambdai + i] = carr[i] - zarr[i];
523: Farr[pdipm->off_z + i] = zarr[i]*larr[i] - pdipm->mu;
524: }
525: VecRestoreArrayRead(pdipm->ci,&carr);
526: }
528: VecPlaceArray(pdipm->ci,Farr+pdipm->off_lambdai);
529: VecNorm(pdipm->ci,NORM_2,&cnorm[1]);
530: VecResetArray(pdipm->ci);
532: /* note: pdipm->z is not changed below */
533: VecPlaceArray(pdipm->z,Farr+pdipm->off_z);
534: VecNorm(pdipm->z,NORM_2,&res[1]);
535: VecResetArray(pdipm->z);
537: tao->residual = PetscSqrtReal(res[0]*res[0] + res[1]*res[1]);
538: tao->cnorm = PetscSqrtReal(cnorm[0]*cnorm[0] + cnorm[1]*cnorm[1]);
540: VecRestoreArrayRead(X,&Xarr);
541: VecRestoreArray(F,&Farr);
542: return(0);
543: }
545: /*
546: PDIPMLineSearch - Custom line search used with PDIPM.
548: Collective on TAO
550: Notes:
551: PDIPMLineSearch employs a simple backtracking line-search to keep
552: the slack variables (z) and inequality constraints lagrange multipliers
553: (lambdai) positive, i.e., z,lambdai >=0. It does this by calculating scalars
554: alpha_p and alpha_d to keep z,lambdai non-negative. The decision (x), and the
555: slack variables are updated as X = X + alpha_d*dx. The constraint multipliers
556: are updated as Lambdai = Lambdai + alpha_p*dLambdai. The barrier parameter mu
557: is also updated as mu = mu + z'lambdai/Nci
558: */
559: PetscErrorCode PDIPMLineSearch(SNESLineSearch linesearch,void *ctx)
560: {
562: Tao tao=(Tao)ctx;
563: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
564: SNES snes;
565: Vec X,F,Y,W,G;
566: PetscInt i,iter;
567: PetscReal alpha_p=1.0,alpha_d=1.0,alpha[4];
568: PetscScalar *Xarr,*z,*lambdai,dot;
569: const PetscScalar *dXarr,*dz,*dlambdai;
570: PetscScalar *taosolarr;
573: SNESLineSearchGetSNES(linesearch,&snes);
574: SNESGetIterationNumber(snes,&iter);
576: SNESLineSearchSetReason(linesearch,SNES_LINESEARCH_SUCCEEDED);
577: SNESLineSearchGetVecs(linesearch,&X,&F,&Y,&W,&G);
579: VecGetArray(X,&Xarr);
580: VecGetArrayRead(Y,&dXarr);
581: z = Xarr + pdipm->off_z;
582: dz = dXarr + pdipm->off_z;
583: for (i=0; i < pdipm->nci; i++) {
584: if (z[i] - dz[i] < 0.0) {
585: alpha_p = PetscMin(alpha_p,0.9999*z[i]/dz[i]);
586: }
587: }
589: lambdai = Xarr + pdipm->off_lambdai;
590: dlambdai = dXarr + pdipm->off_lambdai;
592: for (i=0; i<pdipm->nci; i++) {
593: if (lambdai[i] - dlambdai[i] < 0.0) {
594: alpha_d = PetscMin(0.9999*lambdai[i]/dlambdai[i],alpha_d);
595: }
596: }
598: alpha[0] = alpha_p;
599: alpha[1] = alpha_d;
600: VecRestoreArrayRead(Y,&dXarr);
601: VecRestoreArray(X,&Xarr);
603: /* alpha = min(alpha) over all processes */
604: MPI_Allreduce(alpha,alpha+2,2,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)tao));
606: alpha_p = alpha[2];
607: alpha_d = alpha[3];
609: VecGetArray(X,&Xarr);
610: VecGetArrayRead(Y,&dXarr);
611: for (i=0; i<pdipm->nx; i++) {
612: Xarr[i] = Xarr[i] - alpha_p * dXarr[i];
613: }
615: for (i=0; i<pdipm->nce; i++) {
616: Xarr[i+pdipm->off_lambdae] = Xarr[i+pdipm->off_lambdae] - alpha_d * dXarr[i+pdipm->off_lambdae];
617: }
619: for (i=0; i<pdipm->nci; i++) {
620: Xarr[i+pdipm->off_lambdai] = Xarr[i+pdipm->off_lambdai] - alpha_d * dXarr[i+pdipm->off_lambdai];
621: }
623: for (i=0; i<pdipm->nci; i++) {
624: Xarr[i+pdipm->off_z] = Xarr[i+pdipm->off_z] - alpha_p * dXarr[i+pdipm->off_z];
625: }
627: VecGetArray(tao->solution,&taosolarr);
628: PetscMemcpy(taosolarr,Xarr,pdipm->nx*sizeof(PetscScalar));
629: VecRestoreArray(tao->solution,&taosolarr);
632: VecRestoreArray(X,&Xarr);
633: VecRestoreArrayRead(Y,&dXarr);
635: /* Evaluate F at X */
636: SNESComputeFunction(snes,X,F);
637: SNESLineSearchComputeNorms(linesearch); /* must call this func, do not know why */
639: /* update mu = mu_update_factor * dot(z,lambdai)/pdipm->nci at updated X */
640: VecDot(pdipm->z,pdipm->lambdai,&dot);
642: /* if (PetscAbsReal(pdipm->gradL) < 0.9*pdipm->mu) */
643: pdipm->mu = pdipm->mu_update_factor * dot/pdipm->Nci;
645: /* Update F; get tao->residual and tao->cnorm */
646: TaoSNESFunction_PDIPM(snes,X,F,(void*)tao);
648: tao->niter++;
649: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
650: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
652: (*tao->ops->convergencetest)(tao,tao->cnvP);
653: if (tao->reason) {
654: SNESSetConvergedReason(snes,SNES_CONVERGED_FNORM_ABS);
655: }
656: return(0);
657: }
659: /*
660: TaoSolve_PDIPM
662: Input Parameter:
663: tao - TAO context
665: Output Parameter:
666: tao - TAO context
667: */
668: PetscErrorCode TaoSolve_PDIPM(Tao tao)
669: {
670: PetscErrorCode ierr;
671: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
672: SNESLineSearch linesearch; /* SNESLineSearch context */
673: Vec dummy;
676: if (!tao->constraints_equality && !tao->constraints_inequality) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_NULL,"Equality and inequality contraints are not set. Either set them or switch to a different algorithm");
678: /* Initialize all variables */
679: TaoPDIPMInitializeSolution(tao);
681: /* Set linesearch */
682: SNESGetLineSearch(pdipm->snes,&linesearch);
683: SNESLineSearchSetType(linesearch,SNESLINESEARCHSHELL);
684: SNESLineSearchShellSetUserFunc(linesearch,PDIPMLineSearch,tao);
685: SNESLineSearchSetFromOptions(linesearch);
687: tao->reason = TAO_CONTINUE_ITERATING;
689: /* -tao_monitor for iteration 0 and check convergence */
690: VecDuplicate(pdipm->X,&dummy);
691: TaoSNESFunction_PDIPM(pdipm->snes,pdipm->X,dummy,(void*)tao);
693: TaoLogConvergenceHistory(tao,pdipm->obj,tao->residual,tao->cnorm,tao->niter);
694: TaoMonitor(tao,tao->niter,pdipm->obj,tao->residual,tao->cnorm,pdipm->mu);
695: VecDestroy(&dummy);
696: (*tao->ops->convergencetest)(tao,tao->cnvP);
697: if (tao->reason) {
698: SNESSetConvergedReason(pdipm->snes,SNES_CONVERGED_FNORM_ABS);
699: }
701: while (tao->reason == TAO_CONTINUE_ITERATING) {
702: SNESConvergedReason reason;
703: SNESSolve(pdipm->snes,NULL,pdipm->X);
705: /* Check SNES convergence */
706: SNESGetConvergedReason(pdipm->snes,&reason);
707: if (reason < 0) {
708: PetscPrintf(PETSC_COMM_WORLD,"SNES solve did not converged due to reason %D\n",reason);
709: }
711: /* Check TAO convergence */
712: if (PetscIsInfOrNanReal(pdipm->obj)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"User-provided compute function generated Inf or NaN");
713: }
714: return(0);
715: }
717: /*
718: TaoSetup_PDIPM - Sets up tao and pdipm
720: Input Parameter:
721: tao - TAO object
723: Output: pdipm - initialized object
724: */
725: PetscErrorCode TaoSetup_PDIPM(Tao tao)
726: {
727: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
729: MPI_Comm comm;
730: PetscMPIInt rank,size;
731: PetscInt row,col,Jcrstart,Jcrend,k,tmp,nc,proc,*nh_all,*ng_all;
732: PetscInt offset,*xa,*xb,i,j,rstart,rend;
733: PetscScalar one=1.0,neg_one=-1.0,*Xarr;
734: const PetscInt *cols,*rranges,*cranges,*aj,*ranges;
735: const PetscScalar *aa;
736: Mat J,jac_equality_trans,jac_inequality_trans;
737: Mat Jce_xfixed_trans,Jci_xb_trans;
738: PetscInt *dnz,*onz,rjstart,nx_all,*nce_all,*Jranges,cols1[2];
741: PetscObjectGetComm((PetscObject)tao,&comm);
742: MPI_Comm_rank(comm,&rank);
743: MPI_Comm_size(comm,&size);
745: /* (1) Setup Bounds and create Tao vectors */
746: TaoPDIPMSetUpBounds(tao);
748: if (!tao->gradient) {
749: VecDuplicate(tao->solution,&tao->gradient);
750: VecDuplicate(tao->solution,&tao->stepdirection);
751: }
753: /* (2) Get sizes */
754: /* Size of vector x - This is set by TaoSetInitialVector */
755: VecGetSize(tao->solution,&pdipm->Nx);
756: VecGetLocalSize(tao->solution,&pdipm->nx);
758: /* Size of equality constraints and vectors */
759: if (tao->constraints_equality) {
760: VecGetSize(tao->constraints_equality,&pdipm->Ng);
761: VecGetLocalSize(tao->constraints_equality,&pdipm->ng);
762: } else {
763: pdipm->ng = pdipm->Ng = 0;
764: }
766: pdipm->nce = pdipm->ng + pdipm->nxfixed;
767: pdipm->Nce = pdipm->Ng + pdipm->Nxfixed;
769: /* Size of inequality constraints and vectors */
770: if (tao->constraints_inequality) {
771: VecGetSize(tao->constraints_inequality,&pdipm->Nh);
772: VecGetLocalSize(tao->constraints_inequality,&pdipm->nh);
773: } else {
774: pdipm->nh = pdipm->Nh = 0;
775: }
777: pdipm->nci = pdipm->nh + pdipm->nxlb + pdipm->nxub + 2*pdipm->nxbox;
778: pdipm->Nci = pdipm->Nh + pdipm->Nxlb + pdipm->Nxub + 2*pdipm->Nxbox;
780: /* Full size of the KKT system to be solved */
781: pdipm->n = pdipm->nx + pdipm->nce + 2*pdipm->nci;
782: pdipm->N = pdipm->Nx + pdipm->Nce + 2*pdipm->Nci;
784: /* list below to TaoView_PDIPM()? */
785: /* PetscPrintf(PETSC_COMM_SELF,"[%d] nce %d = ng %d + nxfixed %d\n",rank,pdipm->nce,pdipm->ng,pdipm->nxfixed); */
786: /* PetscPrintf(PETSC_COMM_SELF,"[%d] nci %d = nh %d + nxlb %d + nxub %d + 2*nxbox %d\n",rank,pdipm->nci,pdipm->nh,pdipm->nxlb,pdipm->nxub,pdipm->nxbox); */
787: /* PetscPrintf(PETSC_COMM_SELF,"[%d] n %d = nx %d + nce %d + 2*nci %d\n",rank,pdipm->n,pdipm->nx,pdipm->nce,pdipm->nci); */
789: /* (3) Offsets for subvectors */
790: pdipm->off_lambdae = pdipm->nx;
791: pdipm->off_lambdai = pdipm->off_lambdae + pdipm->nce;
792: pdipm->off_z = pdipm->off_lambdai + pdipm->nci;
794: /* (4) Create vectors and subvectors */
795: /* Ce and Ci vectors */
796: VecCreate(comm,&pdipm->ce);
797: VecSetSizes(pdipm->ce,pdipm->nce,pdipm->Nce);
798: VecSetFromOptions(pdipm->ce);
800: VecCreate(comm,&pdipm->ci);
801: VecSetSizes(pdipm->ci,pdipm->nci,pdipm->Nci);
802: VecSetFromOptions(pdipm->ci);
804: /* X=[x; lambdae; lambdai; z] for the big KKT system */
805: VecCreate(comm,&pdipm->X);
806: VecSetSizes(pdipm->X,pdipm->n,pdipm->N);
807: VecSetFromOptions(pdipm->X);
809: /* Subvectors; they share local arrays with X */
810: VecGetArray(pdipm->X,&Xarr);
811: /* x shares local array with X.x */
812: if (pdipm->Nx) {
813: VecCreateMPIWithArray(comm,1,pdipm->nx,pdipm->Nx,Xarr,&pdipm->x);
814: }
816: /* lambdae shares local array with X.lambdae */
817: if (pdipm->Nce) {
818: VecCreateMPIWithArray(comm,1,pdipm->nce,pdipm->Nce,Xarr+pdipm->off_lambdae,&pdipm->lambdae);
819: }
821: /* tao->DE shares local array with X.lambdae_g */
822: if (pdipm->Ng) {
823: VecCreateMPIWithArray(comm,1,pdipm->ng,pdipm->Ng,Xarr+pdipm->off_lambdae,&tao->DE);
825: VecCreate(comm,&pdipm->lambdae_xfixed);
826: VecSetSizes(pdipm->lambdae_xfixed,pdipm->nxfixed,PETSC_DECIDE);
827: VecSetFromOptions(pdipm->lambdae_xfixed);
828: }
830: if (pdipm->Nci) {
831: /* lambdai shares local array with X.lambdai */
832: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_lambdai,&pdipm->lambdai);
834: /* z for slack variables; it shares local array with X.z */
835: VecCreateMPIWithArray(comm,1,pdipm->nci,pdipm->Nci,Xarr+pdipm->off_z,&pdipm->z);
836: }
838: /* tao->DI which shares local array with X.lambdai_h */
839: if (pdipm->Nh) {
840: VecCreateMPIWithArray(comm,1,pdipm->nh,pdipm->Nh,Xarr+pdipm->off_lambdai,&tao->DI);
841: }
843: VecCreate(comm,&pdipm->lambdai_xb);
844: VecSetSizes(pdipm->lambdai_xb,(pdipm->nci - pdipm->nh),PETSC_DECIDE);
845: VecSetFromOptions(pdipm->lambdai_xb);
847: VecRestoreArray(pdipm->X,&Xarr);
849: /* (5) Create Jacobians Jce_xfixed and Jci */
850: /* (5.1) PDIPM Jacobian of equality bounds cebound(x) = J_nxfixed */
851: if (pdipm->Nxfixed) {
852: /* Create Jce_xfixed */
853: MatCreate(comm,&pdipm->Jce_xfixed);
854: MatSetSizes(pdipm->Jce_xfixed,pdipm->nxfixed,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
855: MatSetFromOptions(pdipm->Jce_xfixed);
856: MatSeqAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL);
857: MatMPIAIJSetPreallocation(pdipm->Jce_xfixed,1,NULL,1,NULL);
859: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,&Jcrend);
860: ISGetIndices(pdipm->isxfixed,&cols);
861: k = 0;
862: for (row = Jcrstart; row < Jcrend; row++) {
863: MatSetValues(pdipm->Jce_xfixed,1,&row,1,cols+k,&one,INSERT_VALUES);
864: k++;
865: }
866: ISRestoreIndices(pdipm->isxfixed, &cols);
867: MatAssemblyBegin(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
868: MatAssemblyEnd(pdipm->Jce_xfixed,MAT_FINAL_ASSEMBLY);
869: }
871: /* (5.2) PDIPM inequality Jacobian Jci = [tao->jacobian_inequality; ...] */
872: MatCreate(comm,&pdipm->Jci_xb);
873: MatSetSizes(pdipm->Jci_xb,pdipm->nci-pdipm->nh,pdipm->nx,PETSC_DECIDE,pdipm->Nx);
874: MatSetFromOptions(pdipm->Jci_xb);
875: MatSeqAIJSetPreallocation(pdipm->Jci_xb,1,NULL);
876: MatMPIAIJSetPreallocation(pdipm->Jci_xb,1,NULL,1,NULL);
878: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,&Jcrend);
879: offset = Jcrstart;
880: if (pdipm->Nxub) {
881: /* Add xub to Jci_xb */
882: ISGetIndices(pdipm->isxub,&cols);
883: k = 0;
884: for (row = offset; row < offset + pdipm->nxub; row++) {
885: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
886: k++;
887: }
888: ISRestoreIndices(pdipm->isxub, &cols);
889: }
891: if (pdipm->Nxlb) {
892: /* Add xlb to Jci_xb */
893: ISGetIndices(pdipm->isxlb,&cols);
894: k = 0;
895: offset += pdipm->nxub;
896: for (row = offset; row < offset + pdipm->nxlb; row++) {
897: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&one,INSERT_VALUES);
898: k++;
899: }
900: ISRestoreIndices(pdipm->isxlb, &cols);
901: }
903: /* Add xbox to Jci_xb */
904: if (pdipm->Nxbox) {
905: ISGetIndices(pdipm->isxbox,&cols);
906: k = 0;
907: offset += pdipm->nxlb;
908: for (row = offset; row < offset + pdipm->nxbox; row++) {
909: MatSetValues(pdipm->Jci_xb,1,&row,1,cols+k,&neg_one,INSERT_VALUES);
910: tmp = row + pdipm->nxbox;
911: MatSetValues(pdipm->Jci_xb,1,&tmp,1,cols+k,&one,INSERT_VALUES);
912: k++;
913: }
914: ISRestoreIndices(pdipm->isxbox, &cols);
915: }
917: MatAssemblyBegin(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
918: MatAssemblyEnd(pdipm->Jci_xb,MAT_FINAL_ASSEMBLY);
919: /* MatView(pdipm->Jci_xb,PETSC_VIEWER_STDOUT_WORLD); */
921: /* (6) Set up ISs for PC Fieldsplit */
922: if (pdipm->solve_reduced_kkt) {
923: PetscMalloc2(pdipm->nx+pdipm->nce,&xa,2*pdipm->nci,&xb);
924: for(i=0; i < pdipm->nx + pdipm->nce; i++) xa[i] = i;
925: for(i=0; i < 2*pdipm->nci; i++) xb[i] = pdipm->off_lambdai + i;
927: ISCreateGeneral(comm,pdipm->nx+pdipm->nce,xa,PETSC_OWN_POINTER,&pdipm->is1);
928: ISCreateGeneral(comm,2*pdipm->nci,xb,PETSC_OWN_POINTER,&pdipm->is2);
929: }
931: /* (7) Gather offsets from all processes */
932: PetscMalloc1(size,&pdipm->nce_all);
934: /* Get rstart of KKT matrix */
935: MPI_Scan(&pdipm->n,&rstart,1,MPIU_INT,MPI_SUM,comm);
936: rstart -= pdipm->n;
938: MPI_Allgather(&pdipm->nce,1,MPIU_INT,pdipm->nce_all,1,MPIU_INT,comm);
940: PetscMalloc3(size,&ng_all,size,&nh_all,size,&Jranges);
941: MPI_Allgather(&rstart,1,MPIU_INT,Jranges,1,MPIU_INT,comm);
942: MPI_Allgather(&pdipm->nh,1,MPIU_INT,nh_all,1,MPIU_INT,comm);
943: MPI_Allgather(&pdipm->ng,1,MPIU_INT,ng_all,1,MPIU_INT,comm);
945: MatGetOwnershipRanges(tao->hessian,&rranges);
946: MatGetOwnershipRangesColumn(tao->hessian,&cranges);
948: if (pdipm->Ng) {
949: TaoComputeJacobianEquality(tao,tao->solution,tao->jacobian_equality,tao->jacobian_equality_pre);
950: MatTranspose(tao->jacobian_equality,MAT_INITIAL_MATRIX,&pdipm->jac_equality_trans);
951: }
952: if (pdipm->Nh) {
953: TaoComputeJacobianInequality(tao,tao->solution,tao->jacobian_inequality,tao->jacobian_inequality_pre);
954: MatTranspose(tao->jacobian_inequality,MAT_INITIAL_MATRIX,&pdipm->jac_inequality_trans);
955: }
957: /* Count dnz,onz for preallocation of KKT matrix */
958: jac_equality_trans = pdipm->jac_equality_trans;
959: jac_inequality_trans = pdipm->jac_inequality_trans;
960: nce_all = pdipm->nce_all;
962: if (pdipm->Nxfixed) {
963: MatTranspose(pdipm->Jce_xfixed,MAT_INITIAL_MATRIX,&Jce_xfixed_trans);
964: }
965: MatTranspose(pdipm->Jci_xb,MAT_INITIAL_MATRIX,&Jci_xb_trans);
967: MatPreallocateInitialize(comm,pdipm->n,pdipm->n,dnz,onz);
969: /* 1st row block of KKT matrix: [Wxx; gradCe'; -gradCi'; 0] */
970: TaoPDIPMEvaluateFunctionsAndJacobians(tao,pdipm->x);
971: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
973: /* Insert tao->hessian */
974: MatGetOwnershipRange(tao->hessian,&rjstart,NULL);
975: for (i=0; i<pdipm->nx; i++){
976: row = rstart + i;
978: MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
979: proc = 0;
980: for (j=0; j < nc; j++) {
981: while (aj[j] >= cranges[proc+1]) proc++;
982: col = aj[j] - cranges[proc] + Jranges[proc];
983: MatPreallocateSet(row,1,&col,dnz,onz);
984: }
985: MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);
987: if (pdipm->ng) {
988: /* Insert grad g' */
989: MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
990: MatGetOwnershipRanges(tao->jacobian_equality,&ranges);
991: proc = 0;
992: for (j=0; j < nc; j++) {
993: /* find row ownership of */
994: while (aj[j] >= ranges[proc+1]) proc++;
995: nx_all = rranges[proc+1] - rranges[proc];
996: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;
997: MatPreallocateSet(row,1,&col,dnz,onz);
998: }
999: MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);
1000: }
1002: /* Insert Jce_xfixed^T' */
1003: if (pdipm->nxfixed) {
1004: MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1005: MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);
1006: proc = 0;
1007: for (j=0; j < nc; j++) {
1008: /* find row ownership of */
1009: while (aj[j] >= ranges[proc+1]) proc++;
1010: nx_all = rranges[proc+1] - rranges[proc];
1011: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];
1012: MatPreallocateSet(row,1,&col,dnz,onz);
1013: }
1014: MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);
1015: }
1017: if (pdipm->nh) {
1018: /* Insert -grad h' */
1019: MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1020: MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);
1021: proc = 0;
1022: for (j=0; j < nc; j++) {
1023: /* find row ownership of */
1024: while (aj[j] >= ranges[proc+1]) proc++;
1025: nx_all = rranges[proc+1] - rranges[proc];
1026: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];
1027: MatPreallocateSet(row,1,&col,dnz,onz);
1028: }
1029: MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);
1030: }
1032: /* Insert Jci_xb^T' */
1033: MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1034: MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);
1035: proc = 0;
1036: for (j=0; j < nc; j++) {
1037: /* find row ownership of */
1038: while (aj[j] >= ranges[proc+1]) proc++;
1039: nx_all = rranges[proc+1] - rranges[proc];
1040: col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];
1041: MatPreallocateSet(row,1,&col,dnz,onz);
1042: }
1043: MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);
1044: }
1046: /* 2nd Row block of KKT matrix: [grad Ce, 0, 0, 0] */
1047: if (pdipm->Ng) {
1048: MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);
1049: for (i=0; i < pdipm->ng; i++){
1050: row = rstart + pdipm->off_lambdae + i;
1052: MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1053: proc = 0;
1054: for (j=0; j < nc; j++) {
1055: while (aj[j] >= cranges[proc+1]) proc++;
1056: col = aj[j] - cranges[proc] + Jranges[proc];
1057: MatPreallocateSet(row,1,&col,dnz,onz); /* grad g */
1058: }
1059: MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);
1060: }
1061: }
1062: /* Jce_xfixed */
1063: if (pdipm->Nxfixed) {
1064: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1065: for (i=0; i < (pdipm->nce - pdipm->ng); i++ ){
1066: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1068: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1069: if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1071: proc = 0;
1072: j = 0;
1073: while (cols[j] >= cranges[proc+1]) proc++;
1074: col = cols[j] - cranges[proc] + Jranges[proc];
1075: MatPreallocateSet(row,1,&col,dnz,onz);
1076: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);
1077: }
1078: }
1080: /* 3rd Row block of KKT matrix: [ gradCi, 0, 0, -I] */
1081: if (pdipm->Nh) {
1082: MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);
1083: for (i=0; i < pdipm->nh; i++){
1084: row = rstart + pdipm->off_lambdai + i;
1086: MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1087: proc = 0;
1088: for (j=0; j < nc; j++) {
1089: while (aj[j] >= cranges[proc+1]) proc++;
1090: col = aj[j] - cranges[proc] + Jranges[proc];
1091: MatPreallocateSet(row,1,&col,dnz,onz); /* grad h */
1092: }
1093: MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);
1094: }
1095: /* -I */
1096: for (i=0; i < pdipm->nh; i++){
1097: row = rstart + pdipm->off_lambdai + i;
1098: col = rstart + pdipm->off_z + i;
1099: MatPreallocateSet(row,1,&col,dnz,onz);
1100: }
1101: }
1103: /* Jci_xb */
1104: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1105: for (i=0; i < (pdipm->nci - pdipm->nh); i++ ){
1106: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1108: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1109: if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");
1110: proc = 0;
1111: for (j=0; j < nc; j++) {
1112: while (cols[j] >= cranges[proc+1]) proc++;
1113: col = cols[j] - cranges[proc] + Jranges[proc];
1114: MatPreallocateSet(row,1,&col,dnz,onz);
1115: }
1116: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);
1117: /* -I */
1118: col = rstart + pdipm->off_z + pdipm->nh + i;
1119: MatPreallocateSet(row,1,&col,dnz,onz);
1120: }
1122: /* 4-th Row block of KKT matrix: Z and Ci */
1123: for (i=0; i < pdipm->nci; i++) {
1124: row = rstart + pdipm->off_z + i;
1125: cols1[0] = rstart + pdipm->off_lambdai + i;
1126: cols1[1] = row;
1127: MatPreallocateSet(row,2,cols1,dnz,onz);
1128: }
1130: /* diagonal entry */
1131: for (i=0; i<pdipm->n; i++) dnz[i]++; /* diagonal entry */
1133: /* Create KKT matrix */
1134: MatCreate(comm,&J);
1135: MatSetSizes(J,pdipm->n,pdipm->n,PETSC_DECIDE,PETSC_DECIDE);
1136: MatSetFromOptions(J);
1137: MatSeqAIJSetPreallocation(J,0,dnz);
1138: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1139: /* MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); */
1140: MatPreallocateFinalize(dnz,onz);
1141: pdipm->K = J;
1143: /* (8) Set up nonlinear solver SNES */
1144: SNESSetFunction(pdipm->snes,NULL,TaoSNESFunction_PDIPM,(void*)tao);
1145: SNESSetJacobian(pdipm->snes,J,J,TaoSNESJacobian_PDIPM,(void*)tao);
1147: if (pdipm->solve_reduced_kkt) {
1148: PC pc;
1149: KSPGetPC(tao->ksp,&pc);
1150: PCSetType(pc,PCFIELDSPLIT);
1151: PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);
1152: PCFieldSplitSetIS(pc,"2",pdipm->is2);
1153: PCFieldSplitSetIS(pc,"1",pdipm->is1);
1154: }
1155: SNESSetFromOptions(pdipm->snes);
1157: /* (9) Insert constant entries to K */
1158: /* Set 0.0 to diagonal of K, so that the solver does not complain *about missing diagonal value */
1159: MatGetOwnershipRange(J,&rstart,&rend);
1160: for (i=rstart; i<rend; i++){
1161: MatSetValue(J,i,i,0.0,INSERT_VALUES);
1162: }
1164: /* Row block of K: [ grad Ce, 0, 0, 0] */
1165: if (pdipm->Nxfixed) {
1166: MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);
1167: for (i=0; i < (pdipm->nce - pdipm->ng); i++ ){
1168: row = rstart + pdipm->off_lambdae + pdipm->ng + i;
1170: MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1171: proc = 0;
1172: for (j=0; j < nc; j++) {
1173: while (cols[j] >= cranges[proc+1]) proc++;
1174: col = cols[j] - cranges[proc] + Jranges[proc];
1175: MatSetValue(J,row,col,aa[j],INSERT_VALUES); /* grad Ce */
1176: MatSetValue(J,col,row,aa[j],INSERT_VALUES); /* grad Ce' */
1177: }
1178: MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,&aa);
1179: }
1180: }
1182: /* Row block of K: [ grad Ci, 0, 0, -I] */
1183: MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);
1184: for (i=0; i < (pdipm->nci - pdipm->nh); i++ ){
1185: row = rstart + pdipm->off_lambdai + pdipm->nh + i;
1187: MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1188: proc = 0;
1189: for (j=0; j < nc; j++) {
1190: while (cols[j] >= cranges[proc+1]) proc++;
1191: col = cols[j] - cranges[proc] + Jranges[proc];
1192: MatSetValue(J,col,row,-aa[j],INSERT_VALUES);
1193: MatSetValue(J,row,col,aa[j],INSERT_VALUES);
1194: }
1195: MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,&aa);
1197: col = rstart + pdipm->off_z + pdipm->nh + i;
1198: MatSetValue(J,row,col,-1,INSERT_VALUES);
1199: }
1201: for (i=0; i < pdipm->nh; i++){
1202: row = rstart + pdipm->off_lambdai + i;
1203: col = rstart + pdipm->off_z + i;
1204: MatSetValue(J,row,col,-1,INSERT_VALUES);
1205: }
1207: if (pdipm->Nxfixed) {
1208: MatDestroy(&Jce_xfixed_trans);
1209: }
1210: MatDestroy(&Jci_xb_trans);
1211: PetscFree3(ng_all,nh_all,Jranges);
1212: return(0);
1213: }
1215: /*
1216: TaoDestroy_PDIPM - Destroys the pdipm object
1218: Input:
1219: full pdipm
1221: Output:
1222: Destroyed pdipm
1223: */
1224: PetscErrorCode TaoDestroy_PDIPM(Tao tao)
1225: {
1226: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1230: /* Freeing Vectors assocaiated with KKT (X) */
1231: VecDestroy(&pdipm->x); /* Solution x */
1232: VecDestroy(&pdipm->lambdae); /* Equality constraints lagrangian multiplier*/
1233: VecDestroy(&pdipm->lambdai); /* Inequality constraints lagrangian multiplier*/
1234: VecDestroy(&pdipm->z); /* Slack variables */
1235: VecDestroy(&pdipm->X); /* Big KKT system vector [x; lambdae; lambdai; z] */
1237: /* work vectors */
1238: VecDestroy(&pdipm->lambdae_xfixed);
1239: VecDestroy(&pdipm->lambdai_xb);
1241: /* Legrangian equality and inequality Vec */
1242: VecDestroy(&pdipm->ce); /* Vec of equality constraints */
1243: VecDestroy(&pdipm->ci); /* Vec of inequality constraints */
1245: /* Matrices */
1246: MatDestroy(&pdipm->Jce_xfixed);
1247: MatDestroy(&pdipm->Jci_xb); /* Jacobian of inequality constraints Jci = [tao->jacobian_inequality ; J(nxub); J(nxlb); J(nxbx)] */
1248: MatDestroy(&pdipm->K);
1250: /* Index Sets */
1251: if (pdipm->Nxub) {
1252: ISDestroy(&pdipm->isxub); /* Finite upper bound only -inf < x < ub */
1253: }
1255: if (pdipm->Nxlb) {
1256: ISDestroy(&pdipm->isxlb); /* Finite lower bound only lb <= x < inf */
1257: }
1259: if (pdipm->Nxfixed) {
1260: ISDestroy(&pdipm->isxfixed); /* Fixed variables lb = x = ub */
1261: }
1263: if (pdipm->Nxbox) {
1264: ISDestroy(&pdipm->isxbox); /* Boxed variables lb <= x <= ub */
1265: }
1267: if (pdipm->Nxfree) {
1268: ISDestroy(&pdipm->isxfree); /* Free variables -inf <= x <= inf */
1269: }
1271: if (pdipm->solve_reduced_kkt) {
1272: ISDestroy(&pdipm->is1);
1273: ISDestroy(&pdipm->is2);
1274: }
1276: /* SNES */
1277: SNESDestroy(&pdipm->snes); /* Nonlinear solver */
1278: PetscFree(pdipm->nce_all);
1279: MatDestroy(&pdipm->jac_equality_trans);
1280: MatDestroy(&pdipm->jac_inequality_trans);
1282: /* Destroy pdipm */
1283: PetscFree(tao->data); /* Holding locations of pdipm */
1285: /* Destroy Dual */
1286: VecDestroy(&tao->DE); /* equality dual */
1287: VecDestroy(&tao->DI); /* dinequality dual */
1288: return(0);
1289: }
1291: PetscErrorCode TaoSetFromOptions_PDIPM(PetscOptionItems *PetscOptionsObject,Tao tao)
1292: {
1293: TAO_PDIPM *pdipm = (TAO_PDIPM*)tao->data;
1297: PetscOptionsHead(PetscOptionsObject,"PDIPM method for constrained optimization");
1298: 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);
1299: 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);
1300: PetscOptionsBool("-tao_pdipm_solve_reduced_kkt","Solve reduced KKT system using Schur-complement",NULL,pdipm->solve_reduced_kkt,&pdipm->solve_reduced_kkt,NULL);
1301: PetscOptionsReal("-tao_pdipm_mu_update_factor","Update scalar for barrier parameter (mu) update",NULL,pdipm->mu_update_factor,&pdipm->mu_update_factor,NULL);
1302: PetscOptionsTail();
1303: return(0);
1304: }
1306: /*MC
1307: TAOPDIPM - Barrier-based primal-dual interior point algorithm for generally constrained optimization.
1309: Option Database Keys:
1310: + -tao_pdipm_push_init_lambdai - parameter to push initial dual variables away from bounds (> 0)
1311: . -tao_pdipm_push_init_slack - parameter to push initial slack variables away from bounds (> 0)
1312: - -tao_pdipm_mu_update_factor - update scalar for barrier parameter (mu) update (> 0)
1314: Level: beginner
1315: M*/
1316: PETSC_EXTERN PetscErrorCode TaoCreate_PDIPM(Tao tao)
1317: {
1318: TAO_PDIPM *pdipm;
1322: tao->ops->setup = TaoSetup_PDIPM;
1323: tao->ops->solve = TaoSolve_PDIPM;
1324: tao->ops->setfromoptions = TaoSetFromOptions_PDIPM;
1325: tao->ops->destroy = TaoDestroy_PDIPM;
1327: PetscNewLog(tao,&pdipm);
1328: tao->data = (void*)pdipm;
1330: pdipm->nx = pdipm->Nx = 0;
1331: pdipm->nxfixed = pdipm->Nxfixed = 0;
1332: pdipm->nxlb = pdipm->Nxlb = 0;
1333: pdipm->nxub = pdipm->Nxub = 0;
1334: pdipm->nxbox = pdipm->Nxbox = 0;
1335: pdipm->nxfree = pdipm->Nxfree = 0;
1337: pdipm->ng = pdipm->Ng = pdipm->nce = pdipm->Nce = 0;
1338: pdipm->nh = pdipm->Nh = pdipm->nci = pdipm->Nci = 0;
1339: pdipm->n = pdipm->N = 0;
1340: pdipm->mu = 1.0;
1341: pdipm->mu_update_factor = 0.1;
1343: pdipm->push_init_slack = 1.0;
1344: pdipm->push_init_lambdai = 1.0;
1345: pdipm->solve_reduced_kkt = PETSC_FALSE;
1347: /* Override default settings (unless already changed) */
1348: if (!tao->max_it_changed) tao->max_it = 200;
1349: if (!tao->max_funcs_changed) tao->max_funcs = 500;
1351: SNESCreate(((PetscObject)tao)->comm,&pdipm->snes);
1352: SNESSetOptionsPrefix(pdipm->snes,tao->hdr.prefix);
1353: SNESGetKSP(pdipm->snes,&tao->ksp);
1354: PetscObjectReference((PetscObject)tao->ksp);
1355: return(0);
1356: }