1 /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */
2
3 /* ----------------------------------------------------------------------
4 min f(x) = (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 -->
9 g(x) = 0
10 h(x) >= 0
11 -1 <= x0, x1 <= 2
12 where
13 g(x) = x0^2 + x1 - 2
14 h(x) = [x0^2 - x1
15 1 -(x0^2 - x1)]
16 ---------------------------------------------------------------------- */
17
18 #include <petsctao.h>
19
20 static char help[] = "Solves constrained optimization problem using pdipm.\n\
21 Input parameters include:\n\
22 -tao_type pdipm : sets Tao solver\n\
23 -no_eq : removes the equaility constraints from the problem\n\
24 -no_bound : removes the bound constraints from the problem\n\
25 -init_view : view initial object setup\n\
26 -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
27 -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
28 -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\
29 -tao_monitor_solution : view exact solution at each iteration\n\
30 Note: external package MUMPS is required to run pdipm in parallel. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";
31
32 /*
33 User-defined application context - contains data needed by the application
34 */
35 typedef struct {
36 PetscInt n; /* Global length of x */
37 PetscInt ne; /* Global number of equality constraints */
38 PetscInt ni; /* Global number of inequality constraints */
39 PetscBool noeqflag, noboundflag, initview;
40 Vec x, xl, xu;
41 Vec ce, ci, bl, bu, Xseq;
42 Mat Ae, Ai, H;
43 VecScatter scat;
44 } AppCtx;
45
46 /* -------- User-defined Routines --------- */
47 PetscErrorCode InitializeProblem(AppCtx *);
48 PetscErrorCode DestroyProblem(AppCtx *);
49 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
50 PetscErrorCode FormPDIPMHessian(Tao, Vec, Mat, Mat, void *);
51 PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
52 PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
53 PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
54 PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);
55
main(int argc,char ** argv)56 int main(int argc, char **argv)
57 {
58 Tao tao;
59 KSP ksp;
60 PC pc;
61 AppCtx user; /* application context */
62 Vec x, G, CI, CE;
63 PetscMPIInt size;
64 PetscReal f;
65 PetscBool pdipm;
66
67 PetscFunctionBeginUser;
68 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
69 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
70 PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processes detected. Example written to use at most 2 processes.");
71
72 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n"));
73 PetscCall(InitializeProblem(&user)); /* sets up problem, function below */
74
75 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
76 PetscCall(TaoSetType(tao, TAOALMM));
77 PetscCall(TaoSetSolution(tao, user.x));
78 if (!user.noboundflag) PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu));
79 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
80 PetscCall(TaoSetTolerances(tao, 1.e-4, 0.0, 0.0));
81 PetscCall(TaoSetConstraintTolerances(tao, 1.e-4, 0.0));
82 PetscCall(TaoSetMaximumFunctionEvaluations(tao, 1e6));
83 PetscCall(TaoSetFromOptions(tao));
84
85 if (!user.noeqflag) {
86 PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user));
87 PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user));
88 }
89 PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user));
90 PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user));
91
92 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm));
93 if (pdipm) {
94 /*
95 PDIPM produces an indefinite KKT matrix with some zeros along the diagonal
96 Inverting this indefinite matrix requires PETSc to be configured with MUMPS
97 */
98 PetscCheck(PetscDefined(HAVE_MUMPS), PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)");
99 PetscCall(TaoGetKSP(tao, &ksp));
100 PetscCall(KSPGetPC(ksp, &pc));
101 PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
102 PetscCall(KSPSetType(ksp, KSPPREONLY));
103 PetscCall(PCSetType(pc, PCCHOLESKY));
104 PetscCall(KSPSetFromOptions(ksp));
105 PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user));
106 }
107
108 /* Print out an initial view of the problem */
109 if (user.initview) {
110 PetscCall(TaoSetUp(tao));
111 PetscCall(VecDuplicate(user.x, &G));
112 PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
113 if (pdipm) PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user));
114 PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
115 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
116 PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
119 PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
120 if (pdipm) PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
121 PetscCall(VecDestroy(&G));
122 PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
123 PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
124 PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
126 PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
127 PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
128 PetscCall(VecDestroy(&CI));
129 if (!user.noeqflag) {
130 PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
131 PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
132 PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
133 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
134 PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
135 PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
136 PetscCall(VecDestroy(&CE));
137 }
138 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
139 PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
140 }
141
142 PetscCall(TaoSolve(tao));
143 PetscCall(TaoGetSolution(tao, &x));
144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n"));
145 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
146
147 /* Free objects */
148 PetscCall(DestroyProblem(&user));
149 PetscCall(TaoDestroy(&tao));
150 PetscCall(PetscFinalize());
151 return PETSC_SUCCESS;
152 }
153
InitializeProblem(AppCtx * user)154 PetscErrorCode InitializeProblem(AppCtx *user)
155 {
156 PetscMPIInt size;
157 PetscMPIInt rank;
158 PetscInt nloc, neloc, niloc;
159
160 PetscFunctionBegin;
161 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
162 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
163 user->noeqflag = PETSC_FALSE;
164 user->noboundflag = PETSC_FALSE;
165 PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
166 PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_bound", &user->noboundflag, NULL));
167 user->initview = PETSC_FALSE;
168 PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));
169
170 /* Tell user the correct solution, not an error checking */
171 if (!user->noeqflag) {
172 /* With equality constraint, bound or no bound makes no different to the solution */
173 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
174 } else {
175 if (user->noboundflag) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq, -no_bound): f(2.05655, 3.22938) = -9.05728\n"));
176 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
177 }
178
179 /* create vector x and set initial values */
180 user->n = 2; /* global length */
181 nloc = (size == 1) ? user->n : 1;
182 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
183 PetscCall(VecSetSizes(user->x, nloc, user->n));
184 PetscCall(VecSetFromOptions(user->x));
185 PetscCall(VecSet(user->x, 0.0));
186
187 /* create and set lower and upper bound vectors */
188 PetscCall(VecDuplicate(user->x, &user->xl));
189 PetscCall(VecDuplicate(user->x, &user->xu));
190 PetscCall(VecSet(user->xl, -1.0));
191 PetscCall(VecSet(user->xu, 2.0));
192
193 /* create scater to zero */
194 PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));
195
196 user->ne = 1;
197 user->ni = 2;
198 neloc = (rank == 0) ? user->ne : 0;
199 niloc = (size == 1) ? user->ni : 1;
200
201 if (!user->noeqflag) {
202 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
203 PetscCall(VecSetSizes(user->ce, neloc, user->ne));
204 PetscCall(VecSetFromOptions(user->ce));
205 PetscCall(VecSetUp(user->ce));
206 }
207
208 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
209 PetscCall(VecSetSizes(user->ci, niloc, user->ni));
210 PetscCall(VecSetFromOptions(user->ci));
211 PetscCall(VecSetUp(user->ci));
212
213 /* nexn & nixn matrices for equally and inequalty constraints */
214 if (!user->noeqflag) {
215 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
216 PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
217 PetscCall(MatSetFromOptions(user->Ae));
218 PetscCall(MatSetUp(user->Ae));
219 }
220
221 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
222 PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
223 PetscCall(MatSetFromOptions(user->Ai));
224 PetscCall(MatSetUp(user->Ai));
225
226 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
227 PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
228 PetscCall(MatSetFromOptions(user->H));
229 PetscCall(MatSetUp(user->H));
230 PetscFunctionReturn(PETSC_SUCCESS);
231 }
232
DestroyProblem(AppCtx * user)233 PetscErrorCode DestroyProblem(AppCtx *user)
234 {
235 PetscFunctionBegin;
236 if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
237 PetscCall(MatDestroy(&user->Ai));
238 PetscCall(MatDestroy(&user->H));
239
240 PetscCall(VecDestroy(&user->x));
241 if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
242 PetscCall(VecDestroy(&user->ci));
243 PetscCall(VecDestroy(&user->xl));
244 PetscCall(VecDestroy(&user->xu));
245 PetscCall(VecDestroy(&user->Xseq));
246 PetscCall(VecScatterDestroy(&user->scat));
247 PetscFunctionReturn(PETSC_SUCCESS);
248 }
249
250 /* Evaluate
251 f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
252 G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6
253 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * f,Vec G,PetscCtx ctx)254 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, PetscCtx ctx)
255 {
256 const PetscScalar *x;
257 MPI_Comm comm;
258 PetscMPIInt rank;
259 PetscReal fin;
260 AppCtx *user = (AppCtx *)ctx;
261 Vec Xseq = user->Xseq;
262 VecScatter scat = user->scat;
263
264 PetscFunctionBegin;
265 /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */
266 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
267 PetscCallMPI(MPI_Comm_rank(comm, &rank));
268
269 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
270 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
271
272 fin = 0.0;
273 if (rank == 0) {
274 PetscCall(VecGetArrayRead(Xseq, &x));
275 fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
276 PetscCall(VecRestoreArrayRead(Xseq, &x));
277 }
278 PetscCallMPI(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));
279
280 /* G = 2*X - 6 */
281 PetscCall(VecSet(G, -6.0));
282 PetscCall(VecAXPY(G, 2.0, X));
283 PetscFunctionReturn(PETSC_SUCCESS);
284 }
285
286 /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708
287 H = Wxx = fxx + Jacobian(grad g^T*DE) - Jacobian(grad h^T*DI)]
288 = fxx + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T)
289 = [2 0; 0 2] + [2*de[0] 0; 0 0] - [2*di[0]-2*di[1], 0; 0, 0]
290 = [ 2*(1+de[0]-di[0]+di[1]), 0;
291 0, 2]
292 */
FormPDIPMHessian(Tao tao,Vec x,Mat H,Mat Hpre,PetscCtx ctx)293 PetscErrorCode FormPDIPMHessian(Tao tao, Vec x, Mat H, Mat Hpre, PetscCtx ctx)
294 {
295 AppCtx *user = (AppCtx *)ctx;
296 Vec DE, DI;
297 const PetscScalar *de, *di;
298 PetscInt zero = 0, one = 1;
299 PetscScalar two = 2.0;
300 PetscScalar val = 0.0;
301 Vec Deseq, Diseq;
302 VecScatter Descat, Discat;
303 PetscMPIInt rank;
304 MPI_Comm comm;
305
306 PetscFunctionBegin;
307 PetscCall(TaoGetDualVariables(tao, &DE, &DI));
308
309 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
310 PetscCallMPI(MPI_Comm_rank(comm, &rank));
311
312 if (!user->noeqflag) {
313 PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
314 PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
315 PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
316 }
317 PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
318 PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
319 PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
320
321 if (rank == 0) {
322 if (!user->noeqflag) PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */
323 PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */
324
325 if (!user->noeqflag) {
326 val = 2.0 * (1 + de[0] - di[0] + di[1]);
327 PetscCall(VecRestoreArrayRead(Deseq, &de));
328 PetscCall(VecRestoreArrayRead(Diseq, &di));
329 } else {
330 val = 2.0 * (1 - di[0] + di[1]);
331 }
332 PetscCall(VecRestoreArrayRead(Diseq, &di));
333 PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES));
334 PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES));
335 }
336 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
337 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
338 if (!user->noeqflag) {
339 PetscCall(VecScatterDestroy(&Descat));
340 PetscCall(VecDestroy(&Deseq));
341 }
342 PetscCall(VecScatterDestroy(&Discat));
343 PetscCall(VecDestroy(&Diseq));
344 PetscFunctionReturn(PETSC_SUCCESS);
345 }
346
347 /* Evaluate
348 h = [ x0^2 - x1;
349 1 -(x0^2 - x1)]
350 */
FormInequalityConstraints(Tao tao,Vec X,Vec CI,PetscCtx ctx)351 PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, PetscCtx ctx)
352 {
353 const PetscScalar *x;
354 PetscScalar ci;
355 MPI_Comm comm;
356 PetscMPIInt rank;
357 AppCtx *user = (AppCtx *)ctx;
358 Vec Xseq = user->Xseq;
359 VecScatter scat = user->scat;
360
361 PetscFunctionBegin;
362 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
363 PetscCallMPI(MPI_Comm_rank(comm, &rank));
364
365 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
366 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
367
368 if (rank == 0) {
369 PetscCall(VecGetArrayRead(Xseq, &x));
370 ci = x[0] * x[0] - x[1];
371 PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
372 ci = -x[0] * x[0] + x[1] + 1.0;
373 PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
374 PetscCall(VecRestoreArrayRead(Xseq, &x));
375 }
376 PetscCall(VecAssemblyBegin(CI));
377 PetscCall(VecAssemblyEnd(CI));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
381 /* Evaluate
382 g = [ x0^2 + x1 - 2]
383 */
FormEqualityConstraints(Tao tao,Vec X,Vec CE,PetscCtx ctx)384 PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, PetscCtx ctx)
385 {
386 const PetscScalar *x;
387 PetscScalar ce;
388 MPI_Comm comm;
389 PetscMPIInt rank;
390 AppCtx *user = (AppCtx *)ctx;
391 Vec Xseq = user->Xseq;
392 VecScatter scat = user->scat;
393
394 PetscFunctionBegin;
395 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
396 PetscCallMPI(MPI_Comm_rank(comm, &rank));
397
398 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
399 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
400
401 if (rank == 0) {
402 PetscCall(VecGetArrayRead(Xseq, &x));
403 ce = x[0] * x[0] + x[1] - 2.0;
404 PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
405 PetscCall(VecRestoreArrayRead(Xseq, &x));
406 }
407 PetscCall(VecAssemblyBegin(CE));
408 PetscCall(VecAssemblyEnd(CE));
409 PetscFunctionReturn(PETSC_SUCCESS);
410 }
411
412 /*
413 grad h = [ 2*x0, -1;
414 -2*x0, 1]
415 */
FormInequalityJacobian(Tao tao,Vec X,Mat JI,Mat JIpre,PetscCtx ctx)416 PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, PetscCtx ctx)
417 {
418 AppCtx *user = (AppCtx *)ctx;
419 PetscInt zero = 0, one = 1, cols[2];
420 PetscScalar vals[2];
421 const PetscScalar *x;
422 Vec Xseq = user->Xseq;
423 VecScatter scat = user->scat;
424 MPI_Comm comm;
425 PetscMPIInt rank;
426
427 PetscFunctionBegin;
428 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
429 PetscCallMPI(MPI_Comm_rank(comm, &rank));
430 PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
431 PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
432
433 PetscCall(VecGetArrayRead(Xseq, &x));
434 if (rank == 0) {
435 cols[0] = 0;
436 cols[1] = 1;
437 vals[0] = 2 * x[0];
438 vals[1] = -1.0;
439 PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES));
440 vals[0] = -2 * x[0];
441 vals[1] = 1.0;
442 PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES));
443 }
444 PetscCall(VecRestoreArrayRead(Xseq, &x));
445 PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY));
446 PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY));
447 PetscFunctionReturn(PETSC_SUCCESS);
448 }
449
450 /*
451 grad g = [2*x0, 1.0]
452 */
FormEqualityJacobian(Tao tao,Vec X,Mat JE,Mat JEpre,PetscCtx ctx)453 PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, PetscCtx ctx)
454 {
455 PetscInt zero = 0, cols[2];
456 PetscScalar vals[2];
457 const PetscScalar *x;
458 PetscMPIInt rank;
459 MPI_Comm comm;
460
461 PetscFunctionBegin;
462 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
463 PetscCallMPI(MPI_Comm_rank(comm, &rank));
464
465 if (rank == 0) {
466 PetscCall(VecGetArrayRead(X, &x));
467 cols[0] = 0;
468 cols[1] = 1;
469 vals[0] = 2 * x[0];
470 vals[1] = 1.0;
471 PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
472 PetscCall(VecRestoreArrayRead(X, &x));
473 }
474 PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
475 PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
476 PetscFunctionReturn(PETSC_SUCCESS);
477 }
478
479 /*TEST
480
481 build:
482 requires: !complex !defined(PETSC_USE_CXX)
483
484 test:
485 args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
486 requires: mumps !single
487 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
488
489 test:
490 suffix: pdipm_2
491 requires: mumps
492 nsize: 2
493 args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm
494 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
495
496 test:
497 suffix: 2
498 args: -tao_converged_reason
499 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
500
501 test:
502 suffix: 3
503 args: -tao_converged_reason -no_eq
504 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
505
506 test:
507 suffix: 4
508 args: -tao_converged_reason -tao_almm_type classic
509 requires: !single
510 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
511
512 test:
513 suffix: 5
514 args: -tao_converged_reason -tao_almm_type classic -no_eq
515 requires: !single
516 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
517
518 test:
519 suffix: 6
520 args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
521 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
522
523 test:
524 suffix: 7
525 args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
526 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
527
528 test:
529 suffix: 8
530 nsize: 2
531 args: -tao_converged_reason
532 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
533
534 test:
535 suffix: 9
536 nsize: 2
537 args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
538 requires: cuda
539 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
540
541 test:
542 suffix: 10
543 args: -tao_converged_reason -tao_almm_type classic -no_eq -no_bound
544 requires: !single
545 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
546
547 test:
548 suffix: 11
549 args: -tao_converged_reason -tao_almm_type classic -no_bound
550 requires: !single
551 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
552
553 test:
554 suffix: 12
555 args: -tao_converged_reason -tao_almm_type phr -no_eq -no_bound
556 requires: !single
557 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
558
559 test:
560 suffix: 13
561 args: -tao_converged_reason -tao_almm_type phr -no_bound
562 requires: !single
563 filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
564
565 TEST*/
566