Actual source code: ex1.c

  1: #include <petsctao.h>
  2: #include <petscts.h>

  4: typedef struct _n_aircraft *Aircraft;
  5: struct _n_aircraft {
  6:   TS        ts, quadts;
  7:   Vec       V, W;   /* control variables V and W */
  8:   PetscInt  nsteps; /* number of time steps */
  9:   PetscReal ftime;
 10:   Mat       A, H;
 11:   Mat       Jacp, DRDU, DRDP;
 12:   Vec       U, Lambda[1], Mup[1], Lambda2[1], Mup2[1], Dir;
 13:   Vec       rhshp1[1], rhshp2[1], rhshp3[1], rhshp4[1], inthp1[1], inthp2[1], inthp3[1], inthp4[1];
 14:   PetscReal lv, lw;
 15:   PetscBool mf, eh;
 16: };

 18: PetscErrorCode FormObjFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 19: PetscErrorCode FormObjHessian(Tao, Vec, Mat, Mat, void *);
 20: PetscErrorCode ComputeObjHessianWithSOA(Vec, PetscScalar[], Aircraft);
 21: PetscErrorCode MatrixFreeObjHessian(Tao, Vec, Mat, Mat, void *);
 22: PetscErrorCode MyMatMult(Mat, Vec, Vec);

 24: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
 25: {
 26:   Aircraft           actx = (Aircraft)ctx;
 27:   const PetscScalar *u, *v, *w;
 28:   PetscScalar       *f;
 29:   PetscInt           step;

 32:   TSGetStepNumber(ts, &step);
 33:   VecGetArrayRead(U, &u);
 34:   VecGetArrayRead(actx->V, &v);
 35:   VecGetArrayRead(actx->W, &w);
 36:   VecGetArray(F, &f);
 37:   f[0] = v[step] * PetscCosReal(w[step]);
 38:   f[1] = v[step] * PetscSinReal(w[step]);
 39:   VecRestoreArrayRead(U, &u);
 40:   VecRestoreArrayRead(actx->V, &v);
 41:   VecRestoreArrayRead(actx->W, &w);
 42:   VecRestoreArray(F, &f);
 43:   return 0;
 44: }

 46: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec U, Mat A, void *ctx)
 47: {
 48:   Aircraft           actx = (Aircraft)ctx;
 49:   const PetscScalar *u, *v, *w;
 50:   PetscInt           step, rows[2] = {0, 1}, rowcol[2];
 51:   PetscScalar        Jp[2][2];

 54:   MatZeroEntries(A);
 55:   TSGetStepNumber(ts, &step);
 56:   VecGetArrayRead(U, &u);
 57:   VecGetArrayRead(actx->V, &v);
 58:   VecGetArrayRead(actx->W, &w);

 60:   Jp[0][0] = PetscCosReal(w[step]);
 61:   Jp[0][1] = -v[step] * PetscSinReal(w[step]);
 62:   Jp[1][0] = PetscSinReal(w[step]);
 63:   Jp[1][1] = v[step] * PetscCosReal(w[step]);

 65:   VecRestoreArrayRead(U, &u);
 66:   VecRestoreArrayRead(actx->V, &v);
 67:   VecRestoreArrayRead(actx->W, &w);

 69:   rowcol[0] = 2 * step;
 70:   rowcol[1] = 2 * step + 1;
 71:   MatSetValues(A, 2, rows, 2, rowcol, &Jp[0][0], INSERT_VALUES);

 73:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 74:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 75:   return 0;
 76: }

 78: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 79: {
 81:   return 0;
 82: }

 84: static PetscErrorCode RHSHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 85: {
 87:   return 0;
 88: }

 90: static PetscErrorCode RHSHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 91: {
 93:   return 0;
 94: }

 96: static PetscErrorCode RHSHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 97: {
 98:   Aircraft           actx = (Aircraft)ctx;
 99:   const PetscScalar *v, *w, *vl, *vr, *u;
100:   PetscScalar       *vhv;
101:   PetscScalar        dJpdP[2][2][2] = {{{0}}};
102:   PetscInt           step, i, j, k;

105:   TSGetStepNumber(ts, &step);
106:   VecGetArrayRead(U, &u);
107:   VecGetArrayRead(actx->V, &v);
108:   VecGetArrayRead(actx->W, &w);
109:   VecGetArrayRead(Vl[0], &vl);
110:   VecGetArrayRead(Vr, &vr);
111:   VecSet(VHV[0], 0.0);
112:   VecGetArray(VHV[0], &vhv);

114:   dJpdP[0][0][1] = -PetscSinReal(w[step]);
115:   dJpdP[0][1][0] = -PetscSinReal(w[step]);
116:   dJpdP[0][1][1] = -v[step] * PetscCosReal(w[step]);
117:   dJpdP[1][0][1] = PetscCosReal(w[step]);
118:   dJpdP[1][1][0] = PetscCosReal(w[step]);
119:   dJpdP[1][1][1] = -v[step] * PetscSinReal(w[step]);

121:   for (j = 0; j < 2; j++) {
122:     vhv[2 * step + j] = 0;
123:     for (k = 0; k < 2; k++)
124:       for (i = 0; i < 2; i++) vhv[2 * step + j] += vl[i] * dJpdP[i][j][k] * vr[2 * step + k];
125:   }
126:   VecRestoreArrayRead(U, &u);
127:   VecRestoreArrayRead(Vl[0], &vl);
128:   VecRestoreArrayRead(Vr, &vr);
129:   VecRestoreArray(VHV[0], &vhv);
130:   return 0;
131: }

