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: -init_view : view initial object setup\n\
25: -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
26: -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
27: -tao_cmonitor : convergence monitor with constraint norm \n\
28: -tao_view_solution : view exact solution at each iteration\n\
29: Note: external package MUMPS is required to run pdipm in parallel. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";
31: /*
32: User-defined application context - contains data needed by the application
33: */
34: typedef struct {
35: PetscInt n; /* Global length of x */
36: PetscInt ne; /* Global number of equality constraints */
37: PetscInt ni; /* Global number of inequality constraints */
38: PetscBool noeqflag, initview;
39: Vec x,xl,xu;
40: Vec ce,ci,bl,bu,Xseq;
41: Mat Ae,Ai,H;
42: VecScatter scat;
43: } AppCtx;
45: /* -------- User-defined Routines --------- */
46: PetscErrorCode InitializeProblem(AppCtx *);
47: PetscErrorCode DestroyProblem(AppCtx *);
48: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
49: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
50: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
51: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
52: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
53: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
55: PetscErrorCode main(int argc,char **argv)
56: {
58: Tao tao;
59: KSP ksp;
60: PC pc;
61: AppCtx user; /* application context */
62: Vec x,G,CI,CE;
63: PetscMPIInt size;
64: TaoType type;
65: PetscReal f;
66: PetscBool pdipm;
68: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
69: MPI_Comm_size(PETSC_COMM_WORLD,&size);
70: if (size>2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"More than 2 processors detected. Example written to use max of 2 processors.\n");
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: #if defined(PETSC_HAVE_MUMPS)
99: PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); /* requires mumps to solve pdipm */
100: #else
101: if (size > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Requires an external package that supports parallel PCCHOLESKY, e.g., MUMPS.");
102: #endif
103: KSPSetType(ksp,KSPPREONLY);
104: KSPSetFromOptions(ksp);
106: TaoSetFromOptions(tao);
107: TaoGetType(tao,&type);
108: PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);
109: if (pdipm) {
110: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
111: if (user.initview) {
112: TaoSetUp(tao);
113: VecDuplicate(user.x, &G);
114: FormFunctionGradient(tao, user.x, &f, G, (void*)&user);
115: FormHessian(tao, user.x, user.H, user.H, (void*)&user);
116: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD);
117: PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f);
118: VecView(user.x, PETSC_VIEWER_STDOUT_WORLD);
119: PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f);
120: PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f);
121: VecView(G, PETSC_VIEWER_STDOUT_WORLD);
122: MatView(user.H, PETSC_VIEWER_STDOUT_WORLD);
123: VecDestroy(&G);
124: FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user);
125: MatCreateVecs(user.Ai, NULL, &CI);
126: FormInequalityConstraints(tao, user.x, CI, (void*)&user);
127: PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f);
128: VecView(CI, PETSC_VIEWER_STDOUT_WORLD);
129: MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD);
130: VecDestroy(&CI);
131: if (!user.noeqflag) {
132: FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user);
133: MatCreateVecs(user.Ae, NULL, &CE);
134: FormEqualityConstraints(tao, user.x, CE, (void*)&user);
135: PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f);
136: VecView(CE, PETSC_VIEWER_STDOUT_WORLD);
137: MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD);
138: VecDestroy(&CE);
139: }
140: PetscPrintf(PETSC_COMM_WORLD, "\n");
141: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD);
142: }
143: }
145: TaoSolve(tao);
146: TaoGetSolutionVector(tao,&x);
147: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
149: /* Free objects */
150: DestroyProblem(&user);
151: TaoDestroy(&tao);
152: PetscFinalize();
153: return ierr;
154: }
156: PetscErrorCode InitializeProblem(AppCtx *user)
157: {
159: PetscMPIInt size;
160: PetscMPIInt rank;
161: PetscInt nloc,neloc,niloc;
164: MPI_Comm_size(PETSC_COMM_WORLD,&size);
165: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
166: user->noeqflag = PETSC_FALSE;
167: PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);
168: user->initview = PETSC_FALSE;
169: PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL);
171: if (!user->noeqflag) {
172: /* Tell user the correct solution, not an error checking */
173: PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
174: }
176: /* create vector x and set initial values */
177: user->n = 2; /* global length */
178: nloc = (size==1)?user->n:1;
179: VecCreate(PETSC_COMM_WORLD,&user->x);
180: VecSetSizes(user->x,nloc,user->n);
181: VecSetFromOptions(user->x);
182: VecSet(user->x,0);
184: /* create and set lower and upper bound vectors */
185: VecDuplicate(user->x,&user->xl);
186: VecDuplicate(user->x,&user->xu);
187: VecSet(user->xl,-1.0);
188: VecSet(user->xu,2.0);
190: /* create scater to zero */
191: VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);
193: user->ne = 1;
194: user->ni = 2;
195: neloc = (rank==0)?user->ne:0;
196: niloc = (size==1)?user->ni:1;
198: if (!user->noeqflag) {
199: VecCreate(PETSC_COMM_WORLD,&user->ce); /* a 1x1 vec for equality constraints */
200: VecSetSizes(user->ce,neloc,user->ne);
201: VecSetFromOptions(user->ce);
202: VecSetUp(user->ce);
203: }
205: VecCreate(PETSC_COMM_WORLD,&user->ci); /* a 2x1 vec for inequality constraints */
206: VecSetSizes(user->ci,niloc,user->ni);
207: VecSetFromOptions(user->ci);
208: VecSetUp(user->ci);
210: /* nexn & nixn matricies for equally and inequalty constraints */
211: if (!user->noeqflag) {
212: MatCreate(PETSC_COMM_WORLD,&user->Ae);
213: MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
214: MatSetFromOptions(user->Ae);
215: MatSetUp(user->Ae);
216: }
218: MatCreate(PETSC_COMM_WORLD,&user->Ai);
219: MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
220: MatSetFromOptions(user->Ai);
221: MatSetUp(user->Ai);
223: MatCreate(PETSC_COMM_WORLD,&user->H);
224: MatSetSizes(user->H,nloc,nloc,user->n,user->n);
225: MatSetFromOptions(user->H);
226: MatSetUp(user->H);
227: return(0);
228: }
230: PetscErrorCode DestroyProblem(AppCtx *user)
231: {
235: if (!user->noeqflag) {
236: MatDestroy(&user->Ae);
237: }
238: MatDestroy(&user->Ai);
239: MatDestroy(&user->H);
241: VecDestroy(&user->x);
242: if (!user->noeqflag) {
243: VecDestroy(&user->ce);
244: }
245: VecDestroy(&user->ci);
246: VecDestroy(&user->xl);
247: VecDestroy(&user->xu);
248: VecDestroy(&user->Xseq);
249: VecScatterDestroy(&user->scat);
250: return(0);
251: }
253: /* Evaluate
254: f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
255: G = grad f = [2*(x0 - 2) - 2;
256: 2*(x1 - 2) - 2]
257: */
258: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
259: {
260: PetscScalar g;
261: const PetscScalar *x;
262: MPI_Comm comm;
263: PetscMPIInt rank;
264: PetscErrorCode ierr;
265: PetscReal fin;
266: AppCtx *user=(AppCtx*)ctx;
267: Vec Xseq=user->Xseq;
268: VecScatter scat=user->scat;
271: PetscObjectGetComm((PetscObject)tao,&comm);
272: MPI_Comm_rank(comm,&rank);
274: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
275: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
277: fin = 0.0;
278: if (rank == 0) {
279: VecGetArrayRead(Xseq,&x);
280: fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
281: g = 2.0*(x[0]-2.0) - 2.0;
282: VecSetValue(G,0,g,INSERT_VALUES);
283: g = 2.0*(x[1]-2.0) - 2.0;
284: VecSetValue(G,1,g,INSERT_VALUES);
285: VecRestoreArrayRead(Xseq,&x);
286: }
287: MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
288: VecAssemblyBegin(G);
289: VecAssemblyEnd(G);
290: return(0);
291: }
293: /* Evaluate
294: H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
295: = [ 2*(1+de[0]-di[0]+di[1]), 0;
296: 0, 2]
297: */
298: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
299: {
300: AppCtx *user=(AppCtx*)ctx;
301: Vec DE,DI;
302: const PetscScalar *de,*di;
303: PetscInt zero=0,one=1;
304: PetscScalar two=2.0;
305: PetscScalar val=0.0;
306: Vec Deseq,Diseq;
307: VecScatter Descat,Discat;
308: PetscMPIInt rank;
309: MPI_Comm comm;
310: PetscErrorCode ierr;
313: TaoGetDualVariables(tao,&DE,&DI);
315: PetscObjectGetComm((PetscObject)tao,&comm);
316: MPI_Comm_rank(comm,&rank);
318: if (!user->noeqflag) {
319: VecScatterCreateToZero(DE,&Descat,&Deseq);
320: VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
321: VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
322: }
323: VecScatterCreateToZero(DI,&Discat,&Diseq);
324: VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
325: VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
327: if (rank == 0) {
328: if (!user->noeqflag) {
329: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
330: }
331: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
333: if (!user->noeqflag) {
334: val = 2.0 * (1 + de[0] - di[0] + di[1]);
335: VecRestoreArrayRead(Deseq,&de);
336: VecRestoreArrayRead(Diseq,&di);
337: } else {
338: val = 2.0 * (1 - di[0] + di[1]);
339: }
340: VecRestoreArrayRead(Diseq,&di);
341: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
342: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
343: }
344: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
345: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
346: if (!user->noeqflag) {
347: VecScatterDestroy(&Descat);
348: VecDestroy(&Deseq);
349: }
350: VecScatterDestroy(&Discat);
351: VecDestroy(&Diseq);
352: return(0);
353: }
355: /* Evaluate
356: h = [ x0^2 - x1;
357: 1 -(x0^2 - x1)]
358: */
359: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
360: {
361: const PetscScalar *x;
362: PetscScalar ci;
363: PetscErrorCode ierr;
364: MPI_Comm comm;
365: PetscMPIInt rank;
366: AppCtx *user=(AppCtx*)ctx;
367: Vec Xseq=user->Xseq;
368: VecScatter scat=user->scat;
371: PetscObjectGetComm((PetscObject)tao,&comm);
372: MPI_Comm_rank(comm,&rank);
374: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
375: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
377: if (rank == 0) {
378: VecGetArrayRead(Xseq,&x);
379: ci = x[0]*x[0] - x[1];
380: VecSetValue(CI,0,ci,INSERT_VALUES);
381: ci = -x[0]*x[0] + x[1] + 1.0;
382: VecSetValue(CI,1,ci,INSERT_VALUES);
383: VecRestoreArrayRead(Xseq,&x);
384: }
385: VecAssemblyBegin(CI);
386: VecAssemblyEnd(CI);
387: return(0);
388: }
390: /* Evaluate
391: g = [ x0^2 + x1 - 2]
392: */
393: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
394: {
395: const PetscScalar *x;
396: PetscScalar ce;
397: PetscErrorCode ierr;
398: MPI_Comm comm;
399: PetscMPIInt rank;
400: AppCtx *user=(AppCtx*)ctx;
401: Vec Xseq=user->Xseq;
402: VecScatter scat=user->scat;
405: PetscObjectGetComm((PetscObject)tao,&comm);
406: MPI_Comm_rank(comm,&rank);
408: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
409: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
411: if (rank == 0) {
412: VecGetArrayRead(Xseq,&x);
413: ce = x[0]*x[0] + x[1] - 2.0;
414: VecSetValue(CE,0,ce,INSERT_VALUES);
415: VecRestoreArrayRead(Xseq,&x);
416: }
417: VecAssemblyBegin(CE);
418: VecAssemblyEnd(CE);
419: return(0);
420: }
422: /*
423: grad h = [ 2*x0, -1;
424: -2*x0, 1]
425: */
426: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
427: {
428: AppCtx *user=(AppCtx*)ctx;
429: PetscInt zero=0,one=1,cols[2];
430: PetscScalar vals[2];
431: const PetscScalar *x;
432: PetscErrorCode ierr;
433: Vec Xseq=user->Xseq;
434: VecScatter scat=user->scat;
435: MPI_Comm comm;
436: PetscMPIInt rank;
439: PetscObjectGetComm((PetscObject)tao,&comm);
440: MPI_Comm_rank(comm,&rank);
441: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
442: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
444: VecGetArrayRead(Xseq,&x);
445: if (!rank) {
446: cols[0] = 0; cols[1] = 1;
447: vals[0] = 2*x[0]; vals[1] = -1.0;
448: MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES);
449: vals[0] = -2*x[0]; vals[1] = 1.0;
450: MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES);
451: }
452: VecRestoreArrayRead(Xseq,&x);
453: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
454: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
455: return(0);
456: }
458: /*
459: grad g = [2*x0
460: 1.0 ]
461: */
462: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
463: {
464: PetscInt zero=0,cols[2];
465: PetscScalar vals[2];
466: const PetscScalar *x;
467: PetscMPIInt rank;
468: MPI_Comm comm;
469: PetscErrorCode ierr;
472: PetscObjectGetComm((PetscObject)tao,&comm);
473: MPI_Comm_rank(comm,&rank);
475: if (rank == 0) {
476: VecGetArrayRead(X,&x);
477: cols[0] = 0; cols[1] = 1;
478: vals[0] = 2*x[0]; vals[1] = 1.0;
479: MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES);
480: VecRestoreArrayRead(X,&x);
481: }
482: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
483: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
484: return(0);
485: }
487: /*TEST
489: build:
490: requires: !complex !defined(PETSC_USE_CXX)
492: test:
493: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
494: requires: mumps
496: test:
497: suffix: 2
498: nsize: 2
499: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
500: requires: mumps
502: test:
503: suffix: 3
504: args: -tao_converged_reason -no_eq
505: requires: mumps
507: test:
508: suffix: 4
509: nsize: 2
510: args: -tao_converged_reason -no_eq
511: requires: mumps
513: test:
514: suffix: 5
515: args: -tao_cmonitor -tao_type almm
516: requires: mumps
518: test:
519: suffix: 6
520: args: -tao_cmonitor -tao_type almm -tao_almm_type phr
521: requires: mumps
523: test:
524: suffix: 7
525: nsize: 2
526: requires: mumps
527: args: -tao_cmonitor -tao_type almm
529: test:
530: suffix: 8
531: nsize: 2
532: requires: cuda mumps
533: args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
535: test:
536: suffix: 9
537: nsize: 2
538: args: -tao_cmonitor -tao_type almm -no_eq
539: requires: mumps
541: test:
542: suffix: 10
543: nsize: 2
544: args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
545: requires: mumps
547: TEST*/