Actual source code: ex1.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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;
 30:   PetscErrorCode    ierr;

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

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

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

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

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

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

 75:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 77:   return(0);
 78: }

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

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

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

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

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

117:   dJpdP[0][0][1] = -PetscSinReal(w[step-1]);
118:   dJpdP[0][1][0] = -PetscSinReal(w[step-1]);
119:   dJpdP[0][1][1] = -v[step-1]*PetscCosReal(w[step-1]);
120:   dJpdP[1][0][1] = PetscCosReal(w[step-1]);
121:   dJpdP[1][1][0] = PetscCosReal(w[step-1]);
122:   dJpdP[1][1][1] = -v[step-1]*PetscSinReal(w[step-1]);

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

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

148:   TSGetStepNumber(ts,&step);
149:   VecGetArrayRead(U,&u);
150:   VecGetArrayRead(actx->V,&v);
151:   VecGetArrayRead(actx->W,&w);
152:   VecGetArrayRead(Vr,&vr);
153:   VecGetArray(VHV[0],&vhv);

155:   dRudU[0][0] = 2.0;
156:   dRudU[1][1] = 2.0;

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

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

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

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


188: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,void *ctx)
189: {
190:   Aircraft          actx = (Aircraft)ctx;
191:   PetscScalar       *r;
192:   PetscReal         dx,dy;
193:   const PetscScalar *u;
194:   PetscErrorCode    ierr;

197:   VecGetArrayRead(U,&u);
198:   VecGetArray(R,&r);
199:   dx   = u[0] - actx->lv*t*PetscCosReal(actx->lw);
200:   dy   = u[1] - actx->lv*t*PetscSinReal(actx->lw);
201:   r[0] = dx*dx+dy*dy;
202:   VecRestoreArray(R,&r);
203:   VecRestoreArrayRead(U,&u);
204:   return(0);
205: }

207: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,void *ctx)
208: {
209:   Aircraft          actx = (Aircraft)ctx;
210:   PetscScalar       drdu[2][1];
211:   const PetscScalar *u;
212:   PetscReal         dx,dy;
213:   PetscInt          row[] = {0,1},col[] = {0};
214:   PetscErrorCode    ierr;

217:   VecGetArrayRead(U,&u);
218:   dx      = u[0] - actx->lv*t*PetscCosReal(actx->lw);
219:   dy      = u[1] - actx->lv*t*PetscSinReal(actx->lw);
220:   drdu[0][0] = 2.*dx;
221:   drdu[1][0] = 2.*dy;
222:   MatSetValues(DRDU,2,row,1,col,&drdu[0][0],INSERT_VALUES);
223:   VecRestoreArrayRead(U,&u);
224:   MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
225:   MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
226:   return(0);
227: }

229: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,void *ctx)
230: {

234:   MatZeroEntries(DRDP);
235:   MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
236:   MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
237:   return(0);
238: }

240: int main(int argc,char **argv)
241: {
242:   Vec                P,PL,PU;
243:   struct _n_aircraft aircraft;
244:   PetscMPIInt        size;
245:   Tao                tao;
246:   KSP                ksp;
247:   PC                 pc;
248:   PetscScalar        *u,*p;
249:   PetscInt           i;
250:   PetscErrorCode     ierr;

252:   /* Initialize program */
253:   PetscInitialize(&argc,&argv,NULL,NULL);if(ierr) return ierr;
254:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
255:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

257:   /* Parameter settings */
258:   aircraft.ftime = 1.;   /* time interval in hour */
259:   aircraft.nsteps = 10; /* number of steps */
260:   aircraft.lv = 2.0; /* leader speed in kmph */
261:   aircraft.lw = PETSC_PI/4.; /* leader heading angle */

263:   PetscOptionsGetReal(NULL,NULL,"-ftime",&aircraft.ftime,NULL);
264:   PetscOptionsGetInt(NULL,NULL,"-nsteps",&aircraft.nsteps,NULL);
265:   PetscOptionsHasName(NULL,NULL,"-matrixfree",&aircraft.mf);
266:   PetscOptionsHasName(NULL,NULL,"-exacthessian",&aircraft.eh);

268:   /* Create TAO solver and set desired solution method */
269:   TaoCreate(PETSC_COMM_WORLD,&tao);
270:   TaoSetType(tao,TAOBQNLS);

272:   /* Create necessary matrix and vectors, solve same ODE on every process */
273:   MatCreate(PETSC_COMM_WORLD,&aircraft.A);
274:   MatSetSizes(aircraft.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
275:   MatSetFromOptions(aircraft.A);
276:   MatSetUp(aircraft.A);
277:   MatAssemblyBegin(aircraft.A,MAT_FINAL_ASSEMBLY);
278:   MatAssemblyEnd(aircraft.A,MAT_FINAL_ASSEMBLY);
279:   MatShift(aircraft.A,1);
280:   MatShift(aircraft.A,-1);

282:   MatCreate(PETSC_COMM_WORLD,&aircraft.Jacp);
283:   MatSetSizes(aircraft.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2*aircraft.nsteps);
284:   MatSetFromOptions(aircraft.Jacp);
285:   MatSetUp(aircraft.Jacp);
286:   MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,1,NULL,&aircraft.DRDP);
287:   MatSetUp(aircraft.DRDP);
288:   MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&aircraft.DRDU);
289:   MatSetUp(aircraft.DRDU);

