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:   User              user = (User)ctx;
 53:   PetscScalar       *f;
 54:   const PetscScalar *u;

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

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

 75:   VecGetArrayRead(U,&u);

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

 89:   VecRestoreArrayRead(U,&u);
 90:   return 0;
 91: }

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

102:   VecGetArrayRead(U,&u);
103:   VecGetArrayRead(Vl[0],&vl);
104:   VecGetArrayRead(Vr,&vr);
105:   VecGetArray(VHV[0],&vhv);

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

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

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

132:   VecGetArrayRead(U,&u);
133:   VecGetArrayRead(Vl[0],&vl);
134:   VecGetArrayRead(Vr,&vr);
135:   VecGetArray(VHV[0],&vhv);

137:   dJdP[1][0][0] = -(1.+2.*u[0]*u[1]);
138:   dJdP[1][1][0] = 1.-u[0]*u[0];
139:   for (j=0; j<2; j++) {
140:     vhv[j] = 0;
141:     for (k=0; k<1; k++)
142:       for (i=0; i<2; i++)
143:         vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
144:   }

146:   VecRestoreArrayRead(U,&u);
147:   VecRestoreArrayRead(Vl[0],&vl);
148:   VecRestoreArrayRead(Vr,&vr);
149:   VecRestoreArray(VHV[0],&vhv);
150:   return 0;
151: }

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

161:   VecGetArrayRead(U,&u);
162:   VecGetArrayRead(Vl[0],&vl);
163:   VecGetArrayRead(Vr,&vr);
164:   VecGetArray(VHV[0],&vhv);

166:   dJdU[1][0][0] = -1.-2.*u[1]*u[0];
167:   dJdU[1][0][1] = 1.-u[0]*u[0];
168:   for (j=0; j<1; j++) {
169:     vhv[j] = 0;
170:     for (k=0; k<2; k++)
171:       for (i=0; i<2; i++)
172:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
173:   }

175:   VecRestoreArrayRead(U,&u);
176:   VecRestoreArrayRead(Vl[0],&vl);
177:   VecRestoreArrayRead(Vr,&vr);
178:   VecRestoreArray(VHV[0],&vhv);
179:   return 0;
180: }

182: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
183: {
185:   return 0;
186: }

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

190: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
191: {
192:   User              user = (User)ctx;
193:   PetscScalar       *f;
194:   const PetscScalar *u,*udot;

197:   VecGetArrayRead(U,&u);
198:   VecGetArrayRead(Udot,&udot);
199:   VecGetArray(F,&f);

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

204:   VecRestoreArrayRead(U,&u);
205:   VecRestoreArrayRead(Udot,&udot);
206:   VecRestoreArray(F,&f);
207:   return 0;
208: }

210: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
211: {
212:   User              user     = (User)ctx;
213:   PetscInt          rowcol[] = {0,1};
214:   PetscScalar       J[2][2];
215:   const PetscScalar *u;

218:   VecGetArrayRead(U,&u);

220:   J[0][0] = a;     J[0][1] =  -1.0;
221:   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]);
222:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
223:   VecRestoreArrayRead(U,&u);
224:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
225:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
226:   if (A != B) {
227:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
228:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
229:   }

231:   VecRestoreArrayRead(U,&u);
232:   return 0;
233: }

235: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
236: {
237:   PetscInt          row[] = {0,1},col[]={0};
238:   PetscScalar       J[2][1];
239:   const PetscScalar *u;

242:   VecGetArrayRead(U,&u);

244:   J[0][0] = 0;
245:   J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
246:   MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
247:   VecRestoreArrayRead(U,&u);
248:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
249:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

251:   VecRestoreArrayRead(U,&u);
252:   return 0;
253: }

255: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
256: {
257:   const PetscScalar *vl,*vr,*u;
258:   PetscScalar       *vhv;
259:   PetscScalar       dJdU[2][2][2]={{{0}}};
260:   PetscInt          i,j,k;
261:   User              user = (User)ctx;

264:   VecGetArrayRead(U,&u);
265:   VecGetArrayRead(Vl[0],&vl);
266:   VecGetArrayRead(Vr,&vr);
267:   VecGetArray(VHV[0],&vhv);

269:   dJdU[1][0][0] = 2.*user->mu*u[1];
270:   dJdU[1][1][0] = 2.*user->mu*u[0];
271:   dJdU[1][0][1] = 2.*user->mu*u[0];
272:   for (j=0; j<2; j++) {
273:     vhv[j] = 0;
274:     for (k=0; k<2; k++)
275:       for (i=0; i<2; i++)
276:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
277:   }

279:   VecRestoreArrayRead(U,&u);
280:   VecRestoreArrayRead(Vl[0],&vl);
281:   VecRestoreArrayRead(Vr,&vr);
282:   VecRestoreArray(VHV[0],&vhv);
283:   return 0;
284: }