133: /* Vl in NULL,updates to VHV must be added */
134: static PetscErrorCode IntegrandHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
135: {
136:   Aircraft           actx = (Aircraft)ctx;
137:   const PetscScalar *v, *w, *vr, *u;
138:   PetscScalar       *vhv;
139:   PetscScalar        dRudU[2][2] = {{0}};
140:   PetscInt           step, j, k;

143:   TSGetStepNumber(ts, &step);
144:   VecGetArrayRead(U, &u);
145:   VecGetArrayRead(actx->V, &v);
146:   VecGetArrayRead(actx->W, &w);
147:   VecGetArrayRead(Vr, &vr);
148:   VecGetArray(VHV[0], &vhv);

150:   dRudU[0][0] = 2.0;
151:   dRudU[1][1] = 2.0;

153:   for (j = 0; j < 2; j++) {
154:     vhv[j] = 0;
155:     for (k = 0; k < 2; k++) vhv[j] += dRudU[j][k] * vr[k];
156:   }
157:   VecRestoreArrayRead(U, &u);
158:   VecRestoreArrayRead(Vr, &vr);
159:   VecRestoreArray(VHV[0], &vhv);
160:   return 0;
161: }

163: static PetscErrorCode IntegrandHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
164: {
166:   return 0;
167: }

169: static PetscErrorCode IntegrandHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
170: {
172:   return 0;
173: }

175: static PetscErrorCode IntegrandHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
176: {
178:   return 0;
179: }

181: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
182: {
183:   Aircraft           actx = (Aircraft)ctx;
184:   PetscScalar       *r;
185:   PetscReal          dx, dy;
186:   const PetscScalar *u;

188:   VecGetArrayRead(U, &u);
189:   VecGetArray(R, &r);
190:   dx   = u[0] - actx->lv * t * PetscCosReal(actx->lw);
191:   dy   = u[1] - actx->lv * t * PetscSinReal(actx->lw);
192:   r[0] = dx * dx + dy * dy;
193:   VecRestoreArray(R, &r);
194:   VecRestoreArrayRead(U, &u);
195:   return 0;
196: }

