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