291:   /* Create timestepping solver context */
292:   TSCreate(PETSC_COMM_WORLD,&aircraft.ts);
293:   TSSetType(aircraft.ts,TSRK);
294:   TSSetRHSFunction(aircraft.ts,NULL,RHSFunction,&aircraft);
295:   TSSetRHSJacobian(aircraft.ts,aircraft.A,aircraft.A,TSComputeRHSJacobianConstant,&aircraft);
296:   TSSetRHSJacobianP(aircraft.ts,aircraft.Jacp,RHSJacobianP,&aircraft);
297:   TSSetExactFinalTime(aircraft.ts,TS_EXACTFINALTIME_MATCHSTEP);
298:   TSSetEquationType(aircraft.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */

300:   /* Set initial conditions */
301:   MatCreateVecs(aircraft.A,&aircraft.U,NULL);
302:   TSSetSolution(aircraft.ts,aircraft.U);
303:   VecGetArray(aircraft.U,&u);
304:   u[0] = 1.5;
305:   u[1] = 0;
306:   VecRestoreArray(aircraft.U,&u);
307:   VecCreate(PETSC_COMM_WORLD,&aircraft.V);
308:   VecSetSizes(aircraft.V,PETSC_DECIDE,aircraft.nsteps);
309:   VecSetUp(aircraft.V);
310:   VecDuplicate(aircraft.V,&aircraft.W);
311:   VecSet(aircraft.V,1.);
312:   VecSet(aircraft.W,PETSC_PI/4.);

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

317:   /* Set sensitivity context */
318:   TSCreateQuadratureTS(aircraft.ts,PETSC_FALSE,&aircraft.quadts);
319:   TSSetRHSFunction(aircraft.quadts,NULL,(TSRHSFunction)CostIntegrand,&aircraft);
320:   TSSetRHSJacobian(aircraft.quadts,aircraft.DRDU,aircraft.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&aircraft);
321:   TSSetRHSJacobianP(aircraft.quadts,aircraft.DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&aircraft);
322:   MatCreateVecs(aircraft.A,&aircraft.Lambda[0],NULL);
323:   MatCreateVecs(aircraft.Jacp,&aircraft.Mup[0],NULL);
324:   if (aircraft.eh) {
325:     MatCreateVecs(aircraft.A,&aircraft.rhshp1[0],NULL);
326:     MatCreateVecs(aircraft.A,&aircraft.rhshp2[0],NULL);
327:     MatCreateVecs(aircraft.Jacp,&aircraft.rhshp3[0],NULL);
328:     MatCreateVecs(aircraft.Jacp,&aircraft.rhshp4[0],NULL);
329:     MatCreateVecs(aircraft.DRDU,&aircraft.inthp1[0],NULL);
330:     MatCreateVecs(aircraft.DRDU,&aircraft.inthp2[0],NULL);
331:     MatCreateVecs(aircraft.DRDP,&aircraft.inthp3[0],NULL);
332:     MatCreateVecs(aircraft.DRDP,&aircraft.inthp4[0],NULL);
333:     MatCreateVecs(aircraft.Jacp,&aircraft.Dir,NULL);
334:     TSSetRHSHessianProduct(aircraft.ts,aircraft.rhshp1,RHSHessianProductUU,aircraft.rhshp2,RHSHessianProductUP,aircraft.rhshp3,RHSHessianProductPU,aircraft.rhshp4,RHSHessianProductPP,&aircraft);
335:     TSSetRHSHessianProduct(aircraft.quadts,aircraft.inthp1,IntegrandHessianProductUU,aircraft.inthp2,IntegrandHessianProductUP,aircraft.inthp3,IntegrandHessianProductPU,aircraft.inthp4,IntegrandHessianProductPP,&aircraft);
336:     MatCreateVecs(aircraft.A,&aircraft.Lambda2[0],NULL);
337:     MatCreateVecs(aircraft.Jacp,&aircraft.Mup2[0],NULL);
338:   }
339:   TSSetFromOptions(aircraft.ts);
340:   TSSetMaxTime(aircraft.ts,aircraft.ftime);
341:   TSSetTimeStep(aircraft.ts,aircraft.ftime/aircraft.nsteps);

343:   /* Set initial solution guess */
344:   MatCreateVecs(aircraft.Jacp,&P,NULL);
345:   VecGetArray(P,&p);
346:   for (i=0; i<aircraft.nsteps; i++) {
347:     p[2*i] = 2.0;
348:     p[2*i+1] = PETSC_PI/2.0;
349:   }
350:   VecRestoreArray(P,&p);
351:   VecDuplicate(P,&PU);
352:   VecDuplicate(P,&PL);
353:   VecGetArray(PU,&p);
354:   for (i=0; i<aircraft.nsteps; i++) {
355:     p[2*i] = 2.0;
356:     p[2*i+1] = PETSC_PI;
357:   }
358:   VecRestoreArray(PU,&p);
359:   VecGetArray(PL,&p);
360:   for (i=0; i<aircraft.nsteps; i++) {
361:     p[2*i] = 0.0;
362:     p[2*i+1] = -PETSC_PI;
363:   }
364:   VecRestoreArray(PL,&p);

366:   TaoSetInitialVector(tao,P);
367:   TaoSetVariableBounds(tao,PL,PU);
368:   /* Set routine for function and gradient evaluation */
369:   TaoSetObjectiveAndGradientRoutine(tao,FormObjFunctionGradient,(void *)&aircraft);

371:   if (aircraft.eh) {
372:     if (aircraft.mf) {
373:       MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,2*aircraft.nsteps,(void*)&aircraft,&aircraft.H);
374:       MatShellSetOperation(aircraft.H,MATOP_MULT,(void(*)(void))MyMatMult);
375:       MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
376:       TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,MatrixFreeObjHessian,(void*)&aircraft);
377:     } else {
378:       MatCreateDense(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,2*aircraft.nsteps,2*aircraft.nsteps,NULL,&(aircraft.H));
379:       MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
380:       TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,FormObjHessian,(void *)&aircraft);
381:     }
382:   }

