Actual source code: ex1.c
1: /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
4: min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
5: s.t. x0^2 + x1 - 2 = 0
6: 0 <= x0^2 - x1 <= 1
7: -1 <= x0, x1 <= 2
8: -->
9: g(x) = 0
10: h(x) >= 0
11: -1 <= x0, x1 <= 2
12: where
13: g(x) = x0^2 + x1 - 2
14: h(x) = [x0^2 - x1
15: 1 -(x0^2 - x1)]
16: ---------------------------------------------------------------------- */
18: #include <petsctao.h>
20: static char help[]= "Solves constrained optimiztion problem using pdipm.\n\
21: Input parameters include:\n\
22: -tao_type pdipm : sets Tao solver\n\
23: -no_eq : removes the equaility constraints from the problem\n\
24: -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
25: -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
26: -tao_cmonitor : convergence monitor with constraint norm \n\
27: -tao_view_solution : view exact solution at each itteration\n\
28: Note: external package MUMPS is required to run pdipm. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";
30: /*
31: User-defined application context - contains data needed by the application
32: */
33: typedef struct {
34: PetscInt n; /* Global length of x */
35: PetscInt ne; /* Global number of equality constraints */
36: PetscInt ni; /* Global number of inequality constraints */
37: PetscBool noeqflag;
38: Vec x,xl,xu;
39: Vec ce,ci,bl,bu,Xseq;
40: Mat Ae,Ai,H;
41: VecScatter scat;
42: } AppCtx;
44: /* -------- User-defined Routines --------- */
45: PetscErrorCode InitializeProblem(AppCtx *);
46: PetscErrorCode DestroyProblem(AppCtx *);
47: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
48: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
49: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
50: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
51: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
52: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
54: PetscErrorCode main(int argc,char **argv)
55: {
57: Tao tao;
58: KSP ksp;
59: PC pc;
60: AppCtx user; /* application context */
61: Vec x;
62: PetscMPIInt rank;
63: TaoType type;
64: PetscBool pdipm;
66: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
67: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
68: if (rank>1){
69: SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n");
70: }
72: PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");
73: InitializeProblem(&user); /* sets up problem, function below */
74: TaoCreate(PETSC_COMM_WORLD,&tao);
75: TaoSetType(tao,TAOPDIPM);
76: TaoSetInitialVector(tao,user.x); /* gets solution vector from problem */
77: TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
78: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
80: if (!user.noeqflag){
81: TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
82: }
83: TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);
84: if (!user.noeqflag){
85: TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
86: }
87: TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
88: TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
89: TaoSetConstraintTolerances(tao,1.e-6,1.e-6);
91: TaoGetKSP(tao,&ksp);
92: KSPGetPC(ksp,&pc);
93: PCSetType(pc,PCCHOLESKY);
94: /*
95: This algorithm produces matrices with zeros along the diagonal therefore we use
96: MUMPS which provides solver for indefinite matrices
97: */
98: PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); /* requires mumps to solve pdipm */
99: KSPSetType(ksp,KSPPREONLY);
100: KSPSetFromOptions(ksp);
102: TaoSetFromOptions(tao);
103: TaoGetType(tao,&type);
104: PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);
105: if (pdipm) {
106: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
107: }
109: TaoSolve(tao);
110: TaoGetSolutionVector(tao,&x);
111: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
113: /* Free objects */
114: DestroyProblem(&user);
115: TaoDestroy(&tao);
116: PetscFinalize();
117: return ierr;
118: }
120: PetscErrorCode InitializeProblem(AppCtx *user)
121: {
123: PetscMPIInt size;
124: PetscMPIInt rank;
125: PetscInt nloc,neloc,niloc;
128: MPI_Comm_size(PETSC_COMM_WORLD,&size);
129: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
130: user->noeqflag = PETSC_FALSE;
131: PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);
133: if (!user->noeqflag) {
134: PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
135: }
137: /* create vector x and set initial values */
138: user->n = 2; /* global length */
139: nloc = (rank==0)?user->n:0;
140: VecCreate(PETSC_COMM_WORLD,&user->x);
141: VecSetSizes(user->x,nloc,user->n);
142: VecSetFromOptions(user->x);
143: VecSet(user->x,0);
145: /* create and set lower and upper bound vectors */
146: VecDuplicate(user->x,&user->xl);
147: VecDuplicate(user->x,&user->xu);
148: VecSet(user->xl,-1.0);
149: VecSet(user->xu,2.0);
151: /* create scater to zero */
152: VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);
154: user->ne = 1;
155: user->ni = 2;
156: neloc = (rank==0)?user->ne:0;
157: niloc = (rank==0)?user->ni:0;
159: if (!user->noeqflag){
160: VecCreate(PETSC_COMM_WORLD,&user->ce); /* a 1x1 vec for equality constraints */
161: VecSetSizes(user->ce,neloc,user->ne);
162: VecSetFromOptions(user->ce);
163: VecSetUp(user->ce);
164: }
166: VecCreate(PETSC_COMM_WORLD,&user->ci); /* a 2x1 vec for inequality constraints */
167: VecSetSizes(user->ci,niloc,user->ni);
168: VecSetFromOptions(user->ci);
169: VecSetUp(user->ci);
171: /* nexn & nixn matricies for equaly and inequalty constriants */
172: if (!user->noeqflag){
173: MatCreate(PETSC_COMM_WORLD,&user->Ae);
174: MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
175: MatSetFromOptions(user->Ae);
176: MatSetUp(user->Ae);
177: }
179: MatCreate(PETSC_COMM_WORLD,&user->Ai);
180: MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
181: MatSetFromOptions(user->Ai);
182: MatSetUp(user->Ai);
184: MatCreate(PETSC_COMM_WORLD,&user->H);
185: MatSetSizes(user->H,nloc,nloc,user->n,user->n);
186: MatSetFromOptions(user->H);
187: MatSetUp(user->H);
188: return(0);
189: }
191: PetscErrorCode DestroyProblem(AppCtx *user)
192: {
196: if (!user->noeqflag){
197: MatDestroy(&user->Ae);
198: }
199: MatDestroy(&user->Ai);
200: MatDestroy(&user->H);
202: VecDestroy(&user->x);
203: if (!user->noeqflag){
204: VecDestroy(&user->ce);
205: }
206: VecDestroy(&user->ci);
207: VecDestroy(&user->xl);
208: VecDestroy(&user->xu);
209: VecDestroy(&user->Xseq);
210: VecScatterDestroy(&user->scat);
211: return(0);
212: }
214: /* Evaluate
215: f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
216: G = grad f = [2*(x0 - 2) - 2;
217: 2*(x1 - 2) - 2]
218: */
219: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
220: {
221: PetscScalar g;
222: const PetscScalar *x;
223: MPI_Comm comm;
224: PetscMPIInt rank;
225: PetscErrorCode ierr;
226: PetscReal fin;
227: AppCtx *user=(AppCtx*)ctx;
228: Vec Xseq=user->Xseq;
229: VecScatter scat=user->scat;
232: PetscObjectGetComm((PetscObject)tao,&comm);
233: MPI_Comm_rank(comm,&rank);
235: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
236: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
238: fin = 0.0;
239: if (!rank) {
240: VecGetArrayRead(Xseq,&x);
241: fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
242: g = 2.0*(x[0]-2.0) - 2.0;
243: VecSetValue(G,0,g,INSERT_VALUES);
244: g = 2.0*(x[1]-2.0) - 2.0;
245: VecSetValue(G,1,g,INSERT_VALUES);
246: VecRestoreArrayRead(Xseq,&x);
247: }
248: MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
249: VecAssemblyBegin(G);
250: VecAssemblyEnd(G);
251: return(0);
252: }
254: /* Evaluate
255: H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
256: = [ 2*(1+de[0]-di[0]+di[1]), 0;
257: 0, 2]
258: */
259: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
260: {
261: AppCtx *user=(AppCtx*)ctx;
262: Vec DE,DI;
263: const PetscScalar *de,*di;
264: PetscInt zero=0,one=1;
265: PetscScalar two=2.0;
266: PetscScalar val=0.0;
267: Vec Deseq,Diseq=user->Xseq;
268: VecScatter Descat,Discat=user->scat;
269: PetscMPIInt rank;
270: MPI_Comm comm;
271: PetscErrorCode ierr;
274: TaoGetDualVariables(tao,&DE,&DI);
276: PetscObjectGetComm((PetscObject)tao,&comm);
277: MPI_Comm_rank(comm,&rank);
279: if (!user->noeqflag){
280: VecScatterCreateToZero(DE,&Descat,&Deseq);
281: VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
282: VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
283: }
284: VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
285: VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
287: if (!rank){
288: if (!user->noeqflag){
289: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
290: }
291: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
293: if (!user->noeqflag){
294: val = 2.0 * (1 + de[0] - di[0] + di[1]);
295: VecRestoreArrayRead(Deseq,&de);
296: VecRestoreArrayRead(Diseq,&di);
297: } else {
298: val = 2.0 * (1 - di[0] + di[1]);
299: }
300: VecRestoreArrayRead(Diseq,&di);
301: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
302: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
303: }
304: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
306: if (!user->noeqflag){
307: VecScatterDestroy(&Descat);
308: VecDestroy(&Deseq);
309: }
310: return(0);
311: }
313: /* Evaluate
314: h = [ x0^2 - x1;
315: 1 -(x0^2 - x1)]
316: */
317: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
318: {
319: const PetscScalar *x;
320: PetscScalar ci;
321: PetscErrorCode ierr;
322: MPI_Comm comm;
323: PetscMPIInt rank;
324: AppCtx *user=(AppCtx*)ctx;
325: Vec Xseq=user->Xseq;
326: VecScatter scat=user->scat;
329: PetscObjectGetComm((PetscObject)tao,&comm);
330: MPI_Comm_rank(comm,&rank);
332: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
333: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
335: if (!rank) {
336: VecGetArrayRead(Xseq,&x);
337: ci = x[0]*x[0] - x[1];
338: VecSetValue(CI,0,ci,INSERT_VALUES);
339: ci = -x[0]*x[0] + x[1] + 1.0;
340: VecSetValue(CI,1,ci,INSERT_VALUES);
341: VecRestoreArrayRead(Xseq,&x);
342: }
343: VecAssemblyBegin(CI);
344: VecAssemblyEnd(CI);
345: return(0);
346: }
348: /* Evaluate
349: g = [ x0^2 + x1 - 2]
350: */
351: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
352: {
353: const PetscScalar *x;
354: PetscScalar ce;
355: PetscErrorCode ierr;
356: MPI_Comm comm;
357: PetscMPIInt rank;
358: AppCtx *user=(AppCtx*)ctx;
359: Vec Xseq=user->Xseq;
360: VecScatter scat=user->scat;
363: PetscObjectGetComm((PetscObject)tao,&comm);
364: MPI_Comm_rank(comm,&rank);
366: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
367: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
369: if (!rank) {
370: VecGetArrayRead(Xseq,&x);
371: ce = x[0]*x[0] + x[1] - 2.0;
372: VecSetValue(CE,0,ce,INSERT_VALUES);
373: VecRestoreArrayRead(Xseq,&x);
374: }
375: VecAssemblyBegin(CE);
376: VecAssemblyEnd(CE);
377: return(0);
378: }
380: /*
381: grad h = [ 2*x0, -1;
382: -2*x0, 1]
383: */
384: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
385: {
386: AppCtx *user=(AppCtx*)ctx;
387: PetscInt cols[2],min,max,i;
388: PetscScalar vals[2];
389: const PetscScalar *x;
390: PetscErrorCode ierr;
391: Vec Xseq=user->Xseq;
392: VecScatter scat=user->scat;
393: MPI_Comm comm;
394: PetscMPIInt rank;
397: PetscObjectGetComm((PetscObject)tao,&comm);
398: MPI_Comm_rank(comm,&rank);
399: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
400: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
402: VecGetArrayRead(Xseq,&x);
403: MatGetOwnershipRange(JI,&min,&max);
405: cols[0] = 0; cols[1] = 1;
406: for (i=min;i<max;i++) {
407: if (i==0){
408: vals[0] = 2*x[0]; vals[1] = -1.0;
409: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
410: }
411: if (i==1) {
412: vals[0] = -2*x[0]; vals[1] = 1.0;
413: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
414: }
415: }
416: VecRestoreArrayRead(Xseq,&x);
418: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
419: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
420: return(0);
421: }
423: /*
424: grad g = [2*x0
425: 1.0 ]
426: */
427: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
428: {
429: PetscInt rows[2];
430: PetscScalar vals[2];
431: const PetscScalar *x;
432: PetscMPIInt rank;
433: MPI_Comm comm;
434: PetscErrorCode ierr;
437: PetscObjectGetComm((PetscObject)tao,&comm);
438: MPI_Comm_rank(comm,&rank);
440: if (!rank) {
441: VecGetArrayRead(X,&x);
442: rows[0] = 0; rows[1] = 1;
443: vals[0] = 2*x[0]; vals[1] = 1.0;
444: MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
445: VecRestoreArrayRead(X,&x);
446: }
447: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
448: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
449: return(0);
450: }
453: /*TEST
455: build:
456: requires: !complex !define(PETSC_USE_CXX) mumps
458: test:
459: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
461: test:
462: suffix: 2
463: nsize: 2
464: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
466: test:
467: suffix: 3
468: args: -tao_converged_reason -no_eq
470: test:
471: suffix: 4
472: nsize: 2
473: args: -tao_converged_reason -no_eq
475: test:
476: suffix: 5
477: args: -tao_cmonitor -tao_type almm
479: test:
480: suffix: 6
481: args: -tao_cmonitor -tao_type almm -tao_almm_type phr
483: test:
484: suffix: 7
485: nsize: 2
486: args: -tao_cmonitor -tao_type almm
488: test:
489: suffix: 8
490: nsize: 2
491: requires: cuda
492: args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
494: test:
495: suffix: 9
496: nsize: 2
497: args: -tao_cmonitor -tao_type almm -no_eq
499: test:
500: suffix: 10
501: nsize: 2
502: args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
504: TEST*/