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 optimiztion 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: -init_view : view initial object setup\n\
25: -snes_fd : snes with finite difference Jacobian (needed for pdipm)\n\
26: -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
27: -tao_cmonitor : convergence monitor with constraint norm \n\
28: -tao_view_solution : view exact solution at each iteration\n\
29: 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: User-defined application context - contains data needed by the application
33: */
34: typedef struct {
35: PetscInt n; /* Global length of x */
36: PetscInt ne; /* Global number of equality constraints */
37: PetscInt ni; /* Global number of inequality constraints */
38: PetscBool noeqflag, initview;
39: Vec x, xl, xu;
40: Vec ce, ci, bl, bu, Xseq;
41: Mat Ae, Ai, H;
42: VecScatter scat;
43: } AppCtx;
45: /* -------- User-defined Routines --------- */
46: PetscErrorCode InitializeProblem(AppCtx *);
47: PetscErrorCode DestroyProblem(AppCtx *);
48: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
49: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
50: PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
51: PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
52: PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
53: PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);
55: PetscErrorCode main(int argc, char **argv)
56: {
57: Tao tao;
58: KSP ksp;
59: PC pc;
60: AppCtx user; /* application context */
61: Vec x, G, CI, CE;
62: PetscMPIInt size;
63: TaoType type;
64: PetscReal f;
65: PetscBool pdipm, mumps;
68: PetscInitialize(&argc, &argv, (char *)0, help);
69: MPI_Comm_size(PETSC_COMM_WORLD, &size);
72: PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n");
73: InitializeProblem(&user); /* sets up problem, function below */
75: TaoCreate(PETSC_COMM_WORLD, &tao);
76: TaoSetType(tao, TAOALMM);
77: TaoSetSolution(tao, user.x);
78: TaoSetVariableBounds(tao, user.xl, user.xu);
79: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
80: TaoSetTolerances(tao, 1.e-4, 0.0, 0.0);
81: TaoSetConstraintTolerances(tao, 1.e-4, 0.0);
82: TaoSetMaximumFunctionEvaluations(tao, 1e6);
83: TaoSetFromOptions(tao);
85: if (!user.noeqflag) {
86: TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user);
87: TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user);
88: }
89: TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user);
90: TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user);
92: TaoGetType(tao, &type);
93: PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm);
94: if (pdipm) {
95: /*
96: PDIPM produces an indefinite KKT matrix with some zeros along the diagonal
97: Inverting this indefinite matrix requires PETSc to be configured with MUMPS
98: */
99: PetscHasExternalPackage("mumps", &mumps);
101: TaoGetKSP(tao, &ksp);
102: KSPGetPC(ksp, &pc);
103: PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
104: KSPSetType(ksp, KSPPREONLY);
105: KSPSetFromOptions(ksp);
106: TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
107: }
109: /* Print out an initial view of the problem */
110: if (user.initview) {
111: TaoSetUp(tao);
112: VecDuplicate(user.x, &G);
113: FormFunctionGradient(tao, user.x, &f, G, (void *)&user);
114: FormHessian(tao, user.x, user.H, user.H, (void *)&user);
115: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD);
116: PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n");
117: VecView(user.x, PETSC_VIEWER_STDOUT_WORLD);
118: PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f);
119: PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n");
120: VecView(G, PETSC_VIEWER_STDOUT_WORLD);
121: MatView(user.H, PETSC_VIEWER_STDOUT_WORLD);
122: VecDestroy(&G);
123: FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user);
124: MatCreateVecs(user.Ai, NULL, &CI);
125: FormInequalityConstraints(tao, user.x, CI, (void *)&user);
126: PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n");
127: VecView(CI, PETSC_VIEWER_STDOUT_WORLD);
128: MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD);
129: VecDestroy(&CI);
130: if (!user.noeqflag) {
131: FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user);
132: MatCreateVecs(user.Ae, NULL, &CE);
133: FormEqualityConstraints(tao, user.x, CE, (void *)&user);
134: PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n");
135: VecView(CE, PETSC_VIEWER_STDOUT_WORLD);
136: MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD);
137: VecDestroy(&CE);
138: }
139: PetscPrintf(PETSC_COMM_WORLD, "\n");
140: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD);
141: }
143: TaoSolve(tao);
144: TaoGetSolution(tao, &x);
145: PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n");
146: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
148: /* Free objects */
149: DestroyProblem(&user);
150: TaoDestroy(&tao);
151: PetscFinalize();
152: return 0;
153: }
155: PetscErrorCode InitializeProblem(AppCtx *user)
156: {
157: PetscMPIInt size;
158: PetscMPIInt rank;
159: PetscInt nloc, neloc, niloc;
161: MPI_Comm_size(PETSC_COMM_WORLD, &size);
162: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
163: user->noeqflag = PETSC_FALSE;
164: PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL);
165: user->initview = PETSC_FALSE;
166: PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL);
168: /* Tell user the correct solution, not an error checking */
169: if (!user->noeqflag) {
170: PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n");
171: } else {
172: PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n");
173: }
175: /* create vector x and set initial values */
176: user->n = 2; /* global length */
177: nloc = (size == 1) ? user->n : 1;
178: VecCreate(PETSC_COMM_WORLD, &user->x);
179: VecSetSizes(user->x, nloc, user->n);
180: VecSetFromOptions(user->x);
181: VecSet(user->x, 0.0);
183: /* create and set lower and upper bound vectors */
184: VecDuplicate(user->x, &user->xl);
185: VecDuplicate(user->x, &user->xu);
186: VecSet(user->xl, -1.0);
187: VecSet(user->xu, 2.0);
189: /* create scater to zero */
190: VecScatterCreateToZero(user->x, &user->scat, &user->Xseq);
192: user->ne = 1;
193: user->ni = 2;
194: neloc = (rank == 0) ? user->ne : 0;
195: niloc = (size == 1) ? user->ni : 1;
197: if (!user->noeqflag) {
198: VecCreate(PETSC_COMM_WORLD, &user->ce); /* a 1x1 vec for equality constraints */
199: VecSetSizes(user->ce, neloc, user->ne);
200: VecSetFromOptions(user->ce);
201: VecSetUp(user->ce);
202: }
204: VecCreate(PETSC_COMM_WORLD, &user->ci); /* a 2x1 vec for inequality constraints */
205: VecSetSizes(user->ci, niloc, user->ni);
206: VecSetFromOptions(user->ci);
207: VecSetUp(user->ci);
209: /* nexn & nixn matricies for equally and inequalty constraints */
210: if (!user->noeqflag) {
211: MatCreate(PETSC_COMM_WORLD, &user->Ae);
212: MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n);
213: MatSetFromOptions(user->Ae);
214: MatSetUp(user->Ae);
215: }
217: MatCreate(PETSC_COMM_WORLD, &user->Ai);
218: MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n);
219: MatSetFromOptions(user->Ai);
220: MatSetUp(user->Ai);
222: MatCreate(PETSC_COMM_WORLD, &user->H);
223: MatSetSizes(user->H, nloc, nloc, user->n, user->n);
224: MatSetFromOptions(user->H);
225: MatSetUp(user->H);
226: return 0;
227: }
229: PetscErrorCode DestroyProblem(AppCtx *user)
230: {
231: if (!user->noeqflag) MatDestroy(&user->Ae);
232: MatDestroy(&user->Ai);
233: MatDestroy(&user->H);
235: VecDestroy(&user->x);
236: if (!user->noeqflag) VecDestroy(&user->ce);
237: VecDestroy(&user->ci);
238: VecDestroy(&user->xl);
239: VecDestroy(&user->xu);
240: VecDestroy(&user->Xseq);
241: VecScatterDestroy(&user->scat);
242: return 0;
243: }
245: /* Evaluate
246: f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
247: G = grad f = [2*(x0 - 2) - 2;
248: 2*(x1 - 2) - 2]
249: */
250: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
251: {
252: PetscScalar g;
253: const PetscScalar *x;
254: MPI_Comm comm;
255: PetscMPIInt rank;
256: PetscReal fin;
257: AppCtx *user = (AppCtx *)ctx;
258: Vec Xseq = user->Xseq;
259: VecScatter scat = user->scat;
261: PetscObjectGetComm((PetscObject)tao, &comm);
262: MPI_Comm_rank(comm, &rank);
264: VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
265: VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
267: fin = 0.0;
268: if (rank == 0) {
269: VecGetArrayRead(Xseq, &x);
270: fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
271: g = 2.0 * (x[0] - 2.0) - 2.0;
272: VecSetValue(G, 0, g, INSERT_VALUES);
273: g = 2.0 * (x[1] - 2.0) - 2.0;
274: VecSetValue(G, 1, g, INSERT_VALUES);
275: VecRestoreArrayRead(Xseq, &x);
276: }
277: MPI_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm);
278: VecAssemblyBegin(G);
279: VecAssemblyEnd(G);
280: return 0;
281: }
283: /* Evaluate
284: H = fxx + grad (grad g^T*DI) - grad (grad h^T*DE)]
285: = [ 2*(1+de[0]-di[0]+di[1]), 0;
286: 0, 2]
287: */
288: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
289: {
290: AppCtx *user = (AppCtx *)ctx;
291: Vec DE, DI;
292: const PetscScalar *de, *di;
293: PetscInt zero = 0, one = 1;
294: PetscScalar two = 2.0;
295: PetscScalar val = 0.0;
296: Vec Deseq, Diseq;
297: VecScatter Descat, Discat;
298: PetscMPIInt rank;
299: MPI_Comm comm;
301: TaoGetDualVariables(tao, &DE, &DI);
303: PetscObjectGetComm((PetscObject)tao, &comm);
304: MPI_Comm_rank(comm, &rank);
306: if (!user->noeqflag) {
307: VecScatterCreateToZero(DE, &Descat, &Deseq);
308: VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD);
309: VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD);
310: }
311: VecScatterCreateToZero(DI, &Discat, &Diseq);
312: VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD);
313: VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD);
315: if (rank == 0) {
316: if (!user->noeqflag) { VecGetArrayRead(Deseq, &de); /* places equality constraint dual into array */ }
317: VecGetArrayRead(Diseq, &di); /* places inequality constraint dual into array */
319: if (!user->noeqflag) {
320: val = 2.0 * (1 + de[0] - di[0] + di[1]);
321: VecRestoreArrayRead(Deseq, &de);
322: VecRestoreArrayRead(Diseq, &di);
323: } else {
324: val = 2.0 * (1 - di[0] + di[1]);
325: }
326: VecRestoreArrayRead(Diseq, &di);
327: MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES);
328: MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES);
329: }
330: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
331: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
332: if (!user->noeqflag) {
333: VecScatterDestroy(&Descat);
334: VecDestroy(&Deseq);
335: }
336: VecScatterDestroy(&Discat);
337: VecDestroy(&Diseq);
338: return 0;
339: }
341: /* Evaluate
342: h = [ x0^2 - x1;
343: 1 -(x0^2 - x1)]
344: */
345: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
346: {
347: const PetscScalar *x;
348: PetscScalar ci;
349: MPI_Comm comm;
350: PetscMPIInt rank;
351: AppCtx *user = (AppCtx *)ctx;
352: Vec Xseq = user->Xseq;
353: VecScatter scat = user->scat;
355: PetscObjectGetComm((PetscObject)tao, &comm);
356: MPI_Comm_rank(comm, &rank);
358: VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
359: VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
361: if (rank == 0) {
362: VecGetArrayRead(Xseq, &x);
363: ci = x[0] * x[0] - x[1];
364: VecSetValue(CI, 0, ci, INSERT_VALUES);
365: ci = -x[0] * x[0] + x[1] + 1.0;
366: VecSetValue(CI, 1, ci, INSERT_VALUES);
367: VecRestoreArrayRead(Xseq, &x);
368: }
369: VecAssemblyBegin(CI);
370: VecAssemblyEnd(CI);
371: return 0;
372: }
374: /* Evaluate
375: g = [ x0^2 + x1 - 2]
376: */
377: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx)
378: {
379: const PetscScalar *x;
380: PetscScalar ce;
381: MPI_Comm comm;
382: PetscMPIInt rank;
383: AppCtx *user = (AppCtx *)ctx;
384: Vec Xseq = user->Xseq;
385: VecScatter scat = user->scat;
387: PetscObjectGetComm((PetscObject)tao, &comm);
388: MPI_Comm_rank(comm, &rank);
390: VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
391: VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
393: if (rank == 0) {
394: VecGetArrayRead(Xseq, &x);
395: ce = x[0] * x[0] + x[1] - 2.0;
396: VecSetValue(CE, 0, ce, INSERT_VALUES);
397: VecRestoreArrayRead(Xseq, &x);
398: }
399: VecAssemblyBegin(CE);
400: VecAssemblyEnd(CE);
401: return 0;
402: }
404: /*
405: grad h = [ 2*x0, -1;
406: -2*x0, 1]
407: */
408: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
409: {
410: AppCtx *user = (AppCtx *)ctx;
411: PetscInt zero = 0, one = 1, cols[2];
412: PetscScalar vals[2];
413: const PetscScalar *x;
414: Vec Xseq = user->Xseq;
415: VecScatter scat = user->scat;
416: MPI_Comm comm;
417: PetscMPIInt rank;
419: PetscObjectGetComm((PetscObject)tao, &comm);
420: MPI_Comm_rank(comm, &rank);
421: VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
422: VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD);
424: VecGetArrayRead(Xseq, &x);
425: if (rank == 0) {
426: cols[0] = 0;
427: cols[1] = 1;
428: vals[0] = 2 * x[0];
429: vals[1] = -1.0;
430: MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES);
431: vals[0] = -2 * x[0];
432: vals[1] = 1.0;
433: MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES);
434: }
435: VecRestoreArrayRead(Xseq, &x);
436: MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY);
437: MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY);
438: return 0;
439: }
441: /*
442: grad g = [2*x0
443: 1.0 ]
444: */
445: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
446: {
447: PetscInt zero = 0, cols[2];
448: PetscScalar vals[2];
449: const PetscScalar *x;
450: PetscMPIInt rank;
451: MPI_Comm comm;
453: PetscObjectGetComm((PetscObject)tao, &comm);
454: MPI_Comm_rank(comm, &rank);
456: if (rank == 0) {
457: VecGetArrayRead(X, &x);
458: cols[0] = 0;
459: cols[1] = 1;
460: vals[0] = 2 * x[0];
461: vals[1] = 1.0;
462: MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES);
463: VecRestoreArrayRead(X, &x);
464: }
465: MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY);
466: MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY);
467: return 0;
468: }
470: /*TEST
472: build:
473: requires: !complex !defined(PETSC_USE_CXX)
475: test:
476: args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
477: requires: mumps
478: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
480: test:
481: suffix: 2
482: args: -tao_converged_reason
483: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
485: test:
486: suffix: 3
487: args: -tao_converged_reason -no_eq
488: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
490: test:
491: suffix: 4
492: args: -tao_converged_reason -tao_almm_type classic
493: requires: !single
494: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
496: test:
497: suffix: 5
498: args: -tao_converged_reason -tao_almm_type classic -no_eq
499: requires: !single
500: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
502: test:
503: suffix: 6
504: args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
505: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
507: test:
508: suffix: 7
509: args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
510: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
512: test:
513: suffix: 8
514: nsize: 2
515: args: -tao_converged_reason
516: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
518: test:
519: suffix: 9
520: nsize: 2
521: args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
522: requires: cuda
523: filter: sed -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"
525: TEST*/