286: static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
287: {
288:   const PetscScalar *vl,*vr,*u;
289:   PetscScalar       *vhv;
290:   PetscScalar       dJdP[2][2][1]={{{0}}};
291:   PetscInt          i,j,k;

294:   VecGetArrayRead(U,&u);
295:   VecGetArrayRead(Vl[0],&vl);
296:   VecGetArrayRead(Vr,&vr);
297:   VecGetArray(VHV[0],&vhv);

299:   dJdP[1][0][0] = 1.+2.*u[0]*u[1];
300:   dJdP[1][1][0] = u[0]*u[0]-1.;
301:   for (j=0; j<2; j++) {
302:     vhv[j] = 0;
303:     for (k=0; k<1; k++)
304:       for (i=0; i<2; i++)
305:         vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
306:   }

308:   VecRestoreArrayRead(U,&u);
309:   VecRestoreArrayRead(Vl[0],&vl);
310:   VecRestoreArrayRead(Vr,&vr);
311:   VecRestoreArray(VHV[0],&vhv);
312:   return 0;
313: }

315: static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
316: {
317:   const PetscScalar *vl,*vr,*u;
318:   PetscScalar       *vhv;
319:   PetscScalar       dJdU[2][1][2]={{{0}}};
320:   PetscInt          i,j,k;

323:   VecGetArrayRead(U,&u);
324:   VecGetArrayRead(Vl[0],&vl);
325:   VecGetArrayRead(Vr,&vr);
326:   VecGetArray(VHV[0],&vhv);

328:   dJdU[1][0][0] = 1.+2.*u[1]*u[0];
329:   dJdU[1][0][1] = u[0]*u[0]-1.;
330:   for (j=0; j<1; j++) {
331:     vhv[j] = 0;
332:     for (k=0; k<2; k++)
333:       for (i=0; i<2; i++)
334:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
335:   }

337:   VecRestoreArrayRead(U,&u);
338:   VecRestoreArrayRead(Vl[0],&vl);
339:   VecRestoreArrayRead(Vr,&vr);
340:   VecRestoreArray(VHV[0],&vhv);
341:   return 0;
342: }

344: static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
345: {
347:   return 0;
348: }

350: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
351: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
352: {
353:   PetscErrorCode    ierr;
354:   const PetscScalar *x;
355:   PetscReal         tfinal, dt;
356:   User              user = (User)ctx;
357:   Vec               interpolatedX;

360:   TSGetTimeStep(ts,&dt);
361:   TSGetMaxTime(ts,&tfinal);

363:   while (user->next_output <= t && user->next_output <= tfinal) {
364:     VecDuplicate(X,&interpolatedX);
365:     TSInterpolate(ts,user->next_output,interpolatedX);
366:     VecGetArrayRead(interpolatedX,&x);
367:     PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
368:                        (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
369:                        (double)PetscRealPart(x[1]));
370:     VecRestoreArrayRead(interpolatedX,&x);
371:     VecDestroy(&interpolatedX);
372:     user->next_output += PetscRealConstant(0.1);
373:   }
374:   return 0;
375: }

377: int main(int argc,char **argv)
378: {
379:   Vec                P;
380:   PetscBool          monitor = PETSC_FALSE;
381:   PetscScalar        *x_ptr;
382:   const PetscScalar  *y_ptr;
383:   PetscMPIInt        size;
384:   struct _n_User     user;
385:   Tao                tao;
386:   KSP                ksp;
387:   PC                 pc;

389:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
390:      Initialize program
391:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
392:   PetscInitialize(&argc,&argv,NULL,help);
393:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

396:   /* Create TAO solver and set desired solution method */
397:   TaoCreate(PETSC_COMM_WORLD,&tao);
398:   TaoSetType(tao,TAOBQNLS);

400:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
401:     Set runtime options
402:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
403:   user.next_output  = 0.0;
404:   user.mu           = PetscRealConstant(1.0e3);
405:   user.ftime        = PetscRealConstant(0.5);
406:   user.implicitform = PETSC_TRUE;

408:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
409:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
410:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);

412:   /* Create necessary matrix and vectors, solve same ODE on every process */
413:   MatCreate(PETSC_COMM_WORLD,&user.A);
414:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
415:   MatSetFromOptions(user.A);
416:   MatSetUp(user.A);
417:   MatCreateVecs(user.A,&user.U,NULL);
418:   MatCreateVecs(user.A,&user.Lambda[0],NULL);
419:   MatCreateVecs(user.A,&user.Lambda2[0],NULL);
420:   MatCreateVecs(user.A,&user.Ihp1[0],NULL);
421:   MatCreateVecs(user.A,&user.Ihp2[0],NULL);

