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