384:   /* Check for any TAO command line options */
385:   TaoGetKSP(tao,&ksp);
386:   if (ksp) {
387:     KSPGetPC(ksp,&pc);
388:     PCSetType(pc,PCNONE);
389:   }
390:   TaoSetFromOptions(tao);

392:   TaoSolve(tao);
393:   VecView(P,PETSC_VIEWER_STDOUT_WORLD);

395:   /* Free TAO data structures */
396:   TaoDestroy(&tao);

398:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
399:      Free work space.  All PETSc objects should be destroyed when they
400:      are no longer needed.
401:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
402:   TSDestroy(&aircraft.ts);
403:   MatDestroy(&aircraft.A);
404:   VecDestroy(&aircraft.U);
405:   VecDestroy(&aircraft.V);
406:   VecDestroy(&aircraft.W);
407:   VecDestroy(&P);
408:   VecDestroy(&PU);
409:   VecDestroy(&PL);
410:   MatDestroy(&aircraft.Jacp);
411:   MatDestroy(&aircraft.DRDU);
412:   MatDestroy(&aircraft.DRDP);
413:   VecDestroy(&aircraft.Lambda[0]);
414:   VecDestroy(&aircraft.Mup[0]);
415:   VecDestroy(&P);
416:   if (aircraft.eh) {
417:     VecDestroy(&aircraft.Lambda2[0]);
418:     VecDestroy(&aircraft.Mup2[0]);
419:     VecDestroy(&aircraft.Dir);
420:     VecDestroy(&aircraft.rhshp1[0]);
421:     VecDestroy(&aircraft.rhshp2[0]);
422:     VecDestroy(&aircraft.rhshp3[0]);
423:     VecDestroy(&aircraft.rhshp4[0]);
424:     VecDestroy(&aircraft.inthp1[0]);
425:     VecDestroy(&aircraft.inthp2[0]);
426:     VecDestroy(&aircraft.inthp3[0]);
427:     VecDestroy(&aircraft.inthp4[0]);
428:     MatDestroy(&aircraft.H);
429:   }
430:   PetscFinalize();
431:   return ierr;
432: }