423:   MatCreate(PETSC_COMM_WORLD,&user.Jacp);
424:   MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
425:   MatSetFromOptions(user.Jacp);
426:   MatSetUp(user.Jacp);
427:   MatCreateVecs(user.Jacp,&user.Dir,NULL);
428:   MatCreateVecs(user.Jacp,&user.Mup[0],NULL);
429:   MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);
430:   MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);
431:   MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);

433:   /* Create timestepping solver context */
434:   TSCreate(PETSC_COMM_WORLD,&user.ts);
435:   TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
436:   if (user.implicitform) {
437:     TSSetIFunction(user.ts,NULL,IFunction,&user);
438:     TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
439:     TSSetType(user.ts,TSCN);
440:   } else {
441:     TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
442:     TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
443:     TSSetType(user.ts,TSRK);
444:   }
445:   TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
446:   TSSetMaxTime(user.ts,user.ftime);
447:   TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);

449:   if (monitor) {
450:     TSMonitorSet(user.ts,Monitor,&user,NULL);
451:   }

453:   /* Set ODE initial conditions */
454:   VecGetArray(user.U,&x_ptr);
455:   x_ptr[0] = 2.0;
456:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
457:   VecRestoreArray(user.U,&x_ptr);
458:   TSSetTimeStep(user.ts,PetscRealConstant(0.001));

460:   /* Set runtime options */
461:   TSSetFromOptions(user.ts);

463:   TSSolve(user.ts,user.U);
464:   VecGetArrayRead(user.U,&y_ptr);
465:   user.ob[0] = y_ptr[0];
466:   user.ob[1] = y_ptr[1];
467:   VecRestoreArrayRead(user.U,&y_ptr);

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

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

479:   /* Set initial solution guess */
480:   MatCreateVecs(user.Jacp,&P,NULL);
481:   VecGetArray(P,&x_ptr);
482:   x_ptr[0] = PetscRealConstant(1.2);
483:   VecRestoreArray(P,&x_ptr);
484:   TaoSetSolution(tao,P);

486:   /* Set routine for function and gradient evaluation */
487:   TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
488:   TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);

490:   /* Check for any TAO command line options */
491:   TaoGetKSP(tao,&ksp);
492:   if (ksp) {
493:     KSPGetPC(ksp,&pc);
494:     PCSetType(pc,PCNONE);
495:   }
496:   TaoSetFromOptions(tao);

498:   TaoSolve(tao);

500:   VecView(P,PETSC_VIEWER_STDOUT_WORLD);
501:   /* Free TAO data structures */
502:   TaoDestroy(&tao);

504:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
505:      Free work space.  All PETSc objects should be destroyed when they
506:      are no longer needed.
507:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
508:   MatDestroy(&user.H);
509:   MatDestroy(&user.A);
510:   VecDestroy(&user.U);
511:   MatDestroy(&user.Jacp);
512:   VecDestroy(&user.Lambda[0]);
513:   VecDestroy(&user.Mup[0]);
514:   VecDestroy(&user.Lambda2[0]);
515:   VecDestroy(&user.Mup2[0]);
516:   VecDestroy(&user.Ihp1[0]);
517:   VecDestroy(&user.Ihp2[0]);
518:   VecDestroy(&user.Ihp3[0]);
519:   VecDestroy(&user.Ihp4[0]);
520:   VecDestroy(&user.Dir);
521:   TSDestroy(&user.ts);
522:   VecDestroy(&P);
523:   PetscFinalize();
524:   return 0;
525: }

527: /* ------------------------------------------------------------------ */
528: /*
529:    FormFunctionGradient - Evaluates the function and corresponding gradient.

531:    Input Parameters:
532:    tao - the Tao context
533:    X   - the input vector
534:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

536:    Output Parameters:
537:    f   - the newly evaluated function
538:    G   - the newly evaluated gradient
539: */
540: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
541: {
542:   User              user_ptr = (User)ctx;
543:   TS                ts = user_ptr->ts;
544:   PetscScalar       *x_ptr,*g;
545:   const PetscScalar *y_ptr;

548:   VecGetArrayRead(P,&y_ptr);
549:   user_ptr->mu = y_ptr[0];
550:   VecRestoreArrayRead(P,&y_ptr);

552:   TSSetTime(ts,0.0);
553:   TSSetStepNumber(ts,0);
554:   TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
555:   TSSetFromOptions(ts);
556:   VecGetArray(user_ptr->U,&x_ptr);
557:   x_ptr[0] = 2.0;
558:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
559:   VecRestoreArray(user_ptr->U,&x_ptr);

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

563:   VecGetArrayRead(user_ptr->U,&y_ptr);
564:   *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]);

