Actual source code: ex20opt_p.c

      2: static char help[] = "Solves the van der Pol equation.\n\
      3: Input parameters include:\n";

      5: /*
      6:    Concepts: TS^time-dependent nonlinear problems
      7:    Concepts: TS^van der Pol equation DAE equivalent
      8:    Concepts: TS^Optimization using adjoint sensitivity analysis
      9:    Processors: 1
     10: */
     11: /* ------------------------------------------------------------------------

     13:   Notes:
     14:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
     15:   The nonlinear problem is written in a DAE equivalent form.
     16:   The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
     17:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
     18:   ------------------------------------------------------------------------- */
     19: #include <petsctao.h>
     20: #include <petscts.h>

     22: typedef struct _n_User *User;
     23: struct _n_User {
     24:   TS        ts;
     25:   PetscReal mu;
     26:   PetscReal next_output;

     28:   /* Sensitivity analysis support */
     29:   PetscReal ftime;
     30:   Mat       A;                       /* Jacobian matrix */
     31:   Mat       Jacp;                    /* JacobianP matrix */
     32:   Mat       H;                       /* Hessian matrix for optimization */
     33:   Vec       U,Lambda[1],Mup[1];      /* adjoint variables */
     34:   Vec       Lambda2[1],Mup2[1];      /* second-order adjoint variables */
     35:   Vec       Ihp1[1];                 /* working space for Hessian evaluations */
     36:   Vec       Ihp2[1];                 /* working space for Hessian evaluations */
     37:   Vec       Ihp3[1];                 /* working space for Hessian evaluations */
     38:   Vec       Ihp4[1];                 /* working space for Hessian evaluations */
     39:   Vec       Dir;                     /* direction vector */
     40:   PetscReal ob[2];                   /* observation used by the cost function */
     41:   PetscBool implicitform;            /* implicit ODE? */
     42: };

     44: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
     45: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
     46: PetscErrorCode Adjoint2(Vec,PetscScalar[],User);

     48: /* ----------------------- Explicit form of the ODE  -------------------- */

     50: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
     51: {
     52:   PetscErrorCode    ierr;
     53:   User              user = (User)ctx;
     54:   PetscScalar       *f;
     55:   const PetscScalar *u;

     58:   VecGetArrayRead(U,&u);
     59:   VecGetArray(F,&f);
     60:   f[0] = u[1];
     61:   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
     62:   VecRestoreArrayRead(U,&u);
     63:   VecRestoreArray(F,&f);
     64:   return(0);
     65: }

     67: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
     68: {
     69:   PetscErrorCode    ierr;
     70:   User              user = (User)ctx;
     71:   PetscReal         mu   = user->mu;
     72:   PetscInt          rowcol[] = {0,1};
     73:   PetscScalar       J[2][2];
     74:   const PetscScalar *u;

     77:   VecGetArrayRead(U,&u);

     79:   J[0][0] = 0;
     80:   J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
     81:   J[0][1] = 1.0;
     82:   J[1][1] = mu*(1.0-u[0]*u[0]);
     83:   MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
     84:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     85:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     86:   if (B && A != B) {
     87:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
     88:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
     89:   }

     91:   VecRestoreArrayRead(U,&u);
     92:   return(0);
     93: }

     95: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
     96: {
     97:   const PetscScalar *vl,*vr,*u;
     98:   PetscScalar       *vhv;
     99:   PetscScalar       dJdU[2][2][2]={{{0}}};
    100:   PetscInt          i,j,k;
    101:   User              user = (User)ctx;
    102:   PetscErrorCode    ierr;

    105:   VecGetArrayRead(U,&u);
    106:   VecGetArrayRead(Vl[0],&vl);
    107:   VecGetArrayRead(Vr,&vr);
    108:   VecGetArray(VHV[0],&vhv);

    110:   dJdU[1][0][0] = -2.*user->mu*u[1];
    111:   dJdU[1][1][0] = -2.*user->mu*u[0];
    112:   dJdU[1][0][1] = -2.*user->mu*u[0];
    113:   for (j=0; j<2; j++) {
    114:     vhv[j] = 0;
    115:     for (k=0; k<2; k++)
    116:       for (i=0; i<2; i++)
    117:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
    118:   }

    120:   VecRestoreArrayRead(U,&u);
    121:   VecRestoreArrayRead(Vl[0],&vl);
    122:   VecRestoreArrayRead(Vr,&vr);
    123:   VecRestoreArray(VHV[0],&vhv);
    124:   return(0);
    125: }

    127: static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    128: {
    129:   const PetscScalar *vl,*vr,*u;
    130:   PetscScalar       *vhv;
    131:   PetscScalar       dJdP[2][2][1]={{{0}}};
    132:   PetscInt          i,j,k;
    133:   PetscErrorCode    ierr;

    136:   VecGetArrayRead(U,&u);
    137:   VecGetArrayRead(Vl[0],&vl);
    138:   VecGetArrayRead(Vr,&vr);
    139:   VecGetArray(VHV[0],&vhv);

    141:   dJdP[1][0][0] = -(1.+2.*u[0]*u[1]);
    142:   dJdP[1][1][0] = 1.-u[0]*u[0];
    143:   for (j=0; j<2; j++) {
    144:     vhv[j] = 0;
    145:     for (k=0; k<1; k++)
    146:       for (i=0; i<2; i++)
    147:         vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
    148:   }

    150:   VecRestoreArrayRead(U,&u);
    151:   VecRestoreArrayRead(Vl[0],&vl);
    152:   VecRestoreArrayRead(Vr,&vr);
    153:   VecRestoreArray(VHV[0],&vhv);
    154:   return(0);
    155: }

    157: static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    158: {
    159:   const PetscScalar *vl,*vr,*u;
    160:   PetscScalar       *vhv;
    161:   PetscScalar       dJdU[2][1][2]={{{0}}};
    162:   PetscInt          i,j,k;
    163:   PetscErrorCode    ierr;

    166:   VecGetArrayRead(U,&u);
    167:   VecGetArrayRead(Vl[0],&vl);
    168:   VecGetArrayRead(Vr,&vr);
    169:   VecGetArray(VHV[0],&vhv);

    171:   dJdU[1][0][0] = -1.-2.*u[1]*u[0];
    172:   dJdU[1][0][1] = 1.-u[0]*u[0];
    173:   for (j=0; j<1; j++) {
    174:     vhv[j] = 0;
    175:     for (k=0; k<2; k++)
    176:       for (i=0; i<2; i++)
    177:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
    178:   }

    180:   VecRestoreArrayRead(U,&u);
    181:   VecRestoreArrayRead(Vl[0],&vl);
    182:   VecRestoreArrayRead(Vr,&vr);
    183:   VecRestoreArray(VHV[0],&vhv);
    184:   return(0);
    185: }

    187: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    188: {
    190:   return(0);
    191: }

    193: /* ----------------------- Implicit form of the ODE  -------------------- */

    195: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
    196: {
    197:   PetscErrorCode    ierr;
    198:   User              user = (User)ctx;
    199:   PetscScalar       *f;
    200:   const PetscScalar *u,*udot;

    203:   VecGetArrayRead(U,&u);
    204:   VecGetArrayRead(Udot,&udot);
    205:   VecGetArray(F,&f);

    207:   f[0] = udot[0] - u[1];
    208:   f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;

    210:   VecRestoreArrayRead(U,&u);
    211:   VecRestoreArrayRead(Udot,&udot);
    212:   VecRestoreArray(F,&f);
    213:   return(0);
    214: }

    216: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
    217: {
    218:   PetscErrorCode    ierr;
    219:   User              user     = (User)ctx;
    220:   PetscInt          rowcol[] = {0,1};
    221:   PetscScalar       J[2][2];
    222:   const PetscScalar *u;

    225:   VecGetArrayRead(U,&u);

    227:   J[0][0] = a;     J[0][1] =  -1.0;
    228:   J[1][0] = user->mu*(1.0 + 2.0*u[0]*u[1]);   J[1][1] = a - user->mu*(1.0-u[0]*u[0]);
    229:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
    230:   VecRestoreArrayRead(U,&u);
    231:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    232:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    233:   if (A != B) {
    234:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    235:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    236:   }

    238:   VecRestoreArrayRead(U,&u);
    239:   return(0);
    240: }

    242: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
    243: {
    244:   PetscErrorCode    ierr;
    245:   PetscInt          row[] = {0,1},col[]={0};
    246:   PetscScalar       J[2][1];
    247:   const PetscScalar *u;

    250:   VecGetArrayRead(U,&u);

    252:   J[0][0] = 0;
    253:   J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
    254:   MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
    255:   VecRestoreArrayRead(U,&u);
    256:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    257:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

    259:   VecRestoreArrayRead(U,&u);
    260:   return(0);
    261: }

    263: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    264: {
    265:   const PetscScalar *vl,*vr,*u;
    266:   PetscScalar       *vhv;
    267:   PetscScalar       dJdU[2][2][2]={{{0}}};
    268:   PetscInt          i,j,k;
    269:   User              user = (User)ctx;
    270:   PetscErrorCode    ierr;

    273:   VecGetArrayRead(U,&u);
    274:   VecGetArrayRead(Vl[0],&vl);
    275:   VecGetArrayRead(Vr,&vr);
    276:   VecGetArray(VHV[0],&vhv);

    278:   dJdU[1][0][0] = 2.*user->mu*u[1];
    279:   dJdU[1][1][0] = 2.*user->mu*u[0];
    280:   dJdU[1][0][1] = 2.*user->mu*u[0];
    281:   for (j=0; j<2; j++) {
    282:     vhv[j] = 0;
    283:     for (k=0; k<2; k++)
    284:       for (i=0; i<2; i++)
    285:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
    286:   }

    288:   VecRestoreArrayRead(U,&u);
    289:   VecRestoreArrayRead(Vl[0],&vl);
    290:   VecRestoreArrayRead(Vr,&vr);
    291:   VecRestoreArray(VHV[0],&vhv);
    292:   return(0);
    293: }

    295: static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    296: {
    297:   const PetscScalar *vl,*vr,*u;
    298:   PetscScalar       *vhv;
    299:   PetscScalar       dJdP[2][2][1]={{{0}}};
    300:   PetscInt          i,j,k;
    301:   PetscErrorCode    ierr;

    304:   VecGetArrayRead(U,&u);
    305:   VecGetArrayRead(Vl[0],&vl);
    306:   VecGetArrayRead(Vr,&vr);
    307:   VecGetArray(VHV[0],&vhv);

    309:   dJdP[1][0][0] = 1.+2.*u[0]*u[1];
    310:   dJdP[1][1][0] = u[0]*u[0]-1.;
    311:   for (j=0; j<2; j++) {
    312:     vhv[j] = 0;
    313:     for (k=0; k<1; k++)
    314:       for (i=0; i<2; i++)
    315:         vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
    316:   }

    318:   VecRestoreArrayRead(U,&u);
    319:   VecRestoreArrayRead(Vl[0],&vl);
    320:   VecRestoreArrayRead(Vr,&vr);
    321:   VecRestoreArray(VHV[0],&vhv);
    322:   return(0);
    323: }

    325: static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    326: {
    327:   const PetscScalar *vl,*vr,*u;
    328:   PetscScalar       *vhv;
    329:   PetscScalar       dJdU[2][1][2]={{{0}}};
    330:   PetscInt          i,j,k;
    331:   PetscErrorCode    ierr;

    334:   VecGetArrayRead(U,&u);
    335:   VecGetArrayRead(Vl[0],&vl);
    336:   VecGetArrayRead(Vr,&vr);
    337:   VecGetArray(VHV[0],&vhv);

    339:   dJdU[1][0][0] = 1.+2.*u[1]*u[0];
    340:   dJdU[1][0][1] = u[0]*u[0]-1.;
    341:   for (j=0; j<1; j++) {
    342:     vhv[j] = 0;
    343:     for (k=0; k<2; k++)
    344:       for (i=0; i<2; i++)
    345:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
    346:   }

    348:   VecRestoreArrayRead(U,&u);
    349:   VecRestoreArrayRead(Vl[0],&vl);
    350:   VecRestoreArrayRead(Vr,&vr);
    351:   VecRestoreArray(VHV[0],&vhv);
    352:   return(0);
    353: }

    355: static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
    356: {
    358:   return(0);
    359: }

    361: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
    362: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
    363: {
    364:   PetscErrorCode    ierr;
    365:   const PetscScalar *x;
    366:   PetscReal         tfinal, dt;
    367:   User              user = (User)ctx;
    368:   Vec               interpolatedX;

    371:   TSGetTimeStep(ts,&dt);
    372:   TSGetMaxTime(ts,&tfinal);

    374:   while (user->next_output <= t && user->next_output <= tfinal) {
    375:     VecDuplicate(X,&interpolatedX);
    376:     TSInterpolate(ts,user->next_output,interpolatedX);
    377:     VecGetArrayRead(interpolatedX,&x);
    378:     PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
    379:                        (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
    380:                        (double)PetscRealPart(x[1]));
    381:     VecRestoreArrayRead(interpolatedX,&x);
    382:     VecDestroy(&interpolatedX);
    383:     user->next_output += PetscRealConstant(0.1);
    384:   }
    385:   return(0);
    386: }

    388: int main(int argc,char **argv)
    389: {
    390:   Vec                P;
    391:   PetscBool          monitor = PETSC_FALSE;
    392:   PetscScalar        *x_ptr;
    393:   const PetscScalar  *y_ptr;
    394:   PetscMPIInt        size;
    395:   struct _n_User     user;
    396:   PetscErrorCode     ierr;
    397:   Tao                tao;
    398:   KSP                ksp;
    399:   PC                 pc;

    401:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    402:      Initialize program
    403:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    404:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
    405:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    406:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");

    408:   /* Create TAO solver and set desired solution method */
    409:   TaoCreate(PETSC_COMM_WORLD,&tao);
    410:   TaoSetType(tao,TAOBQNLS);

    412:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    413:     Set runtime options
    414:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    415:   user.next_output  = 0.0;
    416:   user.mu           = PetscRealConstant(1.0e3);
    417:   user.ftime        = PetscRealConstant(0.5);
    418:   user.implicitform = PETSC_TRUE;

    420:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
    421:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
    422:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);

    424:   /* Create necessary matrix and vectors, solve same ODE on every process */
    425:   MatCreate(PETSC_COMM_WORLD,&user.A);
    426:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
    427:   MatSetFromOptions(user.A);
    428:   MatSetUp(user.A);
    429:   MatCreateVecs(user.A,&user.U,NULL);
    430:   MatCreateVecs(user.A,&user.Lambda[0],NULL);
    431:   MatCreateVecs(user.A,&user.Lambda2[0],NULL);
    432:   MatCreateVecs(user.A,&user.Ihp1[0],NULL);
    433:   MatCreateVecs(user.A,&user.Ihp2[0],NULL);

    435:   MatCreate(PETSC_COMM_WORLD,&user.Jacp);
    436:   MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
    437:   MatSetFromOptions(user.Jacp);
    438:   MatSetUp(user.Jacp);
    439:   MatCreateVecs(user.Jacp,&user.Dir,NULL);
    440:   MatCreateVecs(user.Jacp,&user.Mup[0],NULL);
    441:   MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);
    442:   MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);
    443:   MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);

    445:   /* Create timestepping solver context */
    446:   TSCreate(PETSC_COMM_WORLD,&user.ts);
    447:   TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
    448:   if (user.implicitform) {
    449:     TSSetIFunction(user.ts,NULL,IFunction,&user);
    450:     TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
    451:     TSSetType(user.ts,TSCN);
    452:   } else {
    453:     TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
    454:     TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
    455:     TSSetType(user.ts,TSRK);
    456:   }
    457:   TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
    458:   TSSetMaxTime(user.ts,user.ftime);
    459:   TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);

    461:   if (monitor) {
    462:     TSMonitorSet(user.ts,Monitor,&user,NULL);
    463:   }

    465:   /* Set ODE initial conditions */
    466:   VecGetArray(user.U,&x_ptr);
    467:   x_ptr[0] = 2.0;
    468:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
    469:   VecRestoreArray(user.U,&x_ptr);
    470:   TSSetTimeStep(user.ts,PetscRealConstant(0.001));

    472:   /* Set runtime options */
    473:   TSSetFromOptions(user.ts);

    475:   TSSolve(user.ts,user.U);
    476:   VecGetArrayRead(user.U,&y_ptr);
    477:   user.ob[0] = y_ptr[0];
    478:   user.ob[1] = y_ptr[1];
    479:   VecRestoreArrayRead(user.U,&y_ptr);

    481:   /* Save trajectory of solution so that TSAdjointSolve() may be used.
    482:      Skip checkpointing for the first TSSolve since no adjoint run follows it.
    483:    */
    484:   TSSetSaveTrajectory(user.ts);

    486:   /* Optimization starts */
    487:   MatCreate(PETSC_COMM_WORLD,&user.H);
    488:   MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);
    489:   MatSetUp(user.H); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */

    491:   /* Set initial solution guess */
    492:   MatCreateVecs(user.Jacp,&P,NULL);
    493:   VecGetArray(P,&x_ptr);
    494:   x_ptr[0] = PetscRealConstant(1.2);
    495:   VecRestoreArray(P,&x_ptr);
    496:   TaoSetInitialVector(tao,P);

    498:   /* Set routine for function and gradient evaluation */
    499:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
    500:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);

    502:   /* Check for any TAO command line options */
    503:   TaoGetKSP(tao,&ksp);
    504:   if (ksp) {
    505:     KSPGetPC(ksp,&pc);
    506:     PCSetType(pc,PCNONE);
    507:   }
    508:   TaoSetFromOptions(tao);

    510:   TaoSolve(tao);

    512:   VecView(P,PETSC_VIEWER_STDOUT_WORLD);
    513:   /* Free TAO data structures */
    514:   TaoDestroy(&tao);

    516:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    517:      Free work space.  All PETSc objects should be destroyed when they
    518:      are no longer needed.
    519:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    520:   MatDestroy(&user.H);
    521:   MatDestroy(&user.A);
    522:   VecDestroy(&user.U);
    523:   MatDestroy(&user.Jacp);
    524:   VecDestroy(&user.Lambda[0]);
    525:   VecDestroy(&user.Mup[0]);
    526:   VecDestroy(&user.Lambda2[0]);
    527:   VecDestroy(&user.Mup2[0]);
    528:   VecDestroy(&user.Ihp1[0]);
    529:   VecDestroy(&user.Ihp2[0]);
    530:   VecDestroy(&user.Ihp3[0]);
    531:   VecDestroy(&user.Ihp4[0]);
    532:   VecDestroy(&user.Dir);
    533:   TSDestroy(&user.ts);
    534:   VecDestroy(&P);
    535:   PetscFinalize();
    536:   return ierr;
    537: }

    539: /* ------------------------------------------------------------------ */
    540: /*
    541:    FormFunctionGradient - Evaluates the function and corresponding gradient.

    543:    Input Parameters:
    544:    tao - the Tao context
    545:    X   - the input vector
    546:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

    548:    Output Parameters:
    549:    f   - the newly evaluated function
    550:    G   - the newly evaluated gradient
    551: */
    552: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
    553: {
    554:   User              user_ptr = (User)ctx;
    555:   TS                ts = user_ptr->ts;
    556:   PetscScalar       *x_ptr,*g;
    557:   const PetscScalar *y_ptr;
    558:   PetscErrorCode    ierr;

    561:   VecGetArrayRead(P,&y_ptr);
    562:   user_ptr->mu = y_ptr[0];
    563:   VecRestoreArrayRead(P,&y_ptr);

    565:   TSSetTime(ts,0.0);
    566:   TSSetStepNumber(ts,0);
    567:   TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
    568:   TSSetFromOptions(ts);
    569:   VecGetArray(user_ptr->U,&x_ptr);
    570:   x_ptr[0] = 2.0;
    571:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
    572:   VecRestoreArray(user_ptr->U,&x_ptr);

    574:   TSSolve(ts,user_ptr->U);

    576:   VecGetArrayRead(user_ptr->U,&y_ptr);
    577:   *f   = (y_ptr[0]-user_ptr->ob[0])*(y_ptr[0]-user_ptr->ob[0])+(y_ptr[1]-user_ptr->ob[1])*(y_ptr[1]-user_ptr->ob[1]);

    579:   /*   Reset initial conditions for the adjoint integration */
    580:   VecGetArray(user_ptr->Lambda[0],&x_ptr);
    581:   x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
    582:   x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
    583:   VecRestoreArrayRead(user_ptr->U,&y_ptr);
    584:   VecRestoreArray(user_ptr->Lambda[0],&x_ptr);

    586:   VecGetArray(user_ptr->Mup[0],&x_ptr);
    587:   x_ptr[0] = 0.0;
    588:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
    589:   TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);

    591:   TSAdjointSolve(ts);

    593:   VecGetArray(user_ptr->Mup[0],&x_ptr);
    594:   VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
    595:   VecGetArray(G,&g);
    596:   g[0] = y_ptr[1]*(-10.0/(81.0*user_ptr->mu*user_ptr->mu)+2.0*292.0/(2187.0*user_ptr->mu*user_ptr->mu*user_ptr->mu))+x_ptr[0];
    597:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
    598:   VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
    599:   VecRestoreArray(G,&g);
    600:   return(0);
    601: }

    603: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
    604: {
    605:   User           user_ptr = (User)ctx;
    606:   PetscScalar    harr[1];
    607:   const PetscInt rows[1] = {0};
    608:   PetscInt       col = 0;

    612:   Adjoint2(P,harr,user_ptr);
    613:   MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);

    615:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
    616:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
    617:   if (H != Hpre) {
    618:     MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
    619:     MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
    620:   }
    621:   return(0);
    622: }

    624: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
    625: {
    626:   TS                ts = ctx->ts;
    627:   const PetscScalar *z_ptr;
    628:   PetscScalar       *x_ptr,*y_ptr,dzdp,dzdp2;
    629:   Mat               tlmsen;
    630:   PetscErrorCode    ierr;

    633:   /* Reset TSAdjoint so that AdjointSetUp will be called again */
    634:   TSAdjointReset(ts);

    636:   /* The directional vector should be 1 since it is one-dimensional */
    637:   VecGetArray(ctx->Dir,&x_ptr);
    638:   x_ptr[0] = 1.;
    639:   VecRestoreArray(ctx->Dir,&x_ptr);

    641:   VecGetArrayRead(P,&z_ptr);
    642:   ctx->mu = z_ptr[0];
    643:   VecRestoreArrayRead(P,&z_ptr);

    645:   dzdp  = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
    646:   dzdp2 = 2.*10.0/(81.0*ctx->mu*ctx->mu*ctx->mu) - 3.0*2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu*ctx->mu);

    648:   TSSetTime(ts,0.0);
    649:   TSSetStepNumber(ts,0);
    650:   TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
    651:   TSSetFromOptions(ts);
    652:   TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);

    654:   MatZeroEntries(ctx->Jacp);
    655:   MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
    656:   MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
    657:   MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);

    659:   TSAdjointSetForward(ts,ctx->Jacp);
    660:   VecGetArray(ctx->U,&y_ptr);
    661:   y_ptr[0] = 2.0;
    662:   y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
    663:   VecRestoreArray(ctx->U,&y_ptr);
    664:   TSSolve(ts,ctx->U);

    666:   /* Set terminal conditions for first- and second-order adjonts */
    667:   VecGetArrayRead(ctx->U,&z_ptr);
    668:   VecGetArray(ctx->Lambda[0],&y_ptr);
    669:   y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
    670:   y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
    671:   VecRestoreArray(ctx->Lambda[0],&y_ptr);
    672:   VecRestoreArrayRead(ctx->U,&z_ptr);
    673:   VecGetArray(ctx->Mup[0],&y_ptr);
    674:   y_ptr[0] = 0.0;
    675:   VecRestoreArray(ctx->Mup[0],&y_ptr);
    676:   TSForwardGetSensitivities(ts,NULL,&tlmsen);
    677:   MatDenseGetColumn(tlmsen,0,&x_ptr);
    678:   VecGetArray(ctx->Lambda2[0],&y_ptr);
    679:   y_ptr[0] = 2.*x_ptr[0];
    680:   y_ptr[1] = 2.*x_ptr[1];
    681:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
    682:   VecGetArray(ctx->Mup2[0],&y_ptr);
    683:   y_ptr[0] = 0.0;
    684:   VecRestoreArray(ctx->Mup2[0],&y_ptr);
    685:   MatDenseRestoreColumn(tlmsen,&x_ptr);
    686:   TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
    687:   if (ctx->implicitform) {
    688:     TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
    689:   } else {
    690:     TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
    691:   }
    692:   TSAdjointSolve(ts);

    694:   VecGetArray(ctx->Lambda[0],&x_ptr);
    695:   VecGetArray(ctx->Lambda2[0],&y_ptr);
    696:   VecGetArrayRead(ctx->Mup2[0],&z_ptr);

    698:   arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];

    700:   VecRestoreArray(ctx->Lambda2[0],&x_ptr);
    701:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
    702:   VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);

    704:   /* Disable second-order adjoint mode */
    705:   TSAdjointReset(ts);
    706:   TSAdjointResetForward(ts);
    707:   return(0);
    708: }

    710: /*TEST
    711:     build:
    712:       requires: !complex !single
    713:     test:
    714:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
    715:       output_file: output/ex20opt_p_1.out

    717:     test:
    718:       suffix: 2
    719:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none
    720:       output_file: output/ex20opt_p_2.out

    722:     test:
    723:       suffix: 3
    724:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
    725:       output_file: output/ex20opt_p_3.out

    727:     test:
    728:       suffix: 4
    729:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none
    730:       output_file: output/ex20opt_p_4.out

    732: TEST*/