198: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, void *ctx)
199: {
200:   Aircraft           actx = (Aircraft)ctx;
201:   PetscScalar        drdu[2][1];
202:   const PetscScalar *u;
203:   PetscReal          dx, dy;
204:   PetscInt           row[] = {0, 1}, col[] = {0};

206:   VecGetArrayRead(U, &u);
207:   dx         = u[0] - actx->lv * t * PetscCosReal(actx->lw);
208:   dy         = u[1] - actx->lv * t * PetscSinReal(actx->lw);
209:   drdu[0][0] = 2. * dx;
210:   drdu[1][0] = 2. * dy;
211:   MatSetValues(DRDU, 2, row, 1, col, &drdu[0][0], INSERT_VALUES);
212:   VecRestoreArrayRead(U, &u);
213:   MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY);
214:   MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY);
215:   return 0;
216: }

218: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, void *ctx)
219: {
220:   MatZeroEntries(DRDP);
221:   MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY);
222:   MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY);
223:   return 0;
224: }

226: int main(int argc, char **argv)
227: {
228:   Vec                P, PL, PU;
229:   struct _n_aircraft aircraft;
230:   PetscMPIInt        size;
231:   Tao                tao;
232:   KSP                ksp;
233:   PC                 pc;
234:   PetscScalar       *u, *p;
235:   PetscInt           i;

237:   /* Initialize program */
239:   PetscInitialize(&argc, &argv, NULL, NULL);
240:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

243:   /* Parameter settings */
244:   aircraft.ftime  = 1.;            /* time interval in hour */
245:   aircraft.nsteps = 10;            /* number of steps */
246:   aircraft.lv     = 2.0;           /* leader speed in kmph */
247:   aircraft.lw     = PETSC_PI / 4.; /* leader heading angle */

249:   PetscOptionsGetReal(NULL, NULL, "-ftime", &aircraft.ftime, NULL);
250:   PetscOptionsGetInt(NULL, NULL, "-nsteps", &aircraft.nsteps, NULL);
251:   PetscOptionsHasName(NULL, NULL, "-matrixfree", &aircraft.mf);
252:   PetscOptionsHasName(NULL, NULL, "-exacthessian", &aircraft.eh);

254:   /* Create TAO solver and set desired solution method */
255:   TaoCreate(PETSC_COMM_WORLD, &tao);
256:   TaoSetType(tao, TAOBQNLS);

258:   /* Create necessary matrix and vectors, solve same ODE on every process */
259:   MatCreate(PETSC_COMM_WORLD, &aircraft.A);
260:   MatSetSizes(aircraft.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
261:   MatSetFromOptions(aircraft.A);
262:   MatSetUp(aircraft.A);
263:   MatAssemblyBegin(aircraft.A, MAT_FINAL_ASSEMBLY);
264:   MatAssemblyEnd(aircraft.A, MAT_FINAL_ASSEMBLY);
265:   MatShift(aircraft.A, 1);
266:   MatShift(aircraft.A, -1);

268:   MatCreate(PETSC_COMM_WORLD, &aircraft.Jacp);
269:   MatSetSizes(aircraft.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 2 * aircraft.nsteps);
270:   MatSetFromOptions(aircraft.Jacp);
271:   MatSetUp(aircraft.Jacp);
272:   MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2 * aircraft.nsteps, 1, NULL, &aircraft.DRDP);
273:   MatSetUp(aircraft.DRDP);
274:   MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &aircraft.DRDU);
275:   MatSetUp(aircraft.DRDU);