434: /*
435:    FormObjFunctionGradient - Evaluates the function and corresponding gradient.

437:    Input Parameters:
438:    tao - the Tao context
439:    P   - the input vector
440:    ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradientRoutine()

442:    Output Parameters:
443:    f   - the newly evaluated function
444:    G   - the newly evaluated gradient
445: */
446: PetscErrorCode FormObjFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
447: {
448:   Aircraft          actx = (Aircraft)ctx;
449:   TS                ts = actx->ts;
450:   Vec               Q;
451:   const PetscScalar *p,*q;
452:   PetscScalar       *u,*v,*w;
453:   PetscInt          i;
454:   PetscErrorCode    ierr;

457:   VecGetArrayRead(P,&p);
458:   VecGetArray(actx->V,&v);
459:   VecGetArray(actx->W,&w);
460:   for (i=0; i<actx->nsteps; i++) {
461:     v[i] = p[2*i];
462:     w[i] = p[2*i+1];
463:   }
464:   VecRestoreArrayRead(P,&p);
465:   VecRestoreArray(actx->V,&v);
466:   VecRestoreArray(actx->W,&w);

468:   TSSetTime(ts,0.0);
469:   TSSetStepNumber(ts,0);
470:   TSSetFromOptions(ts);
471:   TSSetTimeStep(ts,actx->ftime/actx->nsteps);

473:   /* reinitialize system state */
474:   VecGetArray(actx->U,&u);
475:   u[0] = 2.0;
476:   u[1] = 0;
477:   VecRestoreArray(actx->U,&u);

479:   /* reinitialize the integral value */
480:   TSGetCostIntegral(ts,&Q);
481:   VecSet(Q,0.0);

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

485:   /* Reset initial conditions for the adjoint integration */
486:   VecSet(actx->Lambda[0],0.0);
487:   VecSet(actx->Mup[0],0.0);
488:   TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);

490:   TSAdjointSolve(ts);
491:   VecCopy(actx->Mup[0],G);
492:   TSGetCostIntegral(ts,&Q);
493:   VecGetArrayRead(Q,&q);
494:   *f   = q[0];
495:   VecRestoreArrayRead(Q,&q);
496:   return(0);
497: }

499: PetscErrorCode FormObjHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
500: {
501:   Aircraft          actx = (Aircraft)ctx;
502:   const PetscScalar *p;
503:   PetscScalar       *harr,*v,*w,one = 1.0;
504:   PetscInt          ind[1];
505:   PetscInt          *cols,i;
506:   Vec               Dir;
507:   PetscErrorCode    ierr;

510:   /* set up control parameters */
511:   VecGetArrayRead(P,&p);
512:   VecGetArray(actx->V,&v);
513:   VecGetArray(actx->W,&w);
514:   for (i=0; i<actx->nsteps; i++) {
515:     v[i] = p[2*i];
516:     w[i] = p[2*i+1];
517:   }
518:   VecRestoreArrayRead(P,&p);
519:   VecRestoreArray(actx->V,&v);
520:   VecRestoreArray(actx->W,&w);

522:   PetscMalloc1(2*actx->nsteps,&harr);
523:   PetscMalloc1(2*actx->nsteps,&cols);
524:   for (i=0; i<2*actx->nsteps; i++) cols[i] = i;
525:   VecDuplicate(P,&Dir);
526:   for (i=0; i<2*actx->nsteps; i++) {
527:     ind[0] = i;
528:     VecSet(Dir,0.0);
529:     VecSetValues(Dir,1,ind,&one,INSERT_VALUES);
530:     VecAssemblyBegin(Dir);
531:     VecAssemblyEnd(Dir);
532:     ComputeObjHessianWithSOA(Dir,harr,actx);
533:     MatSetValues(H,1,ind,2*actx->nsteps,cols,harr,INSERT_VALUES);
534:     MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
535:     MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
536:     if (H != Hpre) {
537:       MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
538:       MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
539:     }
540:   }
541:   PetscFree(cols);
542:   PetscFree(harr);
543:   VecDestroy(&Dir);
544:   return(0);
545: }