566:   /*   Reset initial conditions for the adjoint integration */
567:   VecGetArray(user_ptr->Lambda[0],&x_ptr);
568:   x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
569:   x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
570:   VecRestoreArrayRead(user_ptr->U,&y_ptr);
571:   VecRestoreArray(user_ptr->Lambda[0],&x_ptr);

573:   VecGetArray(user_ptr->Mup[0],&x_ptr);
574:   x_ptr[0] = 0.0;
575:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
576:   TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);

578:   TSAdjointSolve(ts);

580:   VecGetArray(user_ptr->Mup[0],&x_ptr);
581:   VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
582:   VecGetArray(G,&g);
583:   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];
584:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
585:   VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
586:   VecRestoreArray(G,&g);
587:   return 0;
588: }

590: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
591: {
592:   User           user_ptr = (User)ctx;
593:   PetscScalar    harr[1];
594:   const PetscInt rows[1] = {0};
595:   PetscInt       col = 0;

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

601:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
602:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
603:   if (H != Hpre) {
604:     MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
605:     MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
606:   }
607:   return 0;
608: }

610: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
611: {
612:   TS                ts = ctx->ts;
613:   const PetscScalar *z_ptr;
614:   PetscScalar       *x_ptr,*y_ptr,dzdp,dzdp2;
615:   Mat               tlmsen;

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

621:   /* The directional vector should be 1 since it is one-dimensional */
622:   VecGetArray(ctx->Dir,&x_ptr);
623:   x_ptr[0] = 1.;
624:   VecRestoreArray(ctx->Dir,&x_ptr);

626:   VecGetArrayRead(P,&z_ptr);
627:   ctx->mu = z_ptr[0];
628:   VecRestoreArrayRead(P,&z_ptr);

630:   dzdp  = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
631:   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);

633:   TSSetTime(ts,0.0);
634:   TSSetStepNumber(ts,0);
635:   TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
636:   TSSetFromOptions(ts);
637:   TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);

639:   MatZeroEntries(ctx->Jacp);
640:   MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
641:   MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
642:   MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);

644:   TSAdjointSetForward(ts,ctx->Jacp);
645:   VecGetArray(ctx->U,&y_ptr);
646:   y_ptr[0] = 2.0;
647:   y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
648:   VecRestoreArray(ctx->U,&y_ptr);
649:   TSSolve(ts,ctx->U);

651:   /* Set terminal conditions for first- and second-order adjonts */
652:   VecGetArrayRead(ctx->U,&z_ptr);
653:   VecGetArray(ctx->Lambda[0],&y_ptr);
654:   y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
655:   y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
656:   VecRestoreArray(ctx->Lambda[0],&y_ptr);
657:   VecRestoreArrayRead(ctx->U,&z_ptr);
658:   VecGetArray(ctx->Mup[0],&y_ptr);
659:   y_ptr[0] = 0.0;
660:   VecRestoreArray(ctx->Mup[0],&y_ptr);
661:   TSForwardGetSensitivities(ts,NULL,&tlmsen);
662:   MatDenseGetColumn(tlmsen,0,&x_ptr);
663:   VecGetArray(ctx->Lambda2[0],&y_ptr);
664:   y_ptr[0] = 2.*x_ptr[0];
665:   y_ptr[1] = 2.*x_ptr[1];
666:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
667:   VecGetArray(ctx->Mup2[0],&y_ptr);
668:   y_ptr[0] = 0.0;
669:   VecRestoreArray(ctx->Mup2[0],&y_ptr);
670:   MatDenseRestoreColumn(tlmsen,&x_ptr);
671:   TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
672:   if (ctx->implicitform) {
673:     TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
674:   } else {
675:     TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
676:   }
677:   TSAdjointSolve(ts);

679:   VecGetArray(ctx->Lambda[0],&x_ptr);
680:   VecGetArray(ctx->Lambda2[0],&y_ptr);
681:   VecGetArrayRead(ctx->Mup2[0],&z_ptr);

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

685:   VecRestoreArray(ctx->Lambda2[0],&x_ptr);
686:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
687:   VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);

689:   /* Disable second-order adjoint mode */
690:   TSAdjointReset(ts);
691:   TSAdjointResetForward(ts);
692:   return 0;
693: }

695: /*TEST
696:     build:
697:       requires: !complex !single
698:     test:
699:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
700:       output_file: output/ex20opt_p_1.out

702:     test:
703:       suffix: 2
704:       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
705:       output_file: output/ex20opt_p_2.out

707:     test:
708:       suffix: 3
709:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
710:       output_file: output/ex20opt_p_3.out

712:     test:
713:       suffix: 4
714:       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
715:       output_file: output/ex20opt_p_4.out

717: TEST*/