277:   /* Create timestepping solver context */
278:   TSCreate(PETSC_COMM_WORLD, &aircraft.ts);
279:   TSSetType(aircraft.ts, TSRK);
280:   TSSetRHSFunction(aircraft.ts, NULL, RHSFunction, &aircraft);
281:   TSSetRHSJacobian(aircraft.ts, aircraft.A, aircraft.A, TSComputeRHSJacobianConstant, &aircraft);
282:   TSSetRHSJacobianP(aircraft.ts, aircraft.Jacp, RHSJacobianP, &aircraft);
283:   TSSetExactFinalTime(aircraft.ts, TS_EXACTFINALTIME_MATCHSTEP);
284:   TSSetEquationType(aircraft.ts, TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */

286:   /* Set initial conditions */
287:   MatCreateVecs(aircraft.A, &aircraft.U, NULL);
288:   TSSetSolution(aircraft.ts, aircraft.U);
289:   VecGetArray(aircraft.U, &u);
290:   u[0] = 1.5;
291:   u[1] = 0;
292:   VecRestoreArray(aircraft.U, &u);
293:   VecCreate(PETSC_COMM_WORLD, &aircraft.V);
294:   VecSetSizes(aircraft.V, PETSC_DECIDE, aircraft.nsteps);
295:   VecSetUp(aircraft.V);
296:   VecDuplicate(aircraft.V, &aircraft.W);
297:   VecSet(aircraft.V, 1.);
298:   VecSet(aircraft.W, PETSC_PI / 4.);

300:   /* Save trajectory of solution so that TSAdjointSolve() may be used */
301:   TSSetSaveTrajectory(aircraft.ts);

303:   /* Set sensitivity context */
304:   TSCreateQuadratureTS(aircraft.ts, PETSC_FALSE, &aircraft.quadts);
305:   TSSetRHSFunction(aircraft.quadts, NULL, (TSRHSFunction)CostIntegrand, &aircraft);
306:   TSSetRHSJacobian(aircraft.quadts, aircraft.DRDU, aircraft.DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &aircraft);
307:   TSSetRHSJacobianP(aircraft.quadts, aircraft.DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, &aircraft);
308:   MatCreateVecs(aircraft.A, &aircraft.Lambda[0], NULL);
309:   MatCreateVecs(aircraft.Jacp, &aircraft.Mup[0], NULL);
310:   if (aircraft.eh) {
311:     MatCreateVecs(aircraft.A, &aircraft.rhshp1[0], NULL);
312:     MatCreateVecs(aircraft.A, &aircraft.rhshp2[0], NULL);
313:     MatCreateVecs(aircraft.Jacp, &aircraft.rhshp3[0], NULL);
314:     MatCreateVecs(aircraft.Jacp, &aircraft.rhshp4[0], NULL);
315:     MatCreateVecs(aircraft.DRDU, &aircraft.inthp1[0], NULL);
316:     MatCreateVecs(aircraft.DRDU, &aircraft.inthp2[0], NULL);
317:     MatCreateVecs(aircraft.DRDP, &aircraft.inthp3[0], NULL);
318:     MatCreateVecs(aircraft.DRDP, &aircraft.inthp4[0], NULL);
319:     MatCreateVecs(aircraft.Jacp, &aircraft.Dir, NULL);
320:     TSSetRHSHessianProduct(aircraft.ts, aircraft.rhshp1, RHSHessianProductUU, aircraft.rhshp2, RHSHessianProductUP, aircraft.rhshp3, RHSHessianProductPU, aircraft.rhshp4, RHSHessianProductPP, &aircraft);
321:     TSSetRHSHessianProduct(aircraft.quadts, aircraft.inthp1, IntegrandHessianProductUU, aircraft.inthp2, IntegrandHessianProductUP, aircraft.inthp3, IntegrandHessianProductPU, aircraft.inthp4, IntegrandHessianProductPP, &aircraft);
322:     MatCreateVecs(aircraft.A, &aircraft.Lambda2[0], NULL);
323:     MatCreateVecs(aircraft.Jacp, &aircraft.Mup2[0], NULL);
324:   }
325:   TSSetFromOptions(aircraft.ts);
326:   TSSetMaxTime(aircraft.ts, aircraft.ftime);
327:   TSSetTimeStep(aircraft.ts, aircraft.ftime / aircraft.nsteps);

329:   /* Set initial solution guess */
330:   MatCreateVecs(aircraft.Jacp, &P, NULL);
331:   VecGetArray(P, &p);
332:   for (i = 0; i < aircraft.nsteps; i++) {
333:     p[2 * i]     = 2.0;
334:     p[2 * i + 1] = PETSC_PI / 2.0;
335:   }
336:   VecRestoreArray(P, &p);
337:   VecDuplicate(P, &PU);
338:   VecDuplicate(P, &PL);
339:   VecGetArray(PU, &p);
340:   for (i = 0; i < aircraft.nsteps; i++) {
341:     p[2 * i]     = 2.0;
342:     p[2 * i + 1] = PETSC_PI;
343:   }
344:   VecRestoreArray(PU, &p);
345:   VecGetArray(PL, &p);
346:   for (i = 0; i < aircraft.nsteps; i++) {
347:     p[2 * i]     = 0.0;
348:     p[2 * i + 1] = -PETSC_PI;
349:   }
350:   VecRestoreArray(PL, &p);

352:   TaoSetSolution(tao, P);
353:   TaoSetVariableBounds(tao, PL, PU);
354:   /* Set routine for function and gradient evaluation */
355:   TaoSetObjectiveAndGradient(tao, NULL, FormObjFunctionGradient, (void *)&aircraft);

357:   if (aircraft.eh) {
358:     if (aircraft.mf) {
359:       MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2 * aircraft.nsteps, 2 * aircraft.nsteps, (void *)&aircraft, &aircraft.H);
360:       MatShellSetOperation(aircraft.H, MATOP_MULT, (void (*)(void))MyMatMult);
361:       MatSetOption(aircraft.H, MAT_SYMMETRIC, PETSC_TRUE);
362:       TaoSetHessian(tao, aircraft.H, aircraft.H, MatrixFreeObjHessian, (void *)&aircraft);
363:     } else {
364:       MatCreateDense(MPI_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, 2 * aircraft.nsteps, 2 * aircraft.nsteps, NULL, &(aircraft.H));
365:       MatSetOption(aircraft.H, MAT_SYMMETRIC, PETSC_TRUE);
366:       TaoSetHessian(tao, aircraft.H, aircraft.H, FormObjHessian, (void *)&aircraft);
367:     }
368:   }

