Actual source code: ex1.c
petsc-3.13.4 2020-08-01
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: -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
16: -tao_cmonitor : convergence monitor with constraint norm \n\
17: -tao_view_solution : view exact solution at each itteration\n\
18: 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";
20: /*
21: User-defined Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
22: Section 1.5 Writing Application Codes with PETSc-provided call-back routines, FormFunction(),
23: FormGradient(), and FormHessian().
24: */
26: /*
27: x,d in R^n
28: f in R
29: bin in R^mi
30: beq in R^me
31: Aeq in R^(me x n)
32: Ain in R^(mi x n)
33: H in R^(n x n)
34: min f=(1/2)*x'*H*x + d'*x
35: s.t. Aeq*x == beq
36: Ain*x >= bin
37: */
38: typedef struct {
39: PetscInt n; /* Global length of x */
40: PetscInt ne; /* Global number of equality constraints */
41: PetscInt ni; /* Global number of inequality constraints */
42: Vec x,xl,xu;
43: Vec ce,ci,bl,bu,Xseq;
44: Mat Ae,Ai,H;
45: VecScatter scat;
46: } AppCtx;
48: /* -------- User-defined Routines --------- */
49: PetscErrorCode InitializeProblem(AppCtx *);
50: PetscErrorCode DestroyProblem(AppCtx *);
51: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
52: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
53: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
54: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
55: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
56: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
58: PetscErrorCode main(int argc,char **argv)
59: {
61: Tao tao;
62: KSP ksp;
63: PC pc;
64: AppCtx user; /* Section 1.5 Writing Application Codes with PETSc context */
65: Vec x;
66: PetscMPIInt rank;
68: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
69: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
70: if (rank>1){
71: SETERRQ(PETSC_COMM_SELF,1,"More than 2 processors detected. Example written to use max of 2 processors.\n");
72: }
74: PetscPrintf(PETSC_COMM_WORLD,"---- Constrained Problem -----\n");
75: PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
76: InitializeProblem(&user); /* sets up problem, function below */
77: TaoCreate(PETSC_COMM_WORLD,&tao);
78: TaoSetType(tao,TAOPDIPM);
79: TaoSetInitialVector(tao,user.x); /* gets solution vector from problem */
80: TaoSetVariableBounds(tao,user.xl,user.xu); /* sets lower upper bounds from given solution */
81: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
83: TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
84: TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);
86: TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user); /* equality jacobian */
87: TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user); /* inequality jacobian */
88: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
89: TaoSetTolerances(tao,1.e-6,1.e-6,1.e-6);
90: TaoSetConstraintTolerances(tao,1.e-6,1.e-6);
92: TaoGetKSP(tao,&ksp);
93: KSPGetPC(ksp,&pc);
94: PCSetType(pc,PCLU);
95: /*
96: This algorithm produces matrices with zeros along the diagonal therefore we use
97: SuperLU_DIST which provides shift to the diagonal
98: */
99: PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST); /* requires superlu_dist to solve pdipm */
100: KSPSetType(ksp,KSPPREONLY);
101: KSPSetFromOptions(ksp);
103: TaoSetFromOptions(tao);
105: TaoSolve(tao);
106: TaoGetSolutionVector(tao,&x);
107: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
109: /* Free objects */
110: DestroyProblem(&user);
111: TaoDestroy(&tao);
112: PetscFinalize();
113: return ierr;
114: }
116: PetscErrorCode InitializeProblem(AppCtx *user)
117: {
119: PetscMPIInt size;
120: PetscMPIInt rank;
121: PetscInt nloc,neloc,niloc;
124: MPI_Comm_size(PETSC_COMM_WORLD,&size);
125: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
127: /* create vector x and set initial values */
128: user->n = 2; /* global length */
129: nloc = (rank==0)?user->n:0;
130: VecCreateMPI(PETSC_COMM_WORLD,nloc,user->n,&user->x);
131: VecSet(user->x,0);
133: /* create and set lower and upper bound vectors */
134: VecDuplicate(user->x,&user->xl);
135: VecDuplicate(user->x,&user->xu);
136: VecSet(user->xl,-1.0);
137: VecSet(user->xu,2.0);
139: /* create scater to zero */
140: VecScatterCreateToZero(user->x,&user->scat,&user->Xseq);
142: user->ne = 1;
143: neloc = (rank==0)?user->ne:0;
144: VecCreateMPI(PETSC_COMM_WORLD,neloc,user->ne,&user->ce); /* a 1x1 vec for equality constraints */
146: user->ni = 2;
147: niloc = (rank==0)?user->ni:0;
148: VecCreateMPI(PETSC_COMM_WORLD,niloc,user->ni,&user->ci); /* a 2x1 vec for inequality constraints */
150: /* nexn & nixn matricies for equaly and inequalty constriants */
151: MatCreate(PETSC_COMM_WORLD,&user->Ae);
152: MatCreate(PETSC_COMM_WORLD,&user->Ai);
153: MatCreate(PETSC_COMM_WORLD,&user->H);
155: MatSetSizes(user->Ae,neloc,nloc,user->ne,user->n);
156: MatSetSizes(user->Ai,niloc,nloc,user->ni,user->n);
157: MatSetSizes(user->H,nloc,nloc,user->n,user->n);
159: MatSetUp(user->Ae);
160: MatSetUp(user->Ai);
161: MatSetUp(user->H);
163: MatSetFromOptions(user->Ae);
164: MatSetFromOptions(user->Ai);
165: MatSetFromOptions(user->H);
166: return(0);
167: }
169: PetscErrorCode DestroyProblem(AppCtx *user)
170: {
174: MatDestroy(&user->Ae);
175: MatDestroy(&user->Ai);
176: MatDestroy(&user->H);
178: VecDestroy(&user->x);
179: VecDestroy(&user->ce);
180: VecDestroy(&user->ci);
181: VecDestroy(&user->xl);
182: VecDestroy(&user->xu);
183: VecDestroy(&user->Xseq);
184: VecScatterDestroy(&user->scat);
185: return(0);
186: }
188: /*
189: f(X) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
190: G(X) = fx = [2.0*(x[0]-2.0) - 2.0; 2.0*(x[1]-2.0) - 2.0]
191: */
192: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
193: {
194: PetscScalar g;
195: const PetscScalar *x;
196: MPI_Comm comm;
197: PetscMPIInt rank;
198: PetscErrorCode ierr;
199: PetscReal fin;
200: AppCtx *user=(AppCtx*)ctx;
201: Vec Xseq=user->Xseq;
202: VecScatter scat=user->scat;
205: PetscObjectGetComm((PetscObject)tao,&comm);
206: MPI_Comm_rank(comm,&rank);
208: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
209: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
211: fin = 0.0;
212: if (!rank) {
213: VecGetArrayRead(Xseq,&x);
214: fin = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
215: g = 2.0*(x[0]-2.0) - 2.0;
216: VecSetValue(G,0,g,INSERT_VALUES);
217: g = 2.0*(x[1]-2.0) - 2.0;
218: VecSetValue(G,1,g,INSERT_VALUES);
219: VecRestoreArrayRead(Xseq,&x);
220: }
221: MPI_Allreduce(&fin,f,1,MPIU_REAL,MPIU_SUM,comm);
222: VecAssemblyBegin(G);
223: VecAssemblyEnd(G);
224: return(0);
225: }
227: PetscErrorCode FormHessian(Tao tao, Vec x,Mat H, Mat Hpre, void *ctx)
228: {
229: AppCtx *user=(AppCtx*)ctx;
230: Vec DE,DI;
231: const PetscScalar *de, *di;
232: PetscInt zero=0,one=1;
233: PetscScalar two=2.0;
234: PetscScalar val=0.0;
235: Vec Deseq,Diseq=user->Xseq;
236: VecScatter Descat,Discat=user->scat;
237: PetscMPIInt rank;
238: MPI_Comm comm;
239: PetscErrorCode ierr;
242: TaoGetDualVariables(tao,&DE,&DI);
244: PetscObjectGetComm((PetscObject)tao,&comm);
245: MPI_Comm_rank(comm,&rank);
247: VecScatterCreateToZero(DE,&Descat,&Deseq);
249: VecScatterBegin(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
250: VecScatterBegin(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
251: VecScatterEnd(Descat,DE,Deseq,INSERT_VALUES,SCATTER_FORWARD);
252: VecScatterEnd(Discat,DI,Diseq,INSERT_VALUES,SCATTER_FORWARD);
254: if (!rank){
255: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
256: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
257: val = 2.0 * (1 + de[0] + di[0] - di[1]);
258: VecRestoreArrayRead(Deseq,&de);
259: VecRestoreArrayRead(Diseq,&di);
260: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
261: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
262: }
264: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
267: VecScatterDestroy(&Descat);
268: VecDestroy(&Deseq);
269: return(0);
270: }
272: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
273: {
274: const PetscScalar *x;
275: PetscScalar ci;
276: PetscErrorCode ierr;
277: MPI_Comm comm;
278: PetscMPIInt rank;
279: AppCtx *user=(AppCtx*)ctx;
280: Vec Xseq=user->Xseq;
281: VecScatter scat=user->scat;
284: PetscObjectGetComm((PetscObject)tao,&comm);
285: MPI_Comm_rank(comm,&rank);
287: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
288: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
290: if (!rank) {
291: VecGetArrayRead(Xseq,&x);
292: ci = x[0]*x[0] - x[1];
293: VecSetValue(CI,0,ci,INSERT_VALUES);
294: ci = -x[0]*x[0] + x[1] + 1.0;
295: VecSetValue(CI,1,ci,INSERT_VALUES);
296: VecRestoreArrayRead(Xseq,&x);
297: }
298: VecAssemblyBegin(CI);
299: VecAssemblyEnd(CI);
300: return(0);
301: }
303: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
304: {
305: const PetscScalar *x;
306: PetscScalar ce;
307: PetscErrorCode ierr;
308: MPI_Comm comm;
309: PetscMPIInt rank;
310: AppCtx *user=(AppCtx*)ctx;
311: Vec Xseq=user->Xseq;
312: VecScatter scat=user->scat;
315: PetscObjectGetComm((PetscObject)tao,&comm);
316: MPI_Comm_rank(comm,&rank);
318: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
319: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
321: if (!rank) {
322: VecGetArrayRead(Xseq,&x);
323: ce = x[0]*x[0] + x[1] - 2.0;
324: VecSetValue(CE,0,ce,INSERT_VALUES);
325: VecRestoreArrayRead(Xseq,&x);
326: }
327: VecAssemblyBegin(CE);
328: VecAssemblyEnd(CE);
329: return(0);
330: }
332: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
333: {
334: AppCtx *user=(AppCtx*)ctx;
335: PetscInt cols[2],min,max,i;
336: PetscScalar vals[2];
337: const PetscScalar *x;
338: PetscErrorCode ierr;
339: Vec Xseq=user->Xseq;
340: VecScatter scat=user->scat;
341: MPI_Comm comm;
342: PetscMPIInt rank;
345: PetscObjectGetComm((PetscObject)tao,&comm);
346: MPI_Comm_rank(comm,&rank);
347: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
348: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
350: VecGetArrayRead(Xseq,&x);
351: MatGetOwnershipRange(JI,&min,&max);
353: cols[0] = 0; cols[1] = 1;
354: for (i=min;i<max;i++) {
355: if (i==0){
356: vals[0] = +2*x[0]; vals[1] = -1.0;
357: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
358: }
359: if (i==1) {
360: vals[0] = -2*x[0]; vals[1] = +1.0;
361: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
362: }
363: }
364: VecRestoreArrayRead(Xseq,&x);
366: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
367: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
368: return(0);
369: }
371: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
372: {
373: PetscInt rows[2];
374: PetscScalar vals[2];
375: const PetscScalar *x;
376: PetscMPIInt rank;
377: MPI_Comm comm;
378: PetscErrorCode ierr;
381: PetscObjectGetComm((PetscObject)tao,&comm);
382: MPI_Comm_rank(comm,&rank);
384: if (!rank) {
385: VecGetArrayRead(X,&x);
386: rows[0] = 0; rows[1] = 1;
387: vals[0] = 2*x[0]; vals[1] = 1.0;
388: MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
389: VecRestoreArrayRead(X,&x);
390: }
391: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
392: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
393: return(0);
394: }
397: /*TEST
399: build:
400: requires: !complex !define(PETSC_USE_CXX)
402: test:
403: requires: superlu_dist
404: args: -tao_converged_reason
406: test:
407: suffix: 2
408: requires: superlu_dist
409: nsize: 2
410: args: -tao_converged_reason
412: TEST*/