Actual source code: maros.c
1: /* Program usage: mpiexec -n 1 maros1 [-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
4: TODO Explain maros example
5: ---------------------------------------------------------------------- */
7: #include <petsctao.h>
9: static char help[] = "";
11: /*
12: User-defined application context - contains data needed by the
13: application-provided call-back routines, FormFunction(),
14: FormGradient(), and FormHessian().
15: */
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;
39: /* -------- User-defined Routines --------- */
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 *);
50: PetscErrorCode 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 */
62: /* Initialize TAO,PETSc */
64: PetscInitialize(&argc, &argv, (char *)0, help);
65: MPI_Comm_size(PETSC_COMM_WORLD, &size);
66: /* Specify default parameters for the problem, check for command-line overrides */
67: PetscStrncpy(user.name, "HS21", sizeof(user.name));
68: PetscOptionsGetString(NULL, NULL, "-cutername", user.name, sizeof(user.name), &flg);
70: PetscPrintf(PETSC_COMM_WORLD, "\n---- MAROS Problem %s -----\n", user.name);
71: InitializeProblem(&user);
72: VecDuplicate(user.d, &x);
73: VecDuplicate(user.beq, &ceq);
74: VecDuplicate(user.bin, &cin);
75: VecSet(x, 1.0);
77: TaoCreate(PETSC_COMM_WORLD, &tao);
78: TaoSetType(tao, TAOIPM);
79: TaoSetSolution(tao, x);
80: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
81: TaoSetEqualityConstraintsRoutine(tao, ceq, FormEqualityConstraints, (void *)&user);
82: TaoSetInequalityConstraintsRoutine(tao, cin, FormInequalityConstraints, (void *)&user);
83: TaoSetInequalityBounds(tao, user.bin, NULL);
84: TaoSetJacobianEqualityRoutine(tao, user.Aeq, user.Aeq, FormEqualityJacobian, (void *)&user);
85: TaoSetJacobianInequalityRoutine(tao, user.Ain, user.Ain, FormInequalityJacobian, (void *)&user);
86: TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
87: TaoGetKSP(tao, &ksp);
88: KSPGetPC(ksp, &pc);
89: 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: PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU);
95: KSPSetType(ksp, KSPPREONLY);
96: TaoSetTolerances(tao, 0, 0, 0);
98: TaoSetFromOptions(tao);
99: TaoSolve(tao);
100: TaoGetConvergedReason(tao, &reason);
101: if (reason < 0) {
102: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n", TaoConvergedReasons[reason]);
103: } else {
104: PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n", TaoConvergedReasons[reason]);
105: }
107: DestroyProblem(&user);
108: VecDestroy(&x);
109: VecDestroy(&ceq);
110: VecDestroy(&cin);
111: TaoDestroy(&tao);
113: PetscFinalize();
114: return 0;
115: }
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];
126: comm = PETSC_COMM_WORLD;
127: PetscStrncpy(filebase, user->name, sizeof(filebase));
128: PetscStrlcat(filebase, "/", sizeof(filebase));
129: PetscStrncpy(filename, filebase, sizeof(filename));
130: PetscStrlcat(filename, "f", sizeof(filename));
131: PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
133: VecCreate(comm, &user->d);
134: VecLoad(user->d, loader);
135: PetscViewerDestroy(&loader);
136: VecGetSize(user->d, &nrows);
137: VecSetFromOptions(user->d);
138: user->n = nrows;
140: PetscStrncpy(filename, filebase, sizeof(filename));
141: PetscStrlcat(filename, "H", sizeof(filename));
142: PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
144: MatCreate(comm, &user->H);
145: MatSetSizes(user->H, PETSC_DECIDE, PETSC_DECIDE, nrows, nrows);
146: MatLoad(user->H, loader);
147: PetscViewerDestroy(&loader);
148: MatGetSize(user->H, &nrows, &ncols);
151: MatSetFromOptions(user->H);
153: PetscStrncpy(filename, filebase, sizeof(filename));
154: PetscStrlcat(filename, "Aeq", sizeof(filename));
155: PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
156: MatCreate(comm, &user->Aeq);
157: MatLoad(user->Aeq, loader);
158: PetscViewerDestroy(&loader);
159: MatGetSize(user->Aeq, &nrows, &ncols);
161: MatSetFromOptions(user->Aeq);
162: user->me = nrows;
164: PetscStrncpy(filename, filebase, sizeof(filename));
165: PetscStrlcat(filename, "Beq", sizeof(filename));
166: PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
167: VecCreate(comm, &user->beq);
168: VecLoad(user->beq, loader);
169: PetscViewerDestroy(&loader);
170: VecGetSize(user->beq, &nrows);
172: VecSetFromOptions(user->beq);
174: user->mi = user->n;
175: /* Ain = eye(n,n) */
176: MatCreate(comm, &user->Ain);
177: MatSetType(user->Ain, MATAIJ);
178: MatSetSizes(user->Ain, PETSC_DECIDE, PETSC_DECIDE, user->mi, user->mi);
180: MatMPIAIJSetPreallocation(user->Ain, 1, NULL, 0, NULL);
181: MatSeqAIJSetPreallocation(user->Ain, 1, NULL);
183: for (i = 0; i < user->mi; i++) MatSetValues(user->Ain, 1, &i, 1, &i, &one, INSERT_VALUES);
184: MatAssemblyBegin(user->Ain, MAT_FINAL_ASSEMBLY);
185: MatAssemblyEnd(user->Ain, MAT_FINAL_ASSEMBLY);
186: MatSetFromOptions(user->Ain);
188: /* bin = [0,0 ... 0]' */
189: VecCreate(comm, &user->bin);
190: VecSetType(user->bin, VECMPI);
191: VecSetSizes(user->bin, PETSC_DECIDE, user->mi);
192: VecSet(user->bin, 0.0);
193: VecSetFromOptions(user->bin);
194: user->m = user->me + user->mi;
195: return 0;
196: }
198: PetscErrorCode DestroyProblem(AppCtx *user)
199: {
200: MatDestroy(&user->H);
201: MatDestroy(&user->Aeq);
202: MatDestroy(&user->Ain);
203: VecDestroy(&user->beq);
204: VecDestroy(&user->bin);
205: VecDestroy(&user->d);
206: return 0;
207: }
209: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
210: {
211: AppCtx *user = (AppCtx *)ctx;
212: PetscScalar xtHx;
214: MatMult(user->H, x, g);
215: VecDot(x, g, &xtHx);
216: VecDot(x, user->d, f);
217: *f += 0.5 * xtHx;
218: VecAXPY(g, 1.0, user->d);
219: return 0;
220: }
222: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
223: {
224: return 0;
225: }
227: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
228: {
229: AppCtx *user = (AppCtx *)ctx;
231: MatMult(user->Ain, x, ci);
232: return 0;
233: }
235: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce, void *ctx)
236: {
237: AppCtx *user = (AppCtx *)ctx;
239: MatMult(user->Aeq, x, ce);
240: VecAXPY(ce, -1.0, user->beq);
241: return 0;
242: }
244: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
245: {
246: return 0;
247: }
249: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
250: {
251: return 0;
252: }
254: /*TEST
256: build:
257: requires: !complex
259: test:
260: requires: superlu
261: localrunfiles: HS21
263: TEST*/