370:   /* Check for any TAO command line options */
371:   TaoGetKSP(tao, &ksp);
372:   if (ksp) {
373:     KSPGetPC(ksp, &pc);
374:     PCSetType(pc, PCNONE);
375:   }
376:   TaoSetFromOptions(tao);

378:   TaoSolve(tao);
379:   VecView(P, PETSC_VIEWER_STDOUT_WORLD);

381:   /* Free TAO data structures */
382:   TaoDestroy(&tao);

384:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
385:      Free work space.  All PETSc objects should be destroyed when they
386:      are no longer needed.
387:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
388:   TSDestroy(&aircraft.ts);
389:   MatDestroy(&aircraft.A);
390:   VecDestroy(&aircraft.U);
391:   VecDestroy(&aircraft.V);
392:   VecDestroy(&aircraft.W);
393:   VecDestroy(&P);
394:   VecDestroy(&PU);
395:   VecDestroy(&PL);
396:   MatDestroy(&aircraft.Jacp);
397:   MatDestroy(&aircraft.DRDU);
398:   MatDestroy(&aircraft.DRDP);
399:   VecDestroy(&aircraft.Lambda[0]);
400:   VecDestroy(&aircraft.Mup[0]);
401:   VecDestroy(&P);
402:   if (aircraft.eh) {
403:     VecDestroy(&aircraft.Lambda2[0]);
404:     VecDestroy(&aircraft.Mup2[0]);
405:     VecDestroy(&aircraft.Dir);
406:     VecDestroy(&aircraft.rhshp1[0]);
407:     VecDestroy(&aircraft.rhshp2[0]);
408:     VecDestroy(&aircraft.rhshp3[0]);
409:     VecDestroy(&aircraft.rhshp4[0]);
410:     VecDestroy(&aircraft.inthp1[0]);
411:     VecDestroy(&aircraft.inthp2[0]);
412:     VecDestroy(&aircraft.inthp3[0]);
413:     VecDestroy(&aircraft.inthp4[0]);
414:     MatDestroy(&aircraft.H);
415:   }
416:   PetscFinalize();
417:   return 0;
418: }

