1 /* Program usage: mpiexec -n 1 maros1 [-help] [all TAO options] */
2
3 /* ----------------------------------------------------------------------
4 TODO Explain maros example
5 ---------------------------------------------------------------------- */
6
7 #include <petsctao.h>
8
9 static char help[]="";
10
11 /*T
12 Concepts: TAO^Solving an unconstrained minimization problem
13 Routines: TaoCreate(); TaoSetType();
14 Routines: TaoSetInitialVector();
15 Routines: TaoSetObjectiveAndGradientRoutine();
16 Routines: TaoSetEqualityConstraintsRoutine();
17 Routines: TaoSetInequalityConstraintsRoutine();
18 Routines: TaoSetEqualityJacobianRoutine();
19 Routines: TaoSetInequalityJacobianRoutine();
20 Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
21 Routines: TaoGetKSP(); TaoSolve();
22 Routines: TaoGetConvergedReason();TaoDestroy();
23 Processors: 1
24 T*/
25
26 /*
27 User-defined application context - contains data needed by the
28 application-provided call-back routines, FormFunction(),
29 FormGradient(), and FormHessian().
30 */
31
32 /*
33 x,d in R^n
34 f in R
35 bin in R^mi
36 beq in R^me
37 Aeq in R^(me x n)
38 Ain in R^(mi x n)
39 H in R^(n x n)
40 min f=(1/2)*x'*H*x + d'*x
41 s.t. Aeq*x == beq
42 Ain*x >= bin
43 */
44 typedef struct {
45 char name[32];
46 PetscInt n; /* Length x */
47 PetscInt me; /* number of equality constraints */
48 PetscInt mi; /* number of inequality constraints */
49 PetscInt m; /* me+mi */
50 Mat Aeq,Ain,H;
51 Vec beq,bin,d;
52 } AppCtx;
53
54 /* -------- User-defined Routines --------- */
55
56 PetscErrorCode InitializeProblem(AppCtx*);
57 PetscErrorCode DestroyProblem(AppCtx *);
58 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
59 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
60 PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
61 PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
62 PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
63 PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);
64
main(int argc,char ** argv)65 PetscErrorCode main(int argc,char **argv)
66 {
67 PetscErrorCode ierr; /* used to check for functions returning nonzeros */
68 PetscMPIInt size;
69 Vec x; /* solution */
70 KSP ksp;
71 PC pc;
72 Vec ceq,cin;
73 PetscBool flg; /* A return value when checking for use options */
74 Tao tao; /* Tao solver context */
75 TaoConvergedReason reason;
76 AppCtx user; /* application context */
77
78 /* Initialize TAO,PETSc */
79 ierr = PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
80 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
81 /* Specify default parameters for the problem, check for command-line overrides */
82 ierr = PetscStrncpy(user.name,"HS21",sizeof(user.name));CHKERRQ(ierr);
83 ierr = PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg);CHKERRQ(ierr);
84
85 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);CHKERRQ(ierr);
86 ierr = InitializeProblem(&user);CHKERRQ(ierr);
87 ierr = VecDuplicate(user.d,&x);CHKERRQ(ierr);
88 ierr = VecDuplicate(user.beq,&ceq);CHKERRQ(ierr);
89 ierr = VecDuplicate(user.bin,&cin);CHKERRQ(ierr);
90 ierr = VecSet(x,1.0);CHKERRQ(ierr);
91
92 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
93 ierr = TaoSetType(tao,TAOIPM);CHKERRQ(ierr);
94 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
95 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr);
96 ierr = TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr);
97 ierr = TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr);
98 ierr = TaoSetInequalityBounds(tao,user.bin,NULL);CHKERRQ(ierr);
99 ierr = TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr);
100 ierr = TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr);
101 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr);
102 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
103 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
104 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
105 /*
106 This algorithm produces matrices with zeros along the diagonal therefore we need to use
107 SuperLU which does partial pivoting
108 */
109 ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr);
110 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
111 ierr = TaoSetTolerances(tao,0,0,0);CHKERRQ(ierr);
112
113 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
114 ierr = TaoSolve(tao);CHKERRQ(ierr);
115 ierr = TaoGetConvergedReason(tao,&reason);CHKERRQ(ierr);
116 if (reason < 0) {
117 ierr = PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);CHKERRQ(ierr);
118 } else {
119 ierr = PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);CHKERRQ(ierr);
120 }
121
122 ierr = DestroyProblem(&user);CHKERRQ(ierr);
123 ierr = VecDestroy(&x);CHKERRQ(ierr);
124 ierr = VecDestroy(&ceq);CHKERRQ(ierr);
125 ierr = VecDestroy(&cin);CHKERRQ(ierr);
126 ierr = TaoDestroy(&tao);CHKERRQ(ierr);
127
128 ierr = PetscFinalize();
129 return ierr;
130 }
131
InitializeProblem(AppCtx * user)132 PetscErrorCode InitializeProblem(AppCtx *user)
133 {
134 PetscErrorCode ierr;
135 PetscViewer loader;
136 MPI_Comm comm;
137 PetscInt nrows,ncols,i;
138 PetscScalar one=1.0;
139 char filebase[128];
140 char filename[128];
141
142 PetscFunctionBegin;
143 comm = PETSC_COMM_WORLD;
144 ierr = PetscStrncpy(filebase,user->name,sizeof(filebase));CHKERRQ(ierr);
145 ierr = PetscStrlcat(filebase,"/",sizeof(filebase));CHKERRQ(ierr);
146 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr);
147 ierr = PetscStrlcat(filename,"f",sizeof(filename));CHKERRQ(ierr);
148 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr);
149
150 ierr = VecCreate(comm,&user->d);CHKERRQ(ierr);
151 ierr = VecLoad(user->d,loader);CHKERRQ(ierr);
152 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr);
153 ierr = VecGetSize(user->d,&nrows);CHKERRQ(ierr);
154 ierr = VecSetFromOptions(user->d);CHKERRQ(ierr);
155 user->n = nrows;
156
157 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr);
158 ierr = PetscStrlcat(filename,"H",sizeof(filename));CHKERRQ(ierr);
159 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr);
160
161 ierr = MatCreate(comm,&user->H);CHKERRQ(ierr);
162 ierr = MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);CHKERRQ(ierr);
163 ierr = MatLoad(user->H,loader);CHKERRQ(ierr);
164 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr);
165 ierr = MatGetSize(user->H,&nrows,&ncols);CHKERRQ(ierr);
166 if (nrows != user->n) SETERRQ(comm,0,"H: nrows != n\n");
167 if (ncols != user->n) SETERRQ(comm,0,"H: ncols != n\n");
168 ierr = MatSetFromOptions(user->H);CHKERRQ(ierr);
169
170 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr);
171 ierr = PetscStrlcat(filename,"Aeq",sizeof(filename));CHKERRQ(ierr);
172 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
173 if (ierr) {
174 user->Aeq = NULL;
175 user->me = 0;
176 } else {
177 ierr = MatCreate(comm,&user->Aeq);CHKERRQ(ierr);
178 ierr = MatLoad(user->Aeq,loader);CHKERRQ(ierr);
179 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr);
180 ierr = MatGetSize(user->Aeq,&nrows,&ncols);CHKERRQ(ierr);
181 if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n");
182 ierr = MatSetFromOptions(user->Aeq);CHKERRQ(ierr);
183 user->me = nrows;
184 }
185
186 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr);
187 ierr = PetscStrlcat(filename,"Beq",sizeof(filename));CHKERRQ(ierr);
188 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr);
189 if (ierr) {
190 user->beq = 0;
191 } else {
192 ierr = VecCreate(comm,&user->beq);CHKERRQ(ierr);
193 ierr = VecLoad(user->beq,loader);CHKERRQ(ierr);
194 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr);
195 ierr = VecGetSize(user->beq,&nrows);CHKERRQ(ierr);
196 if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n");
197 ierr = VecSetFromOptions(user->beq);CHKERRQ(ierr);
198 }
199
200 user->mi = user->n;
201 /* Ain = eye(n,n) */
202 ierr = MatCreate(comm,&user->Ain);CHKERRQ(ierr);
203 ierr = MatSetType(user->Ain,MATAIJ);CHKERRQ(ierr);
204 ierr = MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);CHKERRQ(ierr);
205
206 ierr = MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);CHKERRQ(ierr);
207 ierr = MatSeqAIJSetPreallocation(user->Ain,1,NULL);CHKERRQ(ierr);
208
209 for (i=0;i<user->mi;i++) {
210 ierr = MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr);
211 }
212 ierr = MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
213 ierr = MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
214 ierr = MatSetFromOptions(user->Ain);CHKERRQ(ierr);
215
216 /* bin = [0,0 ... 0]' */
217 ierr = VecCreate(comm,&user->bin);CHKERRQ(ierr);
218 ierr = VecSetType(user->bin,VECMPI);CHKERRQ(ierr);
219 ierr = VecSetSizes(user->bin,PETSC_DECIDE,user->mi);CHKERRQ(ierr);
220 ierr = VecSet(user->bin,0.0);CHKERRQ(ierr);
221 ierr = VecSetFromOptions(user->bin);CHKERRQ(ierr);
222 user->m = user->me + user->mi;
223 PetscFunctionReturn(0);
224 }
225
DestroyProblem(AppCtx * user)226 PetscErrorCode DestroyProblem(AppCtx *user)
227 {
228 PetscErrorCode ierr;
229
230 PetscFunctionBegin;
231 ierr = MatDestroy(&user->H);CHKERRQ(ierr);
232 ierr = MatDestroy(&user->Aeq);CHKERRQ(ierr);
233 ierr = MatDestroy(&user->Ain);CHKERRQ(ierr);
234 ierr = VecDestroy(&user->beq);CHKERRQ(ierr);
235 ierr = VecDestroy(&user->bin);CHKERRQ(ierr);
236 ierr = VecDestroy(&user->d);CHKERRQ(ierr);
237 PetscFunctionReturn(0);
238 }
FormFunctionGradient(Tao tao,Vec x,PetscReal * f,Vec g,void * ctx)239 PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
240 {
241 AppCtx *user = (AppCtx*)ctx;
242 PetscScalar xtHx;
243 PetscErrorCode ierr;
244
245 PetscFunctionBegin;
246 ierr = MatMult(user->H,x,g);CHKERRQ(ierr);
247 ierr = VecDot(x,g,&xtHx);CHKERRQ(ierr);
248 ierr = VecDot(x,user->d,f);CHKERRQ(ierr);
249 *f += 0.5*xtHx;
250 ierr = VecAXPY(g,1.0,user->d);CHKERRQ(ierr);
251 PetscFunctionReturn(0);
252 }
253
FormHessian(Tao tao,Vec x,Mat H,Mat Hpre,void * ctx)254 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
255 {
256 PetscFunctionBegin;
257 PetscFunctionReturn(0);
258 }
259
FormInequalityConstraints(Tao tao,Vec x,Vec ci,void * ctx)260 PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
261 {
262 AppCtx *user = (AppCtx*)ctx;
263 PetscErrorCode ierr;
264
265 PetscFunctionBegin;
266 ierr = MatMult(user->Ain,x,ci);CHKERRQ(ierr);
267 PetscFunctionReturn(0);
268 }
269
FormEqualityConstraints(Tao tao,Vec x,Vec ce,void * ctx)270 PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
271 {
272 AppCtx *user = (AppCtx*)ctx;
273 PetscErrorCode ierr;
274
275 PetscFunctionBegin;
276 ierr = MatMult(user->Aeq,x,ce);CHKERRQ(ierr);
277 ierr = VecAXPY(ce,-1.0,user->beq);CHKERRQ(ierr);
278 PetscFunctionReturn(0);
279 }
280
FormInequalityJacobian(Tao tao,Vec x,Mat JI,Mat JIpre,void * ctx)281 PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
282 {
283 PetscFunctionBegin;
284 PetscFunctionReturn(0);
285 }
286
FormEqualityJacobian(Tao tao,Vec x,Mat JE,Mat JEpre,void * ctx)287 PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
288 {
289 PetscFunctionBegin;
290 PetscFunctionReturn(0);
291 }
292
293
294 /*TEST
295
296 build:
297 requires: !complex
298
299 test:
300 requires: superlu
301 localrunfiles: HS21
302
303 TEST*/
304