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*/