Actual source code: ex1.c
petsc-3.14.6 2021-03-30
1: /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
4: min f = (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: ---------------------------------------------------------------------- */
10: #include <petsctao.h>
12: static char help[]= "Solves constrained optimiztion problem using pdipm.\n\
13: Input parameters include:\n\
14: -tao_type pdipm : sets Tao solver\n\
15: -no_eq : removes the equaility constraints from the problem\n\
16: -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
17: -tao_cmonitor : convergence monitor with constraint norm \n\
18: -tao_view_solution : view exact solution at each itteration\n\
19: Note: external package superlu_dist is requried to run either for ipm or pdipm. Additionally This is designed for a maximum of 2 processors, the code will error if size > 2.\n\n";
21: /*
22: User-defined application context - contains data needed by the
23: application-provided call-back routines, FormFunction(),
24: FormGradient(), and FormHessian().
25: */
27: /*
28: x,d in R^n
29: f in R
30: bin in R^mi
31: beq in R^me
32: Aeq in R^(me x n)
33: Ain in R^(mi x n)
34: H in R^(n x n)
35: min f=(1/2)*x'*H*x + d'*x
36: s.t. Aeq*x == beq
37: Ain*x >= bin
38: */
39: typedef struct {
40: PetscInt n; /* Global length of x */
41: PetscInt ne; /* Global number of equality constraints */
42: PetscInt ni; /* Global number of inequality constraints */
43: PetscBool noeqflag;
44: Vec x,xl,xu;
45: Vec ce,ci,bl,bu,Xseq;
46: Mat Ae,Ai,H;
47: VecScatter scat;
48: } AppCtx;
50: /* -------- User-defined Routines --------- */
51: PetscErrorCode InitializeProblem(AppCtx *);
52: PetscErrorCode DestroyProblem(AppCtx *);
53: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
54: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
55: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
56: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
57: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
58: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
60: PetscErrorCode main(int argc,char **argv)
61: {
63: Tao tao;
64: KSP ksp;
65: PC pc;
66: AppCtx user; /* application context */
67: Vec x;
68: PetscMPIInt rank;
70: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
71: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
72: if (rank>1){
73: SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n");
74: }
76: PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");
77: InitializeProblem(&user); /* sets up problem, function below */
78: TaoCreate(PETSC_COMM_WORLD,&tao);
79: TaoSetType(tao,TAOPDIPM);
80: TaoSetInitialVector(tao,user.x); /* gets solution vector from problem */
81: TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
82: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
84: if (!user.noeqflag){
85: TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
86: }
87: TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);
89: if (!user.noeqflag){
90: TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
91: }
92: TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
93: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
94: TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
95: TaoSetConstraintTolerances(tao,1.e-6,1.e-6);
97: TaoGetKSP(tao,&ksp);
98: KSPGetPC(ksp,&pc);
99: PCSetType(pc,PCLU);
100: /*
101: This algorithm produces matrices with zeros along the diagonal therefore we use
102: SuperLU_DIST which provides shift to the diagonal
103: */
104: PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST); /* requires superlu_dist to solve pdipm */
105: KSPSetType(ksp,KSPPREONLY);
106: KSPSetFromOptions(ksp);
108: TaoSetFromOptions(tao);
110: TaoSolve(tao);
111: TaoGetSolutionVector(tao,&x);
112: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
114: /* Free objects */
115: DestroyProblem(&user);
116: TaoDestroy(&tao);
117: PetscFinalize();
118: return ierr;
119: }
121: PetscErrorCode InitializeProblem(AppCtx *user)
122: {
124: PetscMPIInt size;
125: PetscMPIInt rank;
126: PetscInt nloc,neloc,niloc;
129: MPI_Comm_size(PETSC_COMM_WORLD,&size);
130: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
131: user->noeqflag = PETSC_FALSE;
132: 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: VecCreateMPI(PETSC_COMM_WORLD,nloc,user->n,&user->x);
141: VecSet(user->x,0);
143: /* create and set lower and upper bound vectors */
144: VecDuplicate(user->x,&user->xl);
145: VecDuplicate(user->x,&user->xu);
146: VecSet(user->xl,-1.0);
147: VecSet(user->xu,2.0);
149: /* create scater to zero */
150: VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);
152: user->ne = 1;
153: neloc = (rank==0)?user->ne:0;
155: if (!user->noeqflag){
156: VecCreateMPI(PETSC_COMM_WORLD,neloc,user->ne,&user->ce); /* a 1x1 vec for equality constraints */
157: }
158: user->ni = 2;
159: niloc = (rank==0)?user->ni:0;
160: VecCreateMPI(PETSC_COMM_WORLD,niloc,user->ni,&user->ci); /* a 2x1 vec for inequality constraints */
162: /* nexn & nixn matricies for equaly and inequalty constriants */
163: if (!user->noeqflag){
164: MatCreate(PETSC_COMM_WORLD,&user->Ae);
165: MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
166: MatSetUp(user->Ae);
167: MatSetFromOptions(user->Ae);
168: }
170: MatCreate(PETSC_COMM_WORLD,&user->Ai);
171: MatCreate(PETSC_COMM_WORLD,&user->H);
173: MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
174: MatSetSizes(user->H,nloc,nloc,user->n,user->n);
176: MatSetUp(user->Ai);
177: MatSetUp(user->H);
179: MatSetFromOptions(user->Ai);
180: MatSetFromOptions(user->H);
181: return(0);
182: }
184: PetscErrorCode DestroyProblem(AppCtx *user)
185: {
189: if (!user->noeqflag){
190: MatDestroy(&user->Ae);
191: }
192: MatDestroy(&user->Ai);
193: MatDestroy(&user->H);
195: VecDestroy(&user->x);
196: if (!user->noeqflag){
197: VecDestroy(&user->ce);
198: }
199: VecDestroy(&user->ci);
200: VecDestroy(&user->xl);
201: VecDestroy(&user->xu);
202: VecDestroy(&user->Xseq);
203: VecScatterDestroy(&user->scat);
204: return(0);
205: }
207: /*
208: f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
209: G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0]
210: */
211: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
212: {
213: PetscScalar g;
214: const PetscScalar *x;
215: MPI_Comm comm;
216: PetscMPIInt rank;
217: PetscErrorCode ierr;
218: PetscReal fin;
219: AppCtx *user=(AppCtx*)ctx;
220: Vec Xseq=user->Xseq;
221: VecScatter scat=user->scat;
224: PetscObjectGetComm((PetscObject)tao,&comm);
225: MPI_Comm_rank(comm,&rank);
227: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
228: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
230: fin = 0.0;
231: if (!rank) {
232: VecGetArrayRead(Xseq,&x);
233: fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
234: g = 2.0*(x[0]-2.0) - 2.0;
235: VecSetValue(G,0,g,INSERT_VALUES);
236: g = 2.0*(x[1]-2.0) - 2.0;
237: VecSetValue(G,1,g,INSERT_VALUES);
238: VecRestoreArrayRead(Xseq,&x);
239: }
240: MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
241: VecAssemblyBegin(G);
242: VecAssemblyEnd(G);
243: return(0);
244: }
246: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
247: {
248: AppCtx *user=(AppCtx*)ctx;
249: Vec DE,DI;
250: const PetscScalar *de, *di;
251: PetscInt zero=0,one=1;
252: PetscScalar two=2.0;
253: PetscScalar val=0.0;
254: Vec Deseq,Diseq=user->Xseq;
255: VecScatter Descat,Discat=user->scat;
256: PetscMPIInt rank;
257: MPI_Comm comm;
258: PetscErrorCode ierr;
261: TaoGetDualVariables(tao,&DE,&DI);
263: PetscObjectGetComm((PetscObject)tao,&comm);
264: MPI_Comm_rank(comm,&rank);
266: if (!user->noeqflag){
267: VecScatterCreateToZero(DE,&Descat,&Deseq);
268: VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
269: VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
270: }
272: VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
273: VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
275: if (!rank){
276: if (!user->noeqflag){
277: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
278: }
280: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
281: if (!user->noeqflag){
282: val = 2.0 * (1 + de[0] + di[0] - di[1]);
283: VecRestoreArrayRead(Deseq,&de);
284: }else{
285: val = 2.0 * (1 + di[0] - di[1]);
286: }
287: VecRestoreArrayRead(Diseq,&di);
288: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
289: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
290: }
292: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
293: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
294: if (!user->noeqflag){
295: VecScatterDestroy(&Descat);
296: VecDestroy(&Deseq);
297: }
298: return(0);
299: }
301: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
302: {
303: const PetscScalar *x;
304: PetscScalar ci;
305: PetscErrorCode ierr;
306: MPI_Comm comm;
307: PetscMPIInt rank;
308: AppCtx *user=(AppCtx*)ctx;
309: Vec Xseq=user->Xseq;
310: VecScatter scat=user->scat;
313: PetscObjectGetComm((PetscObject)tao,&comm);
314: MPI_Comm_rank(comm,&rank);
316: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
317: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
319: if (!rank) {
320: VecGetArrayRead(Xseq,&x);
321: ci = x[0]*x[0] - x[1];
322: VecSetValue(CI,0,ci,INSERT_VALUES);
323: ci = -x[0]*x[0] + x[1] + 1.0;
324: VecSetValue(CI,1,ci,INSERT_VALUES);
325: VecRestoreArrayRead(Xseq,&x);
326: }
327: VecAssemblyBegin(CI);
328: VecAssemblyEnd(CI);
329: return(0);
330: }
332: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
333: {
334: const PetscScalar *x;
335: PetscScalar ce;
336: PetscErrorCode ierr;
337: MPI_Comm comm;
338: PetscMPIInt rank;
339: AppCtx *user=(AppCtx*)ctx;
340: Vec Xseq=user->Xseq;
341: VecScatter scat=user->scat;
344: PetscObjectGetComm((PetscObject)tao,&comm);
345: MPI_Comm_rank(comm,&rank);
347: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
348: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
350: if (!rank) {
351: VecGetArrayRead(Xseq,&x);
352: ce = x[0]*x[0] + x[1] - 2.0;
353: VecSetValue(CE,0,ce,INSERT_VALUES);
354: VecRestoreArrayRead(Xseq,&x);
355: }
356: VecAssemblyBegin(CE);
357: VecAssemblyEnd(CE);
358: return(0);
359: }
361: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
362: {
363: AppCtx *user=(AppCtx*)ctx;
364: PetscInt cols[2],min,max,i;
365: PetscScalar vals[2];
366: const PetscScalar *x;
367: PetscErrorCode ierr;
368: Vec Xseq=user->Xseq;
369: VecScatter scat=user->scat;
370: MPI_Comm comm;
371: PetscMPIInt rank;
374: PetscObjectGetComm((PetscObject)tao,&comm);
375: MPI_Comm_rank(comm,&rank);
376: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
377: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
379: VecGetArrayRead(Xseq,&x);
380: MatGetOwnershipRange(JI,&min,&max);
382: cols[0] = 0; cols[1] = 1;
383: for (i=min;i<max;i++) {
384: if (i==0){
385: vals[0] = +2*x[0]; vals[1] = -1.0;
386: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
387: }
388: if (i==1) {
389: vals[0] = -2*x[0]; vals[1] = +1.0;
390: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
391: }
392: }
393: VecRestoreArrayRead(Xseq,&x);
395: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
396: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
397: return(0);
398: }
400: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
401: {
402: PetscInt rows[2];
403: PetscScalar vals[2];
404: const PetscScalar *x;
405: PetscMPIInt rank;
406: MPI_Comm comm;
407: PetscErrorCode ierr;
410: PetscObjectGetComm((PetscObject)tao,&comm);
411: MPI_Comm_rank(comm,&rank);
413: if (!rank) {
414: VecGetArrayRead(X,&x);
415: rows[0] = 0; rows[1] = 1;
416: vals[0] = 2*x[0]; vals[1] = 1.0;
417: MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
418: VecRestoreArrayRead(X,&x);
419: }
420: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
421: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
422: return(0);
423: }
426: /*TEST
428: build:
429: requires: !complex !define(PETSC_USE_CXX)
431: test:
432: requires: superlu_dist
433: args: -tao_converged_reason
435: test:
436: suffix: 2
437: requires: superlu_dist
438: nsize: 2
439: args: -tao_converged_reason
441: test:
442: suffix: 3
443: requires: superlu_dist
444: args: -tao_converged_reason -no_eq
446: test:
447: suffix: 4
448: requires: superlu_dist
449: nsize: 2
450: args: -tao_converged_reason -no_eq
452: TEST*/