420: /*
421:    FormObjFunctionGradient - Evaluates the function and corresponding gradient.

423:    Input Parameters:
424:    tao - the Tao context
425:    P   - the input vector
426:    ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradient()

428:    Output Parameters:
429:    f   - the newly evaluated function
430:    G   - the newly evaluated gradient
431: */
432: PetscErrorCode FormObjFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx)
433: {
434:   Aircraft           actx = (Aircraft)ctx;
435:   TS                 ts   = actx->ts;
436:   Vec                Q;
437:   const PetscScalar *p, *q;
438:   PetscScalar       *u, *v, *w;
439:   PetscInt           i;

442:   VecGetArrayRead(P, &p);
443:   VecGetArray(actx->V, &v);
444:   VecGetArray(actx->W, &w);
445:   for (i = 0; i < actx->nsteps; i++) {
446:     v[i] = p[2 * i];
447:     w[i] = p[2 * i + 1];
448:   }
449:   VecRestoreArrayRead(P, &p);
450:   VecRestoreArray(actx->V, &v);
451:   VecRestoreArray(actx->W, &w);

453:   TSSetTime(ts, 0.0);
454:   TSSetStepNumber(ts, 0);
455:   TSSetFromOptions(ts);
456:   TSSetTimeStep(ts, actx->ftime / actx->nsteps);

458:   /* reinitialize system state */
459:   VecGetArray(actx->U, &u);
460:   u[0] = 2.0;
461:   u[1] = 0;
462:   VecRestoreArray(actx->U, &u);

464:   /* reinitialize the integral value */
465:   TSGetCostIntegral(ts, &Q);
466:   VecSet(Q, 0.0);

468:   TSSolve(ts, actx->U);

470:   /* Reset initial conditions for the adjoint integration */
471:   VecSet(actx->Lambda[0], 0.0);
472:   VecSet(actx->Mup[0], 0.0);
473:   TSSetCostGradients(ts, 1, actx->Lambda, actx->Mup);

475:   TSAdjointSolve(ts);
476:   VecCopy(actx->Mup[0], G);
477:   TSGetCostIntegral(ts, &Q);
478:   VecGetArrayRead(Q, &q);
479:   *f = q[0];
480:   VecRestoreArrayRead(Q, &q);
481:   return 0;
482: }

484: PetscErrorCode FormObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
485: {
486:   Aircraft           actx = (Aircraft)ctx;
487:   const PetscScalar *p;
488:   PetscScalar       *harr, *v, *w, one = 1.0;
489:   PetscInt           ind[1];
490:   PetscInt          *cols, i;
491:   Vec                Dir;

494:   /* set up control parameters */
495:   VecGetArrayRead(P, &p);
496:   VecGetArray(actx->V, &v);
497:   VecGetArray(actx->W, &w);
498:   for (i = 0; i < actx->nsteps; i++) {
499:     v[i] = p[2 * i];
500:     w[i] = p[2 * i + 1];
501:   }
502:   VecRestoreArrayRead(P, &p);
503:   VecRestoreArray(actx->V, &v);
504:   VecRestoreArray(actx->W, &w);

506:   PetscMalloc1(2 * actx->nsteps, &harr);
507:   PetscMalloc1(2 * actx->nsteps, &cols);
508:   for (i = 0; i < 2 * actx->nsteps; i++) cols[i] = i;
509:   VecDuplicate(P, &Dir);
510:   for (i = 0; i < 2 * actx->nsteps; i++) {
511:     ind[0] = i;
512:     VecSet(Dir, 0.0);
513:     VecSetValues(Dir, 1, ind, &one, INSERT_VALUES);
514:     VecAssemblyBegin(Dir);
515:     VecAssemblyEnd(Dir);
516:     ComputeObjHessianWithSOA(Dir, harr, actx);
517:     MatSetValues(H, 1, ind, 2 * actx->nsteps, cols, harr, INSERT_VALUES);
518:     MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
519:     MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
520:     if (H != Hpre) {
521:       MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY);
522:       MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY);
523:     }
524:   }
525:   PetscFree(cols);
526:   PetscFree(harr);
527:   VecDestroy(&Dir);
528:   return 0;
529: }

