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