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: {
57: Tao tao;
58: KSP ksp;
59: PC pc;
60: AppCtx user; /* application context */
61: Vec x,G,CI,CE;
62: PetscMPIInt size;
63: TaoType type;
64: PetscReal f;
65: PetscBool pdipm;
67: PetscInitialize(&argc,&argv,(char*)0,help);
68: MPI_Comm_size(PETSC_COMM_WORLD,&size);
71: PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");
72: InitializeProblem(&user); /* sets up problem, function below */
73: TaoCreate(PETSC_COMM_WORLD,&tao);
74: TaoSetType(tao,TAOPDIPM);
75: TaoSetSolution(tao,user.x); /* gets solution vector from problem */
76: TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
77: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user);
79: if (!user.noeqflag) {
80: TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
81: }
82: TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);
83: if (!user.noeqflag) {
84: TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
85: }
86: TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
87: TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
88: TaoSetConstraintTolerances(tao,1.e-6,1.e-6);
90: TaoGetKSP(tao,&ksp);
91: KSPGetPC(ksp,&pc);
92: PCSetType(pc,PCCHOLESKY);
93: /*
94: This algorithm produces matrices with zeros along the diagonal therefore we use
95: MUMPS which provides solver for indefinite matrices
96: */
97: #if defined(PETSC_HAVE_MUMPS)
98: PCFactorSetMatSolverType(pc,MATSOLVERMUMPS); /* requires mumps to solve pdipm */
99: #else
101: #endif
102: KSPSetType(ksp,KSPPREONLY);
103: KSPSetFromOptions(ksp);
105: TaoSetFromOptions(tao);
106: TaoGetType(tao,&type);
107: PetscObjectTypeCompare((PetscObject)tao,TAOPDIPM,&pdipm);
108: if (pdipm) {
109: TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user);
110: if (user.initview) {
111: TaoSetUp(tao);
112: VecDuplicate(user.x, &G);
113: FormFunctionGradient(tao, user.x, &f, G, (void*)&user);
114: FormHessian(tao, user.x, user.H, user.H, (void*)&user);
115: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD);
116: PetscPrintf(PETSC_COMM_WORLD,"\nInitial point X:\n",f);
117: VecView(user.x, PETSC_VIEWER_STDOUT_WORLD);
118: PetscPrintf(PETSC_COMM_WORLD,"\nInitial objective f(x) = %g\n",f);
119: PetscPrintf(PETSC_COMM_WORLD,"\nInitial gradient and Hessian:\n",f);
120: VecView(G, PETSC_VIEWER_STDOUT_WORLD);
121: MatView(user.H, PETSC_VIEWER_STDOUT_WORLD);
122: VecDestroy(&G);
123: FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void*)&user);
124: MatCreateVecs(user.Ai, NULL, &CI);
125: FormInequalityConstraints(tao, user.x, CI, (void*)&user);
126: PetscPrintf(PETSC_COMM_WORLD,"\nInitial inequality constraints and Jacobian:\n",f);
127: VecView(CI, PETSC_VIEWER_STDOUT_WORLD);
128: MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD);
129: VecDestroy(&CI);
130: if (!user.noeqflag) {
131: FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void*)&user);
132: MatCreateVecs(user.Ae, NULL, &CE);
133: FormEqualityConstraints(tao, user.x, CE, (void*)&user);
134: PetscPrintf(PETSC_COMM_WORLD,"\nInitial equality constraints and Jacobian:\n",f);
135: VecView(CE, PETSC_VIEWER_STDOUT_WORLD);
136: MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD);
137: VecDestroy(&CE);
138: }
139: PetscPrintf(PETSC_COMM_WORLD, "\n");
140: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD);
141: }
142: }
144: TaoSolve(tao);
145: TaoGetSolution(tao,&x);
146: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
148: /* Free objects */
149: DestroyProblem(&user);
150: TaoDestroy(&tao);
151: PetscFinalize();
152: return 0;
153: }
155: PetscErrorCode InitializeProblem(AppCtx *user)
156: {
157: PetscMPIInt size;
158: PetscMPIInt rank;
159: PetscInt nloc,neloc,niloc;
161: MPI_Comm_size(PETSC_COMM_WORLD,&size);
162: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
163: user->noeqflag = PETSC_FALSE;
164: PetscOptionsGetBool(NULL,NULL,"-no_eq",&user->noeqflag,NULL);
165: user->initview = PETSC_FALSE;
166: PetscOptionsGetBool(NULL,NULL,"-init_view",&user->initview,NULL);
168: if (!user->noeqflag) {
169: /* Tell user the correct solution, not an error checking */
170: PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
171: }
173: /* create vector x and set initial values */
174: user->n = 2; /* global length */
175: nloc = (size==1)?user->n:1;
176: VecCreate(PETSC_COMM_WORLD,&user->x);
177: VecSetSizes(user->x,nloc,user->n);
178: VecSetFromOptions(user->x);
179: VecSet(user->x,0);
181: /* create and set lower and upper bound vectors */
182: VecDuplicate(user->x,&user->xl);
183: VecDuplicate(user->x,&user->xu);
184: VecSet(user->xl,-1.0);
185: VecSet(user->xu,2.0);
187: /* create scater to zero */
188: VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);
190: user->ne = 1;
191: user->ni = 2;
192: neloc = (rank==0)?user->ne:0;
193: niloc = (size==1)?user->ni:1;
195: if (!user->noeqflag) {
196: VecCreate(PETSC_COMM_WORLD,&user->ce); /* a 1x1 vec for equality constraints */
197: VecSetSizes(user->ce,neloc,user->ne);
198: VecSetFromOptions(user->ce);
199: VecSetUp(user->ce);
200: }
202: VecCreate(PETSC_COMM_WORLD,&user->ci); /* a 2x1 vec for inequality constraints */
203: VecSetSizes(user->ci,niloc,user->ni);
204: VecSetFromOptions(user->ci);
205: VecSetUp(user->ci);
207: /* nexn & nixn matricies for equally and inequalty constraints */
208: if (!user->noeqflag) {
209: MatCreate(PETSC_COMM_WORLD,&user->Ae);
210: MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
211: MatSetFromOptions(user->Ae);
212: MatSetUp(user->Ae);
213: }
215: MatCreate(PETSC_COMM_WORLD,&user->Ai);
216: MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
217: MatSetFromOptions(user->Ai);
218: MatSetUp(user->Ai);
220: MatCreate(PETSC_COMM_WORLD,&user->H);
221: MatSetSizes(user->H,nloc,nloc,user->n,user->n);
222: MatSetFromOptions(user->H);
223: MatSetUp(user->H);
224: return 0;
225: }
227: PetscErrorCode DestroyProblem(AppCtx *user)
228: {
229: if (!user->noeqflag) {
230: MatDestroy(&user->Ae);
231: }
232: MatDestroy(&user->Ai);
233: MatDestroy(&user->H);
235: VecDestroy(&user->x);
236: if (!user->noeqflag) {
237: VecDestroy(&user->ce);
238: }
239: VecDestroy(&user->ci);
240: VecDestroy(&user->xl);
241: VecDestroy(&user->xu);
242: VecDestroy(&user->Xseq);
243: VecScatterDestroy(&user->scat);
244: return 0;
245: }
247: /* Evaluate
248: f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
249: G = grad f = [2*(x0 - 2) - 2;
250: 2*(x1 - 2) - 2]
251: */
252: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
253: {
254: PetscScalar g;
255: const PetscScalar *x;
256: MPI_Comm comm;
257: PetscMPIInt rank;
258: PetscReal fin;
259: AppCtx *user=(AppCtx*)ctx;
260: Vec Xseq=user->Xseq;
261: VecScatter scat=user->scat;
263: PetscObjectGetComm((PetscObject)tao,&comm);
264: MPI_Comm_rank(comm,&rank);
266: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
267: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
269: fin = 0.0;
270: if (rank == 0) {
271: VecGetArrayRead(Xseq,&x);
272: fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
273: g = 2.0*(x[0]-2.0) - 2.0;
274: VecSetValue(G,0,g,INSERT_VALUES);
275: g = 2.0*(x[1]-2.0) - 2.0;
276: VecSetValue(G,1,g,INSERT_VALUES);
277: VecRestoreArrayRead(Xseq,&x);
278: }
279: MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
280: VecAssemblyBegin(G);
281: VecAssemblyEnd(G);
282: return 0;
283: }
285: /* Evaluate
286: H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
287: = [ 2*(1+de[0]-di[0]+di[1]), 0;
288: 0, 2]
289: */
290: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
291: {
292: AppCtx *user=(AppCtx*)ctx;
293: Vec DE,DI;
294: const PetscScalar *de,*di;
295: PetscInt zero=0,one=1;
296: PetscScalar two=2.0;
297: PetscScalar val=0.0;
298: Vec Deseq,Diseq;
299: VecScatter Descat,Discat;
300: PetscMPIInt rank;
301: MPI_Comm comm;
303: TaoGetDualVariables(tao,&DE,&DI);
305: PetscObjectGetComm((PetscObject)tao,&comm);
306: MPI_Comm_rank(comm,&rank);
308: if (!user->noeqflag) {
309: VecScatterCreateToZero(DE,&Descat,&Deseq);
310: VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
311: VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
312: }
313: VecScatterCreateToZero(DI,&Discat,&Diseq);
314: VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
315: VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
317: if (rank == 0) {
318: if (!user->noeqflag) {
319: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
320: }
321: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
323: if (!user->noeqflag) {
324: val = 2.0 * (1 + de[0] - di[0] + di[1]);
325: VecRestoreArrayRead(Deseq,&de);
326: VecRestoreArrayRead(Diseq,&di);
327: } else {
328: val = 2.0 * (1 - di[0] + di[1]);
329: }
330: VecRestoreArrayRead(Diseq,&di);
331: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
332: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
333: }
334: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
335: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
336: if (!user->noeqflag) {
337: VecScatterDestroy(&Descat);
338: VecDestroy(&Deseq);
339: }
340: VecScatterDestroy(&Discat);
341: VecDestroy(&Diseq);
342: return 0;
343: }
345: /* Evaluate
346: h = [ x0^2 - x1;
347: 1 -(x0^2 - x1)]
348: */
349: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
350: {
351: const PetscScalar *x;
352: PetscScalar ci;
353: MPI_Comm comm;
354: PetscMPIInt rank;
355: AppCtx *user=(AppCtx*)ctx;
356: Vec Xseq=user->Xseq;
357: VecScatter scat=user->scat;
359: PetscObjectGetComm((PetscObject)tao,&comm);
360: MPI_Comm_rank(comm,&rank);
362: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
363: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
365: if (rank == 0) {
366: VecGetArrayRead(Xseq,&x);
367: ci = x[0]*x[0] - x[1];
368: VecSetValue(CI,0,ci,INSERT_VALUES);
369: ci = -x[0]*x[0] + x[1] + 1.0;
370: VecSetValue(CI,1,ci,INSERT_VALUES);
371: VecRestoreArrayRead(Xseq,&x);
372: }
373: VecAssemblyBegin(CI);
374: VecAssemblyEnd(CI);
375: return 0;
376: }
378: /* Evaluate
379: g = [ x0^2 + x1 - 2]
380: */
381: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
382: {
383: const PetscScalar *x;
384: PetscScalar ce;
385: MPI_Comm comm;
386: PetscMPIInt rank;
387: AppCtx *user=(AppCtx*)ctx;
388: Vec Xseq=user->Xseq;
389: VecScatter scat=user->scat;
391: PetscObjectGetComm((PetscObject)tao,&comm);
392: MPI_Comm_rank(comm,&rank);
394: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
395: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
397: if (rank == 0) {
398: VecGetArrayRead(Xseq,&x);
399: ce = x[0]*x[0] + x[1] - 2.0;
400: VecSetValue(CE,0,ce,INSERT_VALUES);
401: VecRestoreArrayRead(Xseq,&x);
402: }
403: VecAssemblyBegin(CE);
404: VecAssemblyEnd(CE);
405: return 0;
406: }
408: /*
409: grad h = [ 2*x0, -1;
410: -2*x0, 1]
411: */
412: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
413: {
414: AppCtx *user=(AppCtx*)ctx;
415: PetscInt zero=0,one=1,cols[2];
416: PetscScalar vals[2];
417: const PetscScalar *x;
418: Vec Xseq=user->Xseq;
419: VecScatter scat=user->scat;
420: MPI_Comm comm;
421: PetscMPIInt rank;
423: PetscObjectGetComm((PetscObject)tao,&comm);
424: MPI_Comm_rank(comm,&rank);
425: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
426: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
428: VecGetArrayRead(Xseq,&x);
429: if (!rank) {
430: cols[0] = 0; cols[1] = 1;
431: vals[0] = 2*x[0]; vals[1] = -1.0;
432: MatSetValues(JI,1,&zero,2,cols,vals,INSERT_VALUES);
433: vals[0] = -2*x[0]; vals[1] = 1.0;
434: MatSetValues(JI,1,&one,2,cols,vals,INSERT_VALUES);
435: }
436: VecRestoreArrayRead(Xseq,&x);
437: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
438: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
439: return 0;
440: }
442: /*
443: grad g = [2*x0
444: 1.0 ]
445: */
446: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
447: {
448: PetscInt zero=0,cols[2];
449: PetscScalar vals[2];
450: const PetscScalar *x;
451: PetscMPIInt rank;
452: MPI_Comm comm;
454: PetscObjectGetComm((PetscObject)tao,&comm);
455: MPI_Comm_rank(comm,&rank);
457: if (rank == 0) {
458: VecGetArrayRead(X,&x);
459: cols[0] = 0; cols[1] = 1;
460: vals[0] = 2*x[0]; vals[1] = 1.0;
461: MatSetValues(JE,1,&zero,2,cols,vals,INSERT_VALUES);
462: VecRestoreArrayRead(X,&x);
463: }
464: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
465: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
466: return 0;
467: }
469: /*TEST
471: build:
472: requires: !complex !defined(PETSC_USE_CXX)
474: test:
475: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
476: requires: mumps
478: test:
479: suffix: 2
480: nsize: 2
481: args: -tao_converged_reason -tao_pdipm_kkt_shift_pd
482: requires: mumps
484: test:
485: suffix: 3
486: args: -tao_converged_reason -no_eq
487: requires: mumps
489: test:
490: suffix: 4
491: nsize: 2
492: args: -tao_converged_reason -no_eq
493: requires: mumps
495: test:
496: suffix: 5
497: args: -tao_cmonitor -tao_type almm
498: requires: mumps
500: test:
501: suffix: 6
502: args: -tao_cmonitor -tao_type almm -tao_almm_type phr
503: requires: mumps
505: test:
506: suffix: 7
507: nsize: 2
508: requires: mumps
509: args: -tao_cmonitor -tao_type almm
511: test:
512: suffix: 8
513: nsize: 2
514: requires: cuda mumps
515: args: -tao_cmonitor -tao_type almm -vec_type cuda -mat_type aijcusparse
517: test:
518: suffix: 9
519: nsize: 2
520: args: -tao_cmonitor -tao_type almm -no_eq
521: requires: mumps
523: test:
524: suffix: 10
525: nsize: 2
526: args: -tao_cmonitor -tao_type almm -tao_almm_type phr -no_eq
527: requires: mumps
529: TEST*/