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: 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 */
62: /* Initialize TAO,PETSc */
63: PetscFunctionBeginUser;
64: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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));
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));
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));
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: }
107: PetscCall(DestroyProblem(&user));
108: PetscCall(VecDestroy(&x));
109: PetscCall(VecDestroy(&ceq));
110: PetscCall(VecDestroy(&cin));
111: PetscCall(TaoDestroy(&tao));
113: PetscCall(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: 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));
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;
141: PetscCall(PetscStrncpy(filename, filebase, sizeof(filename)));
142: PetscCall(PetscStrlcat(filename, "H", sizeof(filename)));
143: PetscCall(PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader));
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));
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;
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));
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));
181: PetscCall(MatMPIAIJSetPreallocation(user->Ain, 1, NULL, 0, NULL));
182: PetscCall(MatSeqAIJSetPreallocation(user->Ain, 1, NULL));
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));
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: }
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: }
211: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
212: {
213: AppCtx *user = (AppCtx *)ctx;
214: PetscScalar xtHx;
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: }
225: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
226: {
227: PetscFunctionBegin;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
232: {
233: AppCtx *user = (AppCtx *)ctx;
235: PetscFunctionBegin;
236: PetscCall(MatMult(user->Ain, x, ci));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce, void *ctx)
241: {
242: AppCtx *user = (AppCtx *)ctx;
244: PetscFunctionBegin;
245: PetscCall(MatMult(user->Aeq, x, ce));
246: PetscCall(VecAXPY(ce, -1.0, user->beq));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
251: {
252: PetscFunctionBegin;
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
257: {
258: PetscFunctionBegin;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*TEST
264: build:
265: requires: !complex
267: test:
268: requires: !single superlu
269: localrunfiles: HS21
271: TEST*/