Actual source code: toy.c
petsc-3.6.1 2015-08-06
1: /* Program usage: mpirun -np 1 toy[-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
4: min f=(x1-x2)^2 + (x2-2)^2 -2*x1-2*x2
5: s.t. x1^2 + x2 = 2
6: 0 <= x1^2 - x2 <= 1
7: -1 <= x1,x2 <= 2
8: ---------------------------------------------------------------------- */
10: #include <petsctao.h>
12: static char help[]="";
14: /*
15: User-defined application context - contains data needed by the
16: application-provided call-back routines, FormFunction(),
17: FormGradient(), and FormHessian().
18: */
20: /*
21: x,d in R^n
22: f in R
23: bin in R^mi
24: beq in R^me
25: Aeq in R^(me x n)
26: Ain in R^(mi x n)
27: H in R^(n x n)
28: min f=(1/2)*x'*H*x + d'*x
29: s.t. Aeq*x == beq
30: Ain*x >= bin
31: */
32: typedef struct {
33: PetscInt n; /* Length x */
34: PetscInt ne; /* number of equality constraints */
35: PetscInt ni; /* number of inequality constraints */
36: Vec x,xl,xu;
37: Vec ce,ci,bl,bu;
38: Mat Ae,Ai,H;
39: } AppCtx;
41: /* -------- User-defined Routines --------- */
43: PetscErrorCode InitializeProblem(AppCtx *);
44: PetscErrorCode DestroyProblem(AppCtx *);
45: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
46: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
47: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
48: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
49: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
50: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
56: PetscErrorCode main(int argc,char **argv)
57: {
58: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
59: Tao tao;
60: TaoConvergedReason reason;
61: KSP ksp;
62: PC pc;
63: AppCtx user; /* application context */
65: PetscInitialize(&argc,&argv,(char *)0,help);
66: PetscPrintf(PETSC_COMM_WORLD,"\n---- TOY Problem -----\n");
67: PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");
68: InitializeProblem(&user);
69: TaoCreate(PETSC_COMM_WORLD,&tao);
70: TaoSetType(tao,TAOIPM);
71: TaoSetInitialVector(tao,user.x);
72: TaoSetVariableBounds(tao,user.xl,user.xu);
73: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
75: TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);
76: TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);
78: TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);
79: TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);
80: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
81: TaoSetTolerances(tao,1e-12,0,0,0,0);
83: TaoSetFromOptions(tao);
85: TaoGetKSP(tao,&ksp);
86: KSPGetPC(ksp,&pc);
87: PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);
88: /* TODO -- why didn't that work? */
89: PetscOptionsSetValue("-pc_factor_mat_solver_package","superlu");
90: PCSetType(pc,PCLU);
91: KSPSetType(ksp,KSPPREONLY);
92: KSPSetFromOptions(ksp);
94: TaoSetTolerances(tao,1e-12,0,0,0,0);
96: TaoSolve(tao);
98: /* Analyze solution */
99: TaoGetConvergedReason(tao,&reason);
100: if (reason < 0) {
101: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
102: } else {
103: PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
104: }
106: DestroyProblem(&user);
107: TaoDestroy(&tao);
108: PetscFinalize();
109: return 0;
110: }
114: PetscErrorCode InitializeProblem(AppCtx *user)
115: {
119: user->n = 2;
120: VecCreateSeq(PETSC_COMM_SELF,user->n,&user->x);
121: VecDuplicate(user->x,&user->xl);
122: VecDuplicate(user->x,&user->xu);
123: VecSet(user->x,0.0);
124: VecSet(user->xl,-1.0);
125: VecSet(user->xu,2.0);
127: user->ne = 1;
128: VecCreateSeq(PETSC_COMM_SELF,user->ne,&user->ce);
130: user->ni = 2;
131: VecCreateSeq(PETSC_COMM_SELF,user->ni,&user->ci);
133: MatCreateSeqAIJ(PETSC_COMM_SELF,user->ne,user->n,user->n,NULL,&user->Ae);
134: MatCreateSeqAIJ(PETSC_COMM_SELF,user->ni,user->n,user->n,NULL,&user->Ai);
135: MatSetFromOptions(user->Ae);
136: MatSetFromOptions(user->Ai);
139: MatCreateSeqAIJ(PETSC_COMM_SELF,user->n,user->n,1,NULL,&user->H); MatSetFromOptions(user->H);
141: return(0);
142: }
146: PetscErrorCode DestroyProblem(AppCtx *user)
147: {
151: MatDestroy(&user->Ae);
152: MatDestroy(&user->Ai);
153: MatDestroy(&user->H);
155: VecDestroy(&user->x);
156: VecDestroy(&user->ce);
157: VecDestroy(&user->ci);
158: VecDestroy(&user->xl);
159: VecDestroy(&user->xu);
160: return(0);
161: }
165: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
166: {
167: PetscScalar *g;
168: const PetscScalar *x;
169: PetscErrorCode ierr;
172: VecGetArrayRead(X,&x);
173: VecGetArray(G,&g);
174: *f = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
175: g[0] = 2.0*(x[0]-2.0) - 2.0;
176: g[1] = 2.0*(x[1]-2.0) - 2.0;
177: VecRestoreArrayRead(X,&x);
178: VecRestoreArray(G,&g);
179: return(0);
180: }
184: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
185: {
186: Vec DE,DI;
187: const PetscScalar *de, *di;
188: PetscInt zero=0,one=1;
189: PetscScalar two=2.0;
190: PetscScalar val;
191: PetscErrorCode ierr;
194: TaoGetDualVariables(tao,&DE,&DI);
196: VecGetArrayRead(DE,&de);
197: VecGetArrayRead(DI,&di);
198: val=2.0 * (1 + de[0] + di[0] - di[1]);
199: VecRestoreArrayRead(DE,&de);
200: VecRestoreArrayRead(DI,&di);
202: MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);
203: MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);
205: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
206: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
207: return(0);
208: }
212: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
213: {
214: const PetscScalar *x;
215: PetscScalar *c;
216: PetscErrorCode ierr;
219: VecGetArrayRead(X,&x);
220: VecGetArray(CI,&c);
221: c[0] = x[0]*x[0] - x[1];
222: c[1] = -x[0]*x[0] + x[1] + 1.0;
223: VecRestoreArrayRead(X,&x);
224: VecRestoreArray(CI,&c);
225: return(0);
226: }
230: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE,void *ctx)
231: {
232: PetscScalar *x,*c;
236: VecGetArray(X,&x);
237: VecGetArray(CE,&c);
238: c[0] = x[0]*x[0] + x[1] - 2.0;
239: VecRestoreArray(X,&x);
240: VecRestoreArray(CE,&c);
241: return(0);
242: }
246: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
247: {
248: PetscInt rows[2];
249: PetscInt cols[2];
250: PetscScalar vals[4];
251: const PetscScalar *x;
252: PetscErrorCode ierr;
255: VecGetArrayRead(X,&x);
256: rows[0] = 0; rows[1] = 1;
257: cols[0] = 0; cols[1] = 1;
258: vals[0] = +2*x[0]; vals[1] = -1.0;
259: vals[2] = -2*x[0]; vals[3] = +1.0;
260: VecRestoreArrayRead(X,&x);
261: MatSetValues(JI,2,rows,2,cols,vals,INSERT_VALUES);
262: MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);
263: MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);
265: return(0);
266: }
270: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
271: {
272: PetscInt rows[2];
273: PetscScalar vals[2];
274: const PetscScalar *x;
275: PetscErrorCode ierr;
278: VecGetArrayRead(X,&x);
279: rows[0] = 0; rows[1] = 1;
280: vals[0] = 2*x[0]; vals[1] = 1.0;
281: VecRestoreArrayRead(X,&x);
282: MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);
283: MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);
284: MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);
285: return(0);
286: }