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*/