547: PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
548: {
549:   Aircraft          actx = (Aircraft)ctx;
550:   PetscScalar       *v,*w;
551:   const PetscScalar *p;
552:   PetscInt          i;
553:   PetscErrorCode    ierr;

556:   VecGetArrayRead(P,&p);
557:   VecGetArray(actx->V,&v);
558:   VecGetArray(actx->W,&w);
559:   for (i=0; i<actx->nsteps; i++) {
560:     v[i] = p[2*i];
561:     w[i] = p[2*i+1];
562:   }
563:   VecRestoreArrayRead(P,&p);
564:   VecRestoreArray(actx->V,&v);
565:   VecRestoreArray(actx->W,&w);
566:   return(0);
567: }

569: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
570: {
571:   PetscScalar    *y;
572:   void           *ptr;

576:   MatShellGetContext(H_shell,&ptr);
577:   VecGetArray(Y,&y);
578:   ComputeObjHessianWithSOA(X,y,(Aircraft)ptr);
579:   VecRestoreArray(Y,&y);
580:   return(0);
581: }

583: PetscErrorCode ComputeObjHessianWithSOA(Vec Dir,PetscScalar arr[],Aircraft actx)
584: {
585:   TS                ts = actx->ts;
586:   const PetscScalar *z_ptr;
587:   PetscScalar       *u;
588:   Vec               Q;
589:   PetscInt          i;
590:   PetscErrorCode    ierr;

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

596:   TSSetTime(ts,0.0);
597:   TSSetStepNumber(ts,0);
598:   TSSetFromOptions(ts);
599:   TSSetTimeStep(ts,actx->ftime/actx->nsteps);
600:   TSSetCostHessianProducts(actx->ts,1,actx->Lambda2,actx->Mup2,Dir);

602:   /* reinitialize system state */
603:   VecGetArray(actx->U,&u);
604:   u[0] = 2.0;
605:   u[1] = 0;
606:   VecRestoreArray(actx->U,&u);

608:   /* reinitialize the integral value */
609:   TSGetCostIntegral(ts,&Q);
610:   VecSet(Q,0.0);

612:   /* initialize tlm variable */
613:   MatZeroEntries(actx->Jacp);
614:   MatAssemblyBegin(actx->Jacp,MAT_FINAL_ASSEMBLY);
615:   MatAssemblyEnd(actx->Jacp,MAT_FINAL_ASSEMBLY);
616:   TSAdjointSetForward(ts,actx->Jacp);

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

620:   /* Set terminal conditions for first- and second-order adjonts */
621:   VecSet(actx->Lambda[0],0.0);
622:   VecSet(actx->Mup[0],0.0);
623:   VecSet(actx->Lambda2[0],0.0);
624:   VecSet(actx->Mup2[0],0.0);
625:   TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);

627:   TSGetCostIntegral(ts,&Q);

629:   /* Reset initial conditions for the adjoint integration */
630:   TSAdjointSolve(ts);

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

637:   /* Disable second-order adjoint mode */
638:   TSAdjointReset(ts);
639:   TSAdjointResetForward(ts);
640:   return(0);
641: }

643: /*TEST
644:     build:
645:       requires: !complex !single

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

650:     test:
651:       suffix: 2
652:       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

654:     test:
655:       suffix: 3
656:       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
657: TEST*/