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