531: PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
532: {
533:   Aircraft           actx = (Aircraft)ctx;
534:   PetscScalar       *v, *w;
535:   const PetscScalar *p;
536:   PetscInt           i;

538:   VecGetArrayRead(P, &p);
539:   VecGetArray(actx->V, &v);
540:   VecGetArray(actx->W, &w);
541:   for (i = 0; i < actx->nsteps; i++) {
542:     v[i] = p[2 * i];
543:     w[i] = p[2 * i + 1];
544:   }
545:   VecRestoreArrayRead(P, &p);
546:   VecRestoreArray(actx->V, &v);
547:   VecRestoreArray(actx->W, &w);
548:   return 0;
549: }

551: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
552: {
553:   PetscScalar *y;
554:   void        *ptr;

556:   MatShellGetContext(H_shell, &ptr);
557:   VecGetArray(Y, &y);
558:   ComputeObjHessianWithSOA(X, y, (Aircraft)ptr);
559:   VecRestoreArray(Y, &y);
560:   return 0;
561: }

563: PetscErrorCode ComputeObjHessianWithSOA(Vec Dir, PetscScalar arr[], Aircraft actx)
564: {
565:   TS                 ts = actx->ts;
566:   const PetscScalar *z_ptr;
567:   PetscScalar       *u;
568:   Vec                Q;
569:   PetscInt           i;

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

575:   TSSetTime(ts, 0.0);
576:   TSSetStepNumber(ts, 0);
577:   TSSetFromOptions(ts);
578:   TSSetTimeStep(ts, actx->ftime / actx->nsteps);
579:   TSSetCostHessianProducts(actx->ts, 1, actx->Lambda2, actx->Mup2, Dir);

581:   /* reinitialize system state */
582:   VecGetArray(actx->U, &u);
583:   u[0] = 2.0;
584:   u[1] = 0;
585:   VecRestoreArray(actx->U, &u);

587:   /* reinitialize the integral value */
588:   TSGetCostIntegral(ts, &Q);
589:   VecSet(Q, 0.0);

591:   /* initialize tlm variable */
592:   MatZeroEntries(actx->Jacp);
593:   MatAssemblyBegin(actx->Jacp, MAT_FINAL_ASSEMBLY);
594:   MatAssemblyEnd(actx->Jacp, MAT_FINAL_ASSEMBLY);
595:   TSAdjointSetForward(ts, actx->Jacp);

597:   TSSolve(ts, actx->U);

599:   /* Set terminal conditions for first- and second-order adjonts */
600:   VecSet(actx->Lambda[0], 0.0);
601:   VecSet(actx->Mup[0], 0.0);
602:   VecSet(actx->Lambda2[0], 0.0);
603:   VecSet(actx->Mup2[0], 0.0);
604:   TSSetCostGradients(ts, 1, actx->Lambda, actx->Mup);

606:   TSGetCostIntegral(ts, &Q);

608:   /* Reset initial conditions for the adjoint integration */
609:   TSAdjointSolve(ts);

611:   /* initial condition does not depend on p, so that lambda is not needed to assemble G */
612:   VecGetArrayRead(actx->Mup2[0], &z_ptr);
613:   for (i = 0; i < 2 * actx->nsteps; i++) arr[i] = z_ptr[i];
614:   VecRestoreArrayRead(actx->Mup2[0], &z_ptr);

616:   /* Disable second-order adjoint mode */
617:   TSAdjointReset(ts);
618:   TSAdjointResetForward(ts);
619:   return 0;
620: }

622: /*TEST
623:     build:
624:       requires: !complex !single

626:     test:
627:       args:  -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_gatol 1e-7

629:     test:
630:       suffix: 2
631:       args:  -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian

633:     test:
634:       suffix: 3
635:       args:  -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian -matrixfree
636: TEST*/