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;
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]);
63: Jp[0][1] = -v[step]*PetscSinReal(w[step]);
64: Jp[1][0] = PetscSinReal(w[step]);
65: Jp[1][1] = v[step]*PetscCosReal(w[step]);
67: VecRestoreArrayRead(U,&u);
68: VecRestoreArrayRead(actx->V,&v);
69: VecRestoreArrayRead(actx->W,&w);
71: rowcol[0] = 2*step;
72: rowcol[1] = 2*step+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]);
118: dJpdP[0][1][0] = -PetscSinReal(w[step]);
119: dJpdP[0][1][1] = -v[step]*PetscCosReal(w[step]);
120: dJpdP[1][0][1] = PetscCosReal(w[step]);
121: dJpdP[1][1][0] = PetscCosReal(w[step]);
122: dJpdP[1][1][1] = -v[step]*PetscSinReal(w[step]);
124: for (j=0; j<2; j++) {
125: vhv[2*step+j] = 0;
126: for (k=0; k<2; k++)
127: for (i=0; i<2; i++)
128: vhv[2*step+j] += vl[i]*dJpdP[i][j][k]*vr[2*step+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: }
187: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,void *ctx)
188: {
189: Aircraft actx = (Aircraft)ctx;
190: PetscScalar *r;
191: PetscReal dx,dy;
192: const PetscScalar *u;
193: PetscErrorCode ierr;
196: VecGetArrayRead(U,&u);
197: VecGetArray(R,&r);
198: dx = u[0] - actx->lv*t*PetscCosReal(actx->lw);
199: dy = u[1] - actx->lv*t*PetscSinReal(actx->lw);
200: r[0] = dx*dx+dy*dy;
201: VecRestoreArray(R,&r);
202: VecRestoreArrayRead(U,&u);
203: return(0);
204: }
206: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,void *ctx)
207: {
208: Aircraft actx = (Aircraft)ctx;
209: PetscScalar drdu[2][1];
210: const PetscScalar *u;
211: PetscReal dx,dy;
212: PetscInt row[] = {0,1},col[] = {0};
213: PetscErrorCode ierr;
216: VecGetArrayRead(U,&u);
217: dx = u[0] - actx->lv*t*PetscCosReal(actx->lw);
218: dy = u[1] - actx->lv*t*PetscSinReal(actx->lw);
219: drdu[0][0] = 2.*dx;
220: drdu[1][0] = 2.*dy;
221: MatSetValues(DRDU,2,row,1,col,&drdu[0][0],INSERT_VALUES);
222: VecRestoreArrayRead(U,&u);
223: MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
224: MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
225: return(0);
226: }
228: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,void *ctx)
229: {
233: MatZeroEntries(DRDP);
234: MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
235: MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
236: return(0);
237: }
239: int main(int argc,char **argv)
240: {
241: Vec P,PL,PU;
242: struct _n_aircraft aircraft;
243: PetscMPIInt size;
244: Tao tao;
245: KSP ksp;
246: PC pc;
247: PetscScalar *u,*p;
248: PetscInt i;
249: PetscErrorCode ierr;
251: /* Initialize program */
252: PetscInitialize(&argc,&argv,NULL,NULL);if (ierr) return ierr;
253: MPI_Comm_size(PETSC_COMM_WORLD,&size);
254: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
256: /* Parameter settings */
257: aircraft.ftime = 1.; /* time interval in hour */
258: aircraft.nsteps = 10; /* number of steps */
259: aircraft.lv = 2.0; /* leader speed in kmph */
260: aircraft.lw = PETSC_PI/4.; /* leader heading angle */
262: PetscOptionsGetReal(NULL,NULL,"-ftime",&aircraft.ftime,NULL);
263: PetscOptionsGetInt(NULL,NULL,"-nsteps",&aircraft.nsteps,NULL);
264: PetscOptionsHasName(NULL,NULL,"-matrixfree",&aircraft.mf);
265: PetscOptionsHasName(NULL,NULL,"-exacthessian",&aircraft.eh);
267: /* Create TAO solver and set desired solution method */
268: TaoCreate(PETSC_COMM_WORLD,&tao);
269: TaoSetType(tao,TAOBQNLS);
271: /* Create necessary matrix and vectors, solve same ODE on every process */
272: MatCreate(PETSC_COMM_WORLD,&aircraft.A);
273: MatSetSizes(aircraft.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
274: MatSetFromOptions(aircraft.A);
275: MatSetUp(aircraft.A);
276: MatAssemblyBegin(aircraft.A,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd(aircraft.A,MAT_FINAL_ASSEMBLY);
278: MatShift(aircraft.A,1);
279: MatShift(aircraft.A,-1);
281: MatCreate(PETSC_COMM_WORLD,&aircraft.Jacp);
282: MatSetSizes(aircraft.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2*aircraft.nsteps);
283: MatSetFromOptions(aircraft.Jacp);
284: MatSetUp(aircraft.Jacp);
285: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,1,NULL,&aircraft.DRDP);
286: MatSetUp(aircraft.DRDP);
287: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&aircraft.DRDU);
288: MatSetUp(aircraft.DRDU);
290: /* Create timestepping solver context */
291: TSCreate(PETSC_COMM_WORLD,&aircraft.ts);
292: TSSetType(aircraft.ts,TSRK);
293: TSSetRHSFunction(aircraft.ts,NULL,RHSFunction,&aircraft);
294: TSSetRHSJacobian(aircraft.ts,aircraft.A,aircraft.A,TSComputeRHSJacobianConstant,&aircraft);
295: TSSetRHSJacobianP(aircraft.ts,aircraft.Jacp,RHSJacobianP,&aircraft);
296: TSSetExactFinalTime(aircraft.ts,TS_EXACTFINALTIME_MATCHSTEP);
297: TSSetEquationType(aircraft.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
299: /* Set initial conditions */
300: MatCreateVecs(aircraft.A,&aircraft.U,NULL);
301: TSSetSolution(aircraft.ts,aircraft.U);
302: VecGetArray(aircraft.U,&u);
303: u[0] = 1.5;
304: u[1] = 0;
305: VecRestoreArray(aircraft.U,&u);
306: VecCreate(PETSC_COMM_WORLD,&aircraft.V);
307: VecSetSizes(aircraft.V,PETSC_DECIDE,aircraft.nsteps);
308: VecSetUp(aircraft.V);
309: VecDuplicate(aircraft.V,&aircraft.W);
310: VecSet(aircraft.V,1.);
311: VecSet(aircraft.W,PETSC_PI/4.);
313: /* Save trajectory of solution so that TSAdjointSolve() may be used */
314: TSSetSaveTrajectory(aircraft.ts);
316: /* Set sensitivity context */
317: TSCreateQuadratureTS(aircraft.ts,PETSC_FALSE,&aircraft.quadts);
318: TSSetRHSFunction(aircraft.quadts,NULL,(TSRHSFunction)CostIntegrand,&aircraft);
319: TSSetRHSJacobian(aircraft.quadts,aircraft.DRDU,aircraft.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&aircraft);
320: TSSetRHSJacobianP(aircraft.quadts,aircraft.DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&aircraft);
321: MatCreateVecs(aircraft.A,&aircraft.Lambda[0],NULL);
322: MatCreateVecs(aircraft.Jacp,&aircraft.Mup[0],NULL);
323: if (aircraft.eh) {
324: MatCreateVecs(aircraft.A,&aircraft.rhshp1[0],NULL);
325: MatCreateVecs(aircraft.A,&aircraft.rhshp2[0],NULL);
326: MatCreateVecs(aircraft.Jacp,&aircraft.rhshp3[0],NULL);
327: MatCreateVecs(aircraft.Jacp,&aircraft.rhshp4[0],NULL);
328: MatCreateVecs(aircraft.DRDU,&aircraft.inthp1[0],NULL);
329: MatCreateVecs(aircraft.DRDU,&aircraft.inthp2[0],NULL);
330: MatCreateVecs(aircraft.DRDP,&aircraft.inthp3[0],NULL);
331: MatCreateVecs(aircraft.DRDP,&aircraft.inthp4[0],NULL);
332: MatCreateVecs(aircraft.Jacp,&aircraft.Dir,NULL);
333: TSSetRHSHessianProduct(aircraft.ts,aircraft.rhshp1,RHSHessianProductUU,aircraft.rhshp2,RHSHessianProductUP,aircraft.rhshp3,RHSHessianProductPU,aircraft.rhshp4,RHSHessianProductPP,&aircraft);
334: TSSetRHSHessianProduct(aircraft.quadts,aircraft.inthp1,IntegrandHessianProductUU,aircraft.inthp2,IntegrandHessianProductUP,aircraft.inthp3,IntegrandHessianProductPU,aircraft.inthp4,IntegrandHessianProductPP,&aircraft);
335: MatCreateVecs(aircraft.A,&aircraft.Lambda2[0],NULL);
336: MatCreateVecs(aircraft.Jacp,&aircraft.Mup2[0],NULL);
337: }
338: TSSetFromOptions(aircraft.ts);
339: TSSetMaxTime(aircraft.ts,aircraft.ftime);
340: TSSetTimeStep(aircraft.ts,aircraft.ftime/aircraft.nsteps);
342: /* Set initial solution guess */
343: MatCreateVecs(aircraft.Jacp,&P,NULL);
344: VecGetArray(P,&p);
345: for (i=0; i<aircraft.nsteps; i++) {
346: p[2*i] = 2.0;
347: p[2*i+1] = PETSC_PI/2.0;
348: }
349: VecRestoreArray(P,&p);
350: VecDuplicate(P,&PU);
351: VecDuplicate(P,&PL);
352: VecGetArray(PU,&p);
353: for (i=0; i<aircraft.nsteps; i++) {
354: p[2*i] = 2.0;
355: p[2*i+1] = PETSC_PI;
356: }
357: VecRestoreArray(PU,&p);
358: VecGetArray(PL,&p);
359: for (i=0; i<aircraft.nsteps; i++) {
360: p[2*i] = 0.0;
361: p[2*i+1] = -PETSC_PI;
362: }
363: VecRestoreArray(PL,&p);
365: TaoSetInitialVector(tao,P);
366: TaoSetVariableBounds(tao,PL,PU);
367: /* Set routine for function and gradient evaluation */
368: TaoSetObjectiveAndGradientRoutine(tao,FormObjFunctionGradient,(void *)&aircraft);
370: if (aircraft.eh) {
371: if (aircraft.mf) {
372: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,2*aircraft.nsteps,(void*)&aircraft,&aircraft.H);
373: MatShellSetOperation(aircraft.H,MATOP_MULT,(void(*)(void))MyMatMult);
374: MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
375: TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,MatrixFreeObjHessian,(void*)&aircraft);
376: } else {
377: MatCreateDense(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,2*aircraft.nsteps,2*aircraft.nsteps,NULL,&(aircraft.H));
378: MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
379: TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,FormObjHessian,(void *)&aircraft);
380: }
381: }
383: /* Check for any TAO command line options */
384: TaoGetKSP(tao,&ksp);
385: if (ksp) {
386: KSPGetPC(ksp,&pc);
387: PCSetType(pc,PCNONE);
388: }
389: TaoSetFromOptions(tao);
391: TaoSolve(tao);
392: VecView(P,PETSC_VIEWER_STDOUT_WORLD);
394: /* Free TAO data structures */
395: TaoDestroy(&tao);
397: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
398: Free work space. All PETSc objects should be destroyed when they
399: are no longer needed.
400: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
401: TSDestroy(&aircraft.ts);
402: MatDestroy(&aircraft.A);
403: VecDestroy(&aircraft.U);
404: VecDestroy(&aircraft.V);
405: VecDestroy(&aircraft.W);
406: VecDestroy(&P);
407: VecDestroy(&PU);
408: VecDestroy(&PL);
409: MatDestroy(&aircraft.Jacp);
410: MatDestroy(&aircraft.DRDU);
411: MatDestroy(&aircraft.DRDP);
412: VecDestroy(&aircraft.Lambda[0]);
413: VecDestroy(&aircraft.Mup[0]);
414: VecDestroy(&P);
415: if (aircraft.eh) {
416: VecDestroy(&aircraft.Lambda2[0]);
417: VecDestroy(&aircraft.Mup2[0]);
418: VecDestroy(&aircraft.Dir);
419: VecDestroy(&aircraft.rhshp1[0]);
420: VecDestroy(&aircraft.rhshp2[0]);
421: VecDestroy(&aircraft.rhshp3[0]);
422: VecDestroy(&aircraft.rhshp4[0]);
423: VecDestroy(&aircraft.inthp1[0]);
424: VecDestroy(&aircraft.inthp2[0]);
425: VecDestroy(&aircraft.inthp3[0]);
426: VecDestroy(&aircraft.inthp4[0]);
427: MatDestroy(&aircraft.H);
428: }
429: PetscFinalize();
430: return ierr;
431: }
433: /*
434: FormObjFunctionGradient - Evaluates the function and corresponding gradient.
436: Input Parameters:
437: tao - the Tao context
438: P - the input vector
439: ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradientRoutine()
441: Output Parameters:
442: f - the newly evaluated function
443: G - the newly evaluated gradient
444: */
445: PetscErrorCode FormObjFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
446: {
447: Aircraft actx = (Aircraft)ctx;
448: TS ts = actx->ts;
449: Vec Q;
450: const PetscScalar *p,*q;
451: PetscScalar *u,*v,*w;
452: PetscInt i;
453: PetscErrorCode ierr;
456: VecGetArrayRead(P,&p);
457: VecGetArray(actx->V,&v);
458: VecGetArray(actx->W,&w);
459: for (i=0; i<actx->nsteps; i++) {
460: v[i] = p[2*i];
461: w[i] = p[2*i+1];
462: }
463: VecRestoreArrayRead(P,&p);
464: VecRestoreArray(actx->V,&v);
465: VecRestoreArray(actx->W,&w);
467: TSSetTime(ts,0.0);
468: TSSetStepNumber(ts,0);
469: TSSetFromOptions(ts);
470: TSSetTimeStep(ts,actx->ftime/actx->nsteps);
472: /* reinitialize system state */
473: VecGetArray(actx->U,&u);
474: u[0] = 2.0;
475: u[1] = 0;
476: VecRestoreArray(actx->U,&u);
478: /* reinitialize the integral value */
479: TSGetCostIntegral(ts,&Q);
480: VecSet(Q,0.0);
482: TSSolve(ts,actx->U);
484: /* Reset initial conditions for the adjoint integration */
485: VecSet(actx->Lambda[0],0.0);
486: VecSet(actx->Mup[0],0.0);
487: TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);
489: TSAdjointSolve(ts);
490: VecCopy(actx->Mup[0],G);
491: TSGetCostIntegral(ts,&Q);
492: VecGetArrayRead(Q,&q);
493: *f = q[0];
494: VecRestoreArrayRead(Q,&q);
495: return(0);
496: }
498: PetscErrorCode FormObjHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
499: {
500: Aircraft actx = (Aircraft)ctx;
501: const PetscScalar *p;
502: PetscScalar *harr,*v,*w,one = 1.0;
503: PetscInt ind[1];
504: PetscInt *cols,i;
505: Vec Dir;
506: PetscErrorCode ierr;
509: /* set up control parameters */
510: VecGetArrayRead(P,&p);
511: VecGetArray(actx->V,&v);
512: VecGetArray(actx->W,&w);
513: for (i=0; i<actx->nsteps; i++) {
514: v[i] = p[2*i];
515: w[i] = p[2*i+1];
516: }
517: VecRestoreArrayRead(P,&p);
518: VecRestoreArray(actx->V,&v);
519: VecRestoreArray(actx->W,&w);
521: PetscMalloc1(2*actx->nsteps,&harr);
522: PetscMalloc1(2*actx->nsteps,&cols);
523: for (i=0; i<2*actx->nsteps; i++) cols[i] = i;
524: VecDuplicate(P,&Dir);
525: for (i=0; i<2*actx->nsteps; i++) {
526: ind[0] = i;
527: VecSet(Dir,0.0);
528: VecSetValues(Dir,1,ind,&one,INSERT_VALUES);
529: VecAssemblyBegin(Dir);
530: VecAssemblyEnd(Dir);
531: ComputeObjHessianWithSOA(Dir,harr,actx);
532: MatSetValues(H,1,ind,2*actx->nsteps,cols,harr,INSERT_VALUES);
533: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
534: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
535: if (H != Hpre) {
536: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
537: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
538: }
539: }
540: PetscFree(cols);
541: PetscFree(harr);
542: VecDestroy(&Dir);
543: return(0);
544: }
546: PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
547: {
548: Aircraft actx = (Aircraft)ctx;
549: PetscScalar *v,*w;
550: const PetscScalar *p;
551: PetscInt i;
552: PetscErrorCode ierr;
555: VecGetArrayRead(P,&p);
556: VecGetArray(actx->V,&v);
557: VecGetArray(actx->W,&w);
558: for (i=0; i<actx->nsteps; i++) {
559: v[i] = p[2*i];
560: w[i] = p[2*i+1];
561: }
562: VecRestoreArrayRead(P,&p);
563: VecRestoreArray(actx->V,&v);
564: VecRestoreArray(actx->W,&w);
565: return(0);
566: }
568: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
569: {
570: PetscScalar *y;
571: void *ptr;
575: MatShellGetContext(H_shell,&ptr);
576: VecGetArray(Y,&y);
577: ComputeObjHessianWithSOA(X,y,(Aircraft)ptr);
578: VecRestoreArray(Y,&y);
579: return(0);
580: }
582: PetscErrorCode ComputeObjHessianWithSOA(Vec Dir,PetscScalar arr[],Aircraft actx)
583: {
584: TS ts = actx->ts;
585: const PetscScalar *z_ptr;
586: PetscScalar *u;
587: Vec Q;
588: PetscInt i;
589: PetscErrorCode ierr;
592: /* Reset TSAdjoint so that AdjointSetUp will be called again */
593: TSAdjointReset(ts);
595: TSSetTime(ts,0.0);
596: TSSetStepNumber(ts,0);
597: TSSetFromOptions(ts);
598: TSSetTimeStep(ts,actx->ftime/actx->nsteps);
599: TSSetCostHessianProducts(actx->ts,1,actx->Lambda2,actx->Mup2,Dir);
601: /* reinitialize system state */
602: VecGetArray(actx->U,&u);
603: u[0] = 2.0;
604: u[1] = 0;
605: VecRestoreArray(actx->U,&u);
607: /* reinitialize the integral value */
608: TSGetCostIntegral(ts,&Q);
609: VecSet(Q,0.0);
611: /* initialize tlm variable */
612: MatZeroEntries(actx->Jacp);
613: MatAssemblyBegin(actx->Jacp,MAT_FINAL_ASSEMBLY);
614: MatAssemblyEnd(actx->Jacp,MAT_FINAL_ASSEMBLY);
615: TSAdjointSetForward(ts,actx->Jacp);
617: TSSolve(ts,actx->U);
619: /* Set terminal conditions for first- and second-order adjonts */
620: VecSet(actx->Lambda[0],0.0);
621: VecSet(actx->Mup[0],0.0);
622: VecSet(actx->Lambda2[0],0.0);
623: VecSet(actx->Mup2[0],0.0);
624: TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);
626: TSGetCostIntegral(ts,&Q);
628: /* Reset initial conditions for the adjoint integration */
629: TSAdjointSolve(ts);
631: /* initial condition does not depend on p, so that lambda is not needed to assemble G */
632: VecGetArrayRead(actx->Mup2[0],&z_ptr);
633: for (i=0; i<2*actx->nsteps; i++) arr[i] = z_ptr[i];
634: VecRestoreArrayRead(actx->Mup2[0],&z_ptr);
636: /* Disable second-order adjoint mode */
637: TSAdjointReset(ts);
638: TSAdjointResetForward(ts);
639: return(0);
640: }
642: /*TEST
643: build:
644: requires: !complex !single
646: test:
647: args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_gatol 1e-7
649: test:
650: suffix: 2
651: 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
653: test:
654: suffix: 3
655: 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
656: TEST*/