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++)
125:         vhv[2*step+j] += vl[i]*dJpdP[i][j][k]*vr[2*step+k];
126:   }
127:   VecRestoreArrayRead(U,&u);
128:   VecRestoreArrayRead(Vl[0],&vl);
129:   VecRestoreArrayRead(Vr,&vr);
130:   VecRestoreArray(VHV[0],&vhv);
131:   return 0;
132: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

607:   TSGetCostIntegral(ts,&Q);

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

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

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

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

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

630:     test:
631:       suffix: 2
632:       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

634:     test:
635:       suffix: 3
636:       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
637: TEST*/