Actual source code: ex20opt_ic.c
1: static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n";
3: /*
4: Notes:
5: This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS.
6: The nonlinear problem is written in an ODE equivalent form.
7: The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions.
8: The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details.
9: */
11: #include <petsctao.h>
12: #include <petscts.h>
14: typedef struct _n_User *User;
15: struct _n_User {
16: TS ts;
17: PetscReal mu;
18: PetscReal next_output;
20: /* Sensitivity analysis support */
21: PetscInt steps;
22: PetscReal ftime;
23: Mat A; /* Jacobian matrix for ODE */
24: Mat Jacp; /* JacobianP matrix for ODE*/
25: Mat H; /* Hessian matrix for optimization */
26: Vec U, Lambda[1], Mup[1]; /* first-order adjoint variables */
27: Vec Lambda2[2]; /* second-order adjoint variables */
28: Vec Ihp1[1]; /* working space for Hessian evaluations */
29: Vec Dir; /* direction vector */
30: PetscReal ob[2]; /* observation used by the cost function */
31: PetscBool implicitform; /* implicit ODE? */
32: };
33: PetscErrorCode Adjoint2(Vec, PetscScalar[], User);
35: /* ----------------------- Explicit form of the ODE -------------------- */
37: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
38: {
39: User user = (User)ctx;
40: PetscScalar *f;
41: const PetscScalar *u;
43: PetscFunctionBeginUser;
44: PetscCall(VecGetArrayRead(U, &u));
45: PetscCall(VecGetArray(F, &f));
46: f[0] = u[1];
47: f[1] = user->mu * ((1. - u[0] * u[0]) * u[1] - u[0]);
48: PetscCall(VecRestoreArrayRead(U, &u));
49: PetscCall(VecRestoreArray(F, &f));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
54: {
55: User user = (User)ctx;
56: PetscReal mu = user->mu;
57: PetscInt rowcol[] = {0, 1};
58: PetscScalar J[2][2];
59: const PetscScalar *u;
61: PetscFunctionBeginUser;
62: PetscCall(VecGetArrayRead(U, &u));
63: J[0][0] = 0;
64: J[1][0] = -mu * (2.0 * u[1] * u[0] + 1.);
65: J[0][1] = 1.0;
66: J[1][1] = mu * (1.0 - u[0] * u[0]);
67: PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
68: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
69: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
70: if (A != B) {
71: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
73: }
74: PetscCall(VecRestoreArrayRead(U, &u));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
79: {
80: const PetscScalar *vl, *vr, *u;
81: PetscScalar *vhv;
82: PetscScalar dJdU[2][2][2] = {{{0}}};
83: PetscInt i, j, k;
84: User user = (User)ctx;
86: PetscFunctionBeginUser;
87: PetscCall(VecGetArrayRead(U, &u));
88: PetscCall(VecGetArrayRead(Vl[0], &vl));
89: PetscCall(VecGetArrayRead(Vr, &vr));
90: PetscCall(VecGetArray(VHV[0], &vhv));
92: dJdU[1][0][0] = -2. * user->mu * u[1];
93: dJdU[1][1][0] = -2. * user->mu * u[0];
94: dJdU[1][0][1] = -2. * user->mu * u[0];
95: for (j = 0; j < 2; j++) {
96: vhv[j] = 0;
97: for (k = 0; k < 2; k++)
98: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
99: }
100: PetscCall(VecRestoreArrayRead(U, &u));
101: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
102: PetscCall(VecRestoreArrayRead(Vr, &vr));
103: PetscCall(VecRestoreArray(VHV[0], &vhv));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /* ----------------------- Implicit form of the ODE -------------------- */
109: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
110: {
111: User user = (User)ctx;
112: const PetscScalar *u, *udot;
113: PetscScalar *f;
115: PetscFunctionBeginUser;
116: PetscCall(VecGetArrayRead(U, &u));
117: PetscCall(VecGetArrayRead(Udot, &udot));
118: PetscCall(VecGetArray(F, &f));
119: f[0] = udot[0] - u[1];
120: f[1] = udot[1] - user->mu * ((1.0 - u[0] * u[0]) * u[1] - u[0]);
121: PetscCall(VecRestoreArrayRead(U, &u));
122: PetscCall(VecRestoreArrayRead(Udot, &udot));
123: PetscCall(VecRestoreArray(F, &f));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
128: {
129: User user = (User)ctx;
130: PetscInt rowcol[] = {0, 1};
131: PetscScalar J[2][2];
132: const PetscScalar *u;
134: PetscFunctionBeginUser;
135: PetscCall(VecGetArrayRead(U, &u));
136: J[0][0] = a;
137: J[0][1] = -1.0;
138: J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
139: J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
140: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
141: PetscCall(VecRestoreArrayRead(U, &u));
142: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
143: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
144: if (A != B) {
145: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
146: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
152: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
153: {
154: const PetscScalar *u;
155: PetscReal tfinal, dt;
156: User user = (User)ctx;
157: Vec interpolatedU;
159: PetscFunctionBeginUser;
160: PetscCall(TSGetTimeStep(ts, &dt));
161: PetscCall(TSGetMaxTime(ts, &tfinal));
163: while (user->next_output <= t && user->next_output <= tfinal) {
164: PetscCall(VecDuplicate(U, &interpolatedU));
165: PetscCall(TSInterpolate(ts, user->next_output, interpolatedU));
166: PetscCall(VecGetArrayRead(interpolatedU, &u));
167: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%g] %" PetscInt_FMT " TS %g (dt = %g) X %g %g\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(u[0]), (double)PetscRealPart(u[1])));
168: PetscCall(VecRestoreArrayRead(interpolatedU, &u));
169: PetscCall(VecDestroy(&interpolatedU));
170: user->next_output += 0.1;
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode IHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
176: {
177: const PetscScalar *vl, *vr, *u;
178: PetscScalar *vhv;
179: PetscScalar dJdU[2][2][2] = {{{0}}};
180: PetscInt i, j, k;
181: User user = (User)ctx;
183: PetscFunctionBeginUser;
184: PetscCall(VecGetArrayRead(U, &u));
185: PetscCall(VecGetArrayRead(Vl[0], &vl));
186: PetscCall(VecGetArrayRead(Vr, &vr));
187: PetscCall(VecGetArray(VHV[0], &vhv));
188: dJdU[1][0][0] = 2. * user->mu * u[1];
189: dJdU[1][1][0] = 2. * user->mu * u[0];
190: dJdU[1][0][1] = 2. * user->mu * u[0];
191: for (j = 0; j < 2; j++) {
192: vhv[j] = 0;
193: for (k = 0; k < 2; k++)
194: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
195: }
196: PetscCall(VecRestoreArrayRead(U, &u));
197: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
198: PetscCall(VecRestoreArrayRead(Vr, &vr));
199: PetscCall(VecRestoreArray(VHV[0], &vhv));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /* ------------------ User-defined routines for TAO -------------------------- */
205: static PetscErrorCode FormFunctionGradient(Tao tao, Vec IC, PetscReal *f, Vec G, void *ctx)
206: {
207: User user_ptr = (User)ctx;
208: TS ts = user_ptr->ts;
209: const PetscScalar *x_ptr;
210: PetscScalar *y_ptr;
212: PetscFunctionBeginUser;
213: PetscCall(VecCopy(IC, user_ptr->U)); /* set up the initial condition */
215: PetscCall(TSSetTime(ts, 0.0));
216: PetscCall(TSSetStepNumber(ts, 0));
217: PetscCall(TSSetTimeStep(ts, 0.001)); /* can be overwritten by command line options */
218: PetscCall(TSSetFromOptions(ts));
219: PetscCall(TSSolve(ts, user_ptr->U));
221: PetscCall(VecGetArrayRead(user_ptr->U, &x_ptr));
222: PetscCall(VecGetArray(user_ptr->Lambda[0], &y_ptr));
223: *f = (x_ptr[0] - user_ptr->ob[0]) * (x_ptr[0] - user_ptr->ob[0]) + (x_ptr[1] - user_ptr->ob[1]) * (x_ptr[1] - user_ptr->ob[1]);
224: y_ptr[0] = 2. * (x_ptr[0] - user_ptr->ob[0]);
225: y_ptr[1] = 2. * (x_ptr[1] - user_ptr->ob[1]);
226: PetscCall(VecRestoreArray(user_ptr->Lambda[0], &y_ptr));
227: PetscCall(VecRestoreArrayRead(user_ptr->U, &x_ptr));
229: PetscCall(TSSetCostGradients(ts, 1, user_ptr->Lambda, NULL));
230: PetscCall(TSAdjointSolve(ts));
231: PetscCall(VecCopy(user_ptr->Lambda[0], G));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: static PetscErrorCode FormHessian(Tao tao, Vec U, Mat H, Mat Hpre, void *ctx)
236: {
237: User user_ptr = (User)ctx;
238: PetscScalar harr[2];
239: PetscScalar *x_ptr;
240: const PetscInt rows[2] = {0, 1};
241: PetscInt col;
243: PetscFunctionBeginUser;
244: PetscCall(VecCopy(U, user_ptr->U));
245: PetscCall(VecGetArray(user_ptr->Dir, &x_ptr));
246: x_ptr[0] = 1.;
247: x_ptr[1] = 0.;
248: PetscCall(VecRestoreArray(user_ptr->Dir, &x_ptr));
249: PetscCall(Adjoint2(user_ptr->U, harr, user_ptr));
250: col = 0;
251: PetscCall(MatSetValues(H, 2, rows, 1, &col, harr, INSERT_VALUES));
253: PetscCall(VecCopy(U, user_ptr->U));
254: PetscCall(VecGetArray(user_ptr->Dir, &x_ptr));
255: x_ptr[0] = 0.;
256: x_ptr[1] = 1.;
257: PetscCall(VecRestoreArray(user_ptr->Dir, &x_ptr));
258: PetscCall(Adjoint2(user_ptr->U, harr, user_ptr));
259: col = 1;
260: PetscCall(MatSetValues(H, 2, rows, 1, &col, harr, INSERT_VALUES));
262: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
263: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
264: if (H != Hpre) {
265: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
266: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
267: }
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: static PetscErrorCode MatrixFreeHessian(Tao tao, Vec U, Mat H, Mat Hpre, void *ctx)
272: {
273: User user_ptr = (User)ctx;
275: PetscFunctionBeginUser;
276: PetscCall(VecCopy(U, user_ptr->U));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /* ------------ Routines calculating second-order derivatives -------------- */
282: /*
283: Compute the Hessian-vector product for the cost function using Second-order adjoint
284: */
285: PetscErrorCode Adjoint2(Vec U, PetscScalar arr[], User ctx)
286: {
287: TS ts = ctx->ts;
288: PetscScalar *x_ptr, *y_ptr;
289: Mat tlmsen;
291: PetscFunctionBeginUser;
292: PetscCall(TSAdjointReset(ts));
294: PetscCall(TSSetTime(ts, 0.0));
295: PetscCall(TSSetStepNumber(ts, 0));
296: PetscCall(TSSetTimeStep(ts, 0.001));
297: PetscCall(TSSetFromOptions(ts));
298: PetscCall(TSSetCostHessianProducts(ts, 1, ctx->Lambda2, NULL, ctx->Dir));
299: PetscCall(TSAdjointSetForward(ts, NULL));
300: PetscCall(TSSolve(ts, U));
302: /* Set terminal conditions for first- and second-order adjonts */
303: PetscCall(VecGetArray(U, &x_ptr));
304: PetscCall(VecGetArray(ctx->Lambda[0], &y_ptr));
305: y_ptr[0] = 2. * (x_ptr[0] - ctx->ob[0]);
306: y_ptr[1] = 2. * (x_ptr[1] - ctx->ob[1]);
307: PetscCall(VecRestoreArray(ctx->Lambda[0], &y_ptr));
308: PetscCall(VecRestoreArray(U, &x_ptr));
309: PetscCall(TSForwardGetSensitivities(ts, NULL, &tlmsen));
310: PetscCall(MatDenseGetColumn(tlmsen, 0, &x_ptr));
311: PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
312: y_ptr[0] = 2. * x_ptr[0];
313: y_ptr[1] = 2. * x_ptr[1];
314: PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
315: PetscCall(MatDenseRestoreColumn(tlmsen, &x_ptr));
317: PetscCall(TSSetCostGradients(ts, 1, ctx->Lambda, NULL));
318: if (ctx->implicitform) {
319: PetscCall(TSSetIHessianProduct(ts, ctx->Ihp1, IHessianProductUU, NULL, NULL, NULL, NULL, NULL, NULL, ctx));
320: } else {
321: PetscCall(TSSetRHSHessianProduct(ts, ctx->Ihp1, RHSHessianProductUU, NULL, NULL, NULL, NULL, NULL, NULL, ctx));
322: }
323: PetscCall(TSAdjointSolve(ts));
325: PetscCall(VecGetArray(ctx->Lambda2[0], &x_ptr));
326: arr[0] = x_ptr[0];
327: arr[1] = x_ptr[1];
328: PetscCall(VecRestoreArray(ctx->Lambda2[0], &x_ptr));
330: PetscCall(TSAdjointReset(ts));
331: PetscCall(TSAdjointResetForward(ts));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: PetscErrorCode FiniteDiff(Vec U, PetscScalar arr[], User ctx)
336: {
337: Vec Up, G, Gp;
338: const PetscScalar eps = PetscRealConstant(1e-7);
339: PetscScalar *u;
340: Tao tao = NULL;
341: PetscReal f;
343: PetscFunctionBeginUser;
344: PetscCall(VecDuplicate(U, &Up));
345: PetscCall(VecDuplicate(U, &G));
346: PetscCall(VecDuplicate(U, &Gp));
348: PetscCall(FormFunctionGradient(tao, U, &f, G, ctx));
350: PetscCall(VecCopy(U, Up));
351: PetscCall(VecGetArray(Up, &u));
352: u[0] += eps;
353: PetscCall(VecRestoreArray(Up, &u));
354: PetscCall(FormFunctionGradient(tao, Up, &f, Gp, ctx));
355: PetscCall(VecAXPY(Gp, -1, G));
356: PetscCall(VecScale(Gp, 1. / eps));
357: PetscCall(VecGetArray(Gp, &u));
358: arr[0] = u[0];
359: arr[1] = u[1];
360: PetscCall(VecRestoreArray(Gp, &u));
362: PetscCall(VecCopy(U, Up));
363: PetscCall(VecGetArray(Up, &u));
364: u[1] += eps;
365: PetscCall(VecRestoreArray(Up, &u));
366: PetscCall(FormFunctionGradient(tao, Up, &f, Gp, ctx));
367: PetscCall(VecAXPY(Gp, -1, G));
368: PetscCall(VecScale(Gp, 1. / eps));
369: PetscCall(VecGetArray(Gp, &u));
370: arr[2] = u[0];
371: arr[3] = u[1];
372: PetscCall(VecRestoreArray(Gp, &u));
374: PetscCall(VecDestroy(&G));
375: PetscCall(VecDestroy(&Gp));
376: PetscCall(VecDestroy(&Up));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
381: {
382: User user_ptr;
383: PetscScalar *y_ptr;
385: PetscFunctionBeginUser;
386: PetscCall(MatShellGetContext(mat, &user_ptr));
387: PetscCall(VecCopy(svec, user_ptr->Dir));
388: PetscCall(VecGetArray(y, &y_ptr));
389: PetscCall(Adjoint2(user_ptr->U, y_ptr, user_ptr));
390: PetscCall(VecRestoreArray(y, &y_ptr));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: int main(int argc, char **argv)
395: {
396: PetscBool monitor = PETSC_FALSE, mf = PETSC_TRUE;
397: PetscInt mode = 0;
398: PetscMPIInt size;
399: struct _n_User user;
400: Vec x; /* working vector for TAO */
401: PetscScalar *x_ptr, arr[4];
402: PetscScalar ic1 = 2.2, ic2 = -0.7; /* initial guess for TAO */
403: Tao tao;
404: KSP ksp;
405: PC pc;
407: /* Initialize program */
408: PetscFunctionBeginUser;
409: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
410: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
411: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
413: /* Set runtime options */
414: user.next_output = 0.0;
415: user.mu = 1.0e3;
416: user.steps = 0;
417: user.ftime = 0.5;
418: user.implicitform = PETSC_TRUE;
420: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
421: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
422: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mode", &mode, NULL));
423: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ic1", &ic1, NULL));
424: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ic2", &ic2, NULL));
425: PetscCall(PetscOptionsGetBool(NULL, NULL, "-my_tao_mf", &mf, NULL)); /* matrix-free hessian for optimization */
426: PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &user.implicitform, NULL));
428: /* Create necessary matrix and vectors, solve same ODE on every process */
429: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.A));
430: PetscCall(MatSetSizes(user.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
431: PetscCall(MatSetFromOptions(user.A));
432: PetscCall(MatSetUp(user.A));
433: PetscCall(MatCreateVecs(user.A, &user.U, NULL));
434: PetscCall(MatCreateVecs(user.A, &user.Dir, NULL));
435: PetscCall(MatCreateVecs(user.A, &user.Lambda[0], NULL));
436: PetscCall(MatCreateVecs(user.A, &user.Lambda2[0], NULL));
437: PetscCall(MatCreateVecs(user.A, &user.Ihp1[0], NULL));
439: /* Create timestepping solver context */
440: PetscCall(TSCreate(PETSC_COMM_WORLD, &user.ts));
441: PetscCall(TSSetEquationType(user.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
442: if (user.implicitform) {
443: PetscCall(TSSetIFunction(user.ts, NULL, IFunction, &user));
444: PetscCall(TSSetIJacobian(user.ts, user.A, user.A, IJacobian, &user));
445: PetscCall(TSSetType(user.ts, TSCN));
446: } else {
447: PetscCall(TSSetRHSFunction(user.ts, NULL, RHSFunction, &user));
448: PetscCall(TSSetRHSJacobian(user.ts, user.A, user.A, RHSJacobian, &user));
449: PetscCall(TSSetType(user.ts, TSRK));
450: }
451: PetscCall(TSSetMaxTime(user.ts, user.ftime));
452: PetscCall(TSSetExactFinalTime(user.ts, TS_EXACTFINALTIME_MATCHSTEP));
454: if (monitor) PetscCall(TSMonitorSet(user.ts, Monitor, &user, NULL));
456: /* Set ODE initial conditions */
457: PetscCall(VecGetArray(user.U, &x_ptr));
458: x_ptr[0] = 2.0;
459: x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
460: PetscCall(VecRestoreArray(user.U, &x_ptr));
462: /* Set runtime options */
463: PetscCall(TSSetFromOptions(user.ts));
465: /* Obtain the observation */
466: PetscCall(TSSolve(user.ts, user.U));
467: PetscCall(VecGetArray(user.U, &x_ptr));
468: user.ob[0] = x_ptr[0];
469: user.ob[1] = x_ptr[1];
470: PetscCall(VecRestoreArray(user.U, &x_ptr));
472: PetscCall(VecDuplicate(user.U, &x));
473: PetscCall(VecGetArray(x, &x_ptr));
474: x_ptr[0] = ic1;
475: x_ptr[1] = ic2;
476: PetscCall(VecRestoreArray(x, &x_ptr));
478: /* Save trajectory of solution so that TSAdjointSolve() may be used */
479: PetscCall(TSSetSaveTrajectory(user.ts));
481: /* Compare finite difference and second-order adjoint. */
482: switch (mode) {
483: case 2:
484: PetscCall(FiniteDiff(x, arr, &user));
485: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Finite difference approximation of the Hessian\n"));
486: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g %g\n%g %g\n", (double)arr[0], (double)arr[1], (double)arr[2], (double)arr[3]));
487: break;
488: case 1: /* Compute the Hessian column by column */
489: PetscCall(VecCopy(x, user.U));
490: PetscCall(VecGetArray(user.Dir, &x_ptr));
491: x_ptr[0] = 1.;
492: x_ptr[1] = 0.;
493: PetscCall(VecRestoreArray(user.Dir, &x_ptr));
494: PetscCall(Adjoint2(user.U, arr, &user));
495: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nFirst column of the Hessian\n"));
496: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g\n%g\n", (double)arr[0], (double)arr[1]));
497: PetscCall(VecCopy(x, user.U));
498: PetscCall(VecGetArray(user.Dir, &x_ptr));
499: x_ptr[0] = 0.;
500: x_ptr[1] = 1.;
501: PetscCall(VecRestoreArray(user.Dir, &x_ptr));
502: PetscCall(Adjoint2(user.U, arr, &user));
503: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSecond column of the Hessian\n"));
504: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g\n%g\n", (double)arr[0], (double)arr[1]));
505: break;
506: case 0:
507: /* Create optimization context and set up */
508: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
509: PetscCall(TaoSetType(tao, TAOBLMVM));
510: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
512: if (mf) {
513: PetscCall(MatCreateShell(PETSC_COMM_SELF, 2, 2, 2, 2, (void *)&user, &user.H));
514: PetscCall(MatShellSetOperation(user.H, MATOP_MULT, (void (*)(void))HessianProductMat));
515: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
516: PetscCall(TaoSetHessian(tao, user.H, user.H, MatrixFreeHessian, (void *)&user));
517: } else { /* Create Hessian matrix */
518: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.H));
519: PetscCall(MatSetSizes(user.H, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
520: PetscCall(MatSetUp(user.H));
521: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
522: }
524: /* Not use any preconditioner */
525: PetscCall(TaoGetKSP(tao, &ksp));
526: if (ksp) {
527: PetscCall(KSPGetPC(ksp, &pc));
528: PetscCall(PCSetType(pc, PCNONE));
529: }
531: /* Set initial solution guess */
532: PetscCall(TaoSetSolution(tao, x));
533: PetscCall(TaoSetFromOptions(tao));
534: PetscCall(TaoSolve(tao));
535: PetscCall(TaoDestroy(&tao));
536: PetscCall(MatDestroy(&user.H));
537: break;
538: default:
539: break;
540: }
542: /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */
543: PetscCall(MatDestroy(&user.A));
544: PetscCall(VecDestroy(&user.U));
545: PetscCall(VecDestroy(&user.Lambda[0]));
546: PetscCall(VecDestroy(&user.Lambda2[0]));
547: PetscCall(VecDestroy(&user.Ihp1[0]));
548: PetscCall(VecDestroy(&user.Dir));
549: PetscCall(TSDestroy(&user.ts));
550: PetscCall(VecDestroy(&x));
551: PetscCall(PetscFinalize());
552: return 0;
553: }
555: /*TEST
556: build:
557: requires: !complex !single
559: test:
560: args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125
561: output_file: output/ex20opt_ic_1.out
563: test:
564: suffix: 2
565: args: -ts_type beuler -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none
566: output_file: output/ex20opt_ic_2.out
568: test:
569: suffix: 3
570: args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none
571: output_file: output/ex20opt_ic_3.out
573: test:
574: suffix: 4
575: args: -implicitform 0 -ts_dt 0.01 -ts_max_steps 2 -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view -mode 1 -my_tao_mf
576: TEST*/