1 /* Program usage: mpiexec -n 1 toy[-help] [all TAO options] */
2
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 ---------------------------------------------------------------------- */
9
10 #include <petsctao.h>
11
12 static char help[]="";
13
14 /*
15 User-defined application context - contains data needed by the
16 application-provided call-back routines, FormFunction(),
17 FormGradient(), and FormHessian().
18 */
19
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;
40
41 /* -------- User-defined Routines --------- */
42
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*);
51
52
53
main(int argc,char ** argv)54 PetscErrorCode main(int argc,char **argv)
55 {
56 PetscErrorCode ierr; /* used to check for functions returning nonzeros */
57 Tao tao;
58 KSP ksp;
59 PC pc;
60 AppCtx user; /* application context */
61
62 ierr = PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
63 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n---- TOY Problem -----\n");CHKERRQ(ierr);
64 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution should be f(1,1)=-2\n");CHKERRQ(ierr);
65 ierr = InitializeProblem(&user);CHKERRQ(ierr);
66 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
67 ierr = TaoSetType(tao,TAOIPM);CHKERRQ(ierr);
68 ierr = TaoSetInitialVector(tao,user.x);CHKERRQ(ierr);
69 ierr = TaoSetVariableBounds(tao,user.xl,user.xu);CHKERRQ(ierr);
70 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr);
71
72 ierr = TaoSetEqualityConstraintsRoutine(tao,user.ce,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr);
73 ierr = TaoSetInequalityConstraintsRoutine(tao,user.ci,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr);
74
75 ierr = TaoSetJacobianEqualityRoutine(tao,user.Ae,user.Ae,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr);
76 ierr = TaoSetJacobianInequalityRoutine(tao,user.Ai,user.Ai,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr);
77 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr);
78 /* ierr = TaoSetTolerances(tao,0,0,0);CHKERRQ(ierr); */
79
80 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
81 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
82 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
83 /*
84 This algorithm produces matrices with zeros along the diagonal therefore we need to use
85 SuperLU which does partial pivoting
86 */
87 ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr);
88 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
89 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
90
91 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
92 /* ierr = TaoSetTolerances(tao,0,0,0);CHKERRQ(ierr); */
93 ierr = TaoSolve(tao);CHKERRQ(ierr);
94
95 ierr = DestroyProblem(&user);CHKERRQ(ierr);
96 ierr = TaoDestroy(&tao);CHKERRQ(ierr);
97 ierr = PetscFinalize();
98 return ierr;
99 }
100
InitializeProblem(AppCtx * user)101 PetscErrorCode InitializeProblem(AppCtx *user)
102 {
103 PetscErrorCode ierr;
104
105 PetscFunctionBegin;
106 user->n = 2;
107 ierr = VecCreateSeq(PETSC_COMM_SELF,user->n,&user->x);CHKERRQ(ierr);
108 ierr = VecDuplicate(user->x,&user->xl);CHKERRQ(ierr);
109 ierr = VecDuplicate(user->x,&user->xu);CHKERRQ(ierr);
110 ierr = VecSet(user->x,0.0);CHKERRQ(ierr);
111 ierr = VecSet(user->xl,-1.0);CHKERRQ(ierr);
112 ierr = VecSet(user->xu,2.0);CHKERRQ(ierr);
113
114 user->ne = 1;
115 ierr = VecCreateSeq(PETSC_COMM_SELF,user->ne,&user->ce);CHKERRQ(ierr);
116
117 user->ni = 2;
118 ierr = VecCreateSeq(PETSC_COMM_SELF,user->ni,&user->ci);CHKERRQ(ierr);
119
120 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,user->ne,user->n,user->n,NULL,&user->Ae);CHKERRQ(ierr);
121 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,user->ni,user->n,user->n,NULL,&user->Ai);CHKERRQ(ierr);
122 ierr = MatSetFromOptions(user->Ae);CHKERRQ(ierr);
123 ierr = MatSetFromOptions(user->Ai);CHKERRQ(ierr);
124
125
126 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,user->n,user->n,1,NULL,&user->H);CHKERRQ(ierr);
127 ierr = MatSetFromOptions(user->H);CHKERRQ(ierr);
128
129 PetscFunctionReturn(0);
130 }
131
DestroyProblem(AppCtx * user)132 PetscErrorCode DestroyProblem(AppCtx *user)
133 {
134 PetscErrorCode ierr;
135
136 PetscFunctionBegin;
137 ierr = MatDestroy(&user->Ae);CHKERRQ(ierr);
138 ierr = MatDestroy(&user->Ai);CHKERRQ(ierr);
139 ierr = MatDestroy(&user->H);CHKERRQ(ierr);
140
141 ierr = VecDestroy(&user->x);CHKERRQ(ierr);
142 ierr = VecDestroy(&user->ce);CHKERRQ(ierr);
143 ierr = VecDestroy(&user->ci);CHKERRQ(ierr);
144 ierr = VecDestroy(&user->xl);CHKERRQ(ierr);
145 ierr = VecDestroy(&user->xu);CHKERRQ(ierr);
146 PetscFunctionReturn(0);
147 }
148
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,void * ctx)149 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
150 {
151 PetscScalar *g;
152 const PetscScalar *x;
153 PetscErrorCode ierr;
154
155 PetscFunctionBegin;
156 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
157 ierr = VecGetArray(G,&g);CHKERRQ(ierr);
158 *f = (x[0]-2.0)*(x[0]-2.0) + (x[1]-2.0)*(x[1]-2.0) - 2.0*(x[0]+x[1]);
159 g[0] = 2.0*(x[0]-2.0) - 2.0;
160 g[1] = 2.0*(x[1]-2.0) - 2.0;
161 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
162 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);
163 PetscFunctionReturn(0);
164 }
165
FormHessian(Tao tao,Vec x,Mat H,Mat Hpre,void * ctx)166 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
167 {
168 Vec DE,DI;
169 const PetscScalar *de, *di;
170 PetscInt zero=0,one=1;
171 PetscScalar two=2.0;
172 PetscScalar val;
173 PetscErrorCode ierr;
174
175 PetscFunctionBegin;
176 ierr = TaoGetDualVariables(tao,&DE,&DI);CHKERRQ(ierr);
177
178 ierr = VecGetArrayRead(DE,&de);CHKERRQ(ierr);
179 ierr = VecGetArrayRead(DI,&di);CHKERRQ(ierr);
180 val=2.0 * (1 + de[0] + di[0] - di[1]);
181 ierr = VecRestoreArrayRead(DE,&de);CHKERRQ(ierr);
182 ierr = VecRestoreArrayRead(DI,&di);CHKERRQ(ierr);
183
184 ierr = MatSetValues(H,1,&zero,1,&zero,&val,INSERT_VALUES);CHKERRQ(ierr);
185 ierr = MatSetValues(H,1,&one,1,&one,&two,INSERT_VALUES);CHKERRQ(ierr);
186
187 ierr = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
188 ierr = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
189 PetscFunctionReturn(0);
190 }
191
FormInequalityConstraints(Tao tao,Vec X,Vec CI,void * ctx)192 PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
193 {
194 const PetscScalar *x;
195 PetscScalar *c;
196 PetscErrorCode ierr;
197
198 PetscFunctionBegin;
199 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
200 ierr = VecGetArray(CI,&c);CHKERRQ(ierr);
201 c[0] = x[0]*x[0] - x[1];
202 c[1] = -x[0]*x[0] + x[1] + 1.0;
203 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
204 ierr = VecRestoreArray(CI,&c);CHKERRQ(ierr);
205 PetscFunctionReturn(0);
206 }
207
FormEqualityConstraints(Tao tao,Vec X,Vec CE,void * ctx)208 PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE,void *ctx)
209 {
210 PetscScalar *x,*c;
211 PetscErrorCode ierr;
212
213 PetscFunctionBegin;
214 ierr = VecGetArray(X,&x);CHKERRQ(ierr);
215 ierr = VecGetArray(CE,&c);CHKERRQ(ierr);
216 c[0] = x[0]*x[0] + x[1] - 2.0;
217 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
218 ierr = VecRestoreArray(CE,&c);CHKERRQ(ierr);
219 PetscFunctionReturn(0);
220 }
221
FormInequalityJacobian(Tao tao,Vec X,Mat JI,Mat JIpre,void * ctx)222 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
223 {
224 PetscInt rows[2];
225 PetscInt cols[2];
226 PetscScalar vals[4];
227 const PetscScalar *x;
228 PetscErrorCode ierr;
229
230 PetscFunctionBegin;
231 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
232 rows[0] = 0; rows[1] = 1;
233 cols[0] = 0; cols[1] = 1;
234 vals[0] = +2*x[0]; vals[1] = -1.0;
235 vals[2] = -2*x[0]; vals[3] = +1.0;
236 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
237 ierr = MatSetValues(JI,2,rows,2,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
238 ierr = MatAssemblyBegin(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
239 ierr = MatAssemblyEnd(JI,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
240
241 PetscFunctionReturn(0);
242 }
243
FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,void * ctx)244 PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
245 {
246 PetscInt rows[2];
247 PetscScalar vals[2];
248 const PetscScalar *x;
249 PetscErrorCode ierr;
250
251 PetscFunctionBegin;
252 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
253 rows[0] = 0; rows[1] = 1;
254 vals[0] = 2*x[0]; vals[1] = 1.0;
255 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
256 ierr = MatSetValues(JE,1,rows,2,rows,vals,INSERT_VALUES);CHKERRQ(ierr);
257 ierr = MatAssemblyBegin(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
258 ierr = MatAssemblyEnd(JE,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
259 PetscFunctionReturn(0);
260 }
261
262
263 /*TEST
264
265 build:
266 requires: !complex !define(PETSC_USE_CXX)
267
268 test:
269 requires: superlu
270 args: -tao_smonitor -tao_view -tao_gatol 1.e-5
271
272 TEST*/
273