Actual source code: ex1.c
petsc-3.13.6 2020-09-29
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 Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
23: Section 1.5 Writing Application Codes with PETSc-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; /* Section 1.5 Writing Application Codes with PETSc 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);
172:
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);
274:
276: if (!rank){
277: if (!user->noeqflag){
278: VecGetArrayRead(Deseq,&de); /* places equality constraint dual into array */
279: }
281: VecGetArrayRead(Diseq,&di); /* places inequality constraint dual into array */
282:
283: if (!user->noeqflag){
284: val = 2.0 * (1 + de[0] + di[0] - di[1]);
285: VecRestoreArrayRead(Deseq,&de);
286: }else{
287: val = 2.0 * (1 + di[0] - di[1]);
288: }
289: VecRestoreArrayRead(Diseq,&di);
290: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
291: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
292: }
294: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
296: if (!user->noeqflag){
297: VecScatterDestroy(&Descat);
298: VecDestroy(&Deseq);
299: }
300: return(0);
301: }
303: PetscErrorCode FormInequalityConstraints(Tao tao,Vec X,Vec CI,void *ctx)
304: {
305: const PetscScalar *x;
306: PetscScalar ci;
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: ci = x[0]*x[0] - x[1];
324: VecSetValue(CI,0,ci,INSERT_VALUES);
325: ci = -x[0]*x[0] + x[1] + 1.0;
326: VecSetValue(CI,1,ci,INSERT_VALUES);
327: VecRestoreArrayRead(Xseq,&x);
328: }
329: VecAssemblyBegin(CI);
330: VecAssemblyEnd(CI);
331: return(0);
332: }
334: PetscErrorCode FormEqualityConstraints(Tao tao,Vec X,Vec CE,void *ctx)
335: {
336: const PetscScalar *x;
337: PetscScalar ce;
338: PetscErrorCode ierr;
339: MPI_Comm comm;
340: PetscMPIInt rank;
341: AppCtx *user=(AppCtx*)ctx;
342: Vec Xseq=user->Xseq;
343: VecScatter scat=user->scat;
346: PetscObjectGetComm((PetscObject)tao,&comm);
347: MPI_Comm_rank(comm,&rank);
349: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
350: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
352: if (!rank) {
353: VecGetArrayRead(Xseq,&x);
354: ce = x[0]*x[0] + x[1] - 2.0;
355: VecSetValue(CE,0,ce,INSERT_VALUES);
356: VecRestoreArrayRead(Xseq,&x);
357: }
358: VecAssemblyBegin(CE);
359: VecAssemblyEnd(CE);
360: return(0);
361: }
363: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
364: {
365: AppCtx *user=(AppCtx*)ctx;
366: PetscInt cols[2],min,max,i;
367: PetscScalar vals[2];
368: const PetscScalar *x;
369: PetscErrorCode ierr;
370: Vec Xseq=user->Xseq;
371: VecScatter scat=user->scat;
372: MPI_Comm comm;
373: PetscMPIInt rank;
376: PetscObjectGetComm((PetscObject)tao,&comm);
377: MPI_Comm_rank(comm,&rank);
378: VecScatterBegin(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
379: VecScatterEnd(scat,X,Xseq,INSERT_VALUES,SCATTER_FORWARD);
381: VecGetArrayRead(Xseq,&x);
382: MatGetOwnershipRange(JI,&min,&max);
384: cols[0] = 0; cols[1] = 1;
385: for (i=min;i<max;i++) {
386: if (i==0){
387: vals[0] = +2*x[0]; vals[1] = -1.0;
388: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
389: }
390: if (i==1) {
391: vals[0] = -2*x[0]; vals[1] = +1.0;
392: MatSetValues(JI,1,&i,2,cols,vals,INSERT_VALUES);
393: }
394: }
395: VecRestoreArrayRead(Xseq,&x);
397: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
398: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
399: return(0);
400: }
402: PetscErrorCode FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void *ctx)
403: {
404: PetscInt rows[2];
405: PetscScalar vals[2];
406: const PetscScalar *x;
407: PetscMPIInt rank;
408: MPI_Comm comm;
409: PetscErrorCode ierr;
412: PetscObjectGetComm((PetscObject)tao,&comm);
413: MPI_Comm_rank(comm,&rank);
415: if (!rank) {
416: VecGetArrayRead(X,&x);
417: rows[0] = 0; rows[1] = 1;
418: vals[0] = 2*x[0]; vals[1] = 1.0;
419: MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
420: VecRestoreArrayRead(X,&x);
421: }
422: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
423: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
424: return(0);
425: }
428: /*TEST
430: build:
431: requires: !complex !define(PETSC_USE_CXX)
433: test:
434: requires: superlu_dist
435: args: -tao_converged_reason
437: test:
438: suffix: 2
439: requires: superlu_dist
440: nsize: 2
441: args: -tao_converged_reason
443: test:
444: suffix: 3
445: requires: superlu_dist
446: args: -tao_converged_reason -no_eq
448: test:
449: suffix: 4
450: requires: superlu_dist
451: nsize: 2
452: args: -tao_converged_reason -no_eq
454: TEST*/