Actual source code: ex20opt_ic.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n";
3: /**
4: Concepts: TS^time-dependent nonlinear problems
5: Concepts: TS^van der Pol equation DAE equivalent
6: Concepts: Optimization using adjoint sensitivities
7: Processors: 1
8: */
9: /**
10: Notes:
11: This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS.
12: The nonlinear problem is written in an ODE equivalent form.
13: The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions.
14: The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details.
15: */
17: #include <petsctao.h>
18: #include <petscts.h>
20: typedef struct _n_User *User;
21: struct _n_User {
22: TS ts;
23: PetscReal mu;
24: PetscReal next_output;
26: /* Sensitivity analysis support */
27: PetscInt steps;
28: PetscReal ftime;
29: Mat A; /* Jacobian matrix for ODE */
30: Mat Jacp; /* JacobianP matrix for ODE*/
31: Mat H; /* Hessian matrix for optimization */
32: Vec U,Lambda[1],Mup[1]; /* first-order adjoint variables */
33: Vec Lambda2[2]; /* second-order adjoint variables */
34: Vec Ihp1[1]; /* working space for Hessian evaluations */
35: Vec Dir; /* direction vector */
36: PetscReal ob[2]; /* observation used by the cost function */
37: PetscBool implicitform; /* implicit ODE? */
38: };
39: PetscErrorCode Adjoint2(Vec,PetscScalar[],User);
41: /* ----------------------- Explicit form of the ODE -------------------- */
43: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
44: {
45: PetscErrorCode ierr;
46: User user = (User)ctx;
47: PetscScalar *f;
48: const PetscScalar *u;
51: VecGetArrayRead(U,&u);
52: VecGetArray(F,&f);
53: f[0] = u[1];
54: f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
55: VecRestoreArrayRead(U,&u);
56: VecRestoreArray(F,&f);
57: return(0);
58: }
60: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
61: {
62: PetscErrorCode ierr;
63: User user = (User)ctx;
64: PetscReal mu = user->mu;
65: PetscInt rowcol[] = {0,1};
66: PetscScalar J[2][2];
67: const PetscScalar *u;
70: VecGetArrayRead(U,&u);
71: J[0][0] = 0;
72: J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
73: J[0][1] = 1.0;
74: J[1][1] = mu*(1.0-u[0]*u[0]);
75: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
76: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
78: if (A != B) {
79: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
81: }
82: VecRestoreArrayRead(U,&u);
83: return(0);
84: }
86: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
87: {
88: const PetscScalar *vl,*vr,*u;
89: PetscScalar *vhv;
90: PetscScalar dJdU[2][2][2]={{{0}}};
91: PetscInt i,j,k;
92: User user = (User)ctx;
93: PetscErrorCode ierr;
96: VecGetArrayRead(U,&u);
97: VecGetArrayRead(Vl[0],&vl);
98: VecGetArrayRead(Vr,&vr);
99: VecGetArray(VHV[0],&vhv);
101: dJdU[1][0][0] = -2.*user->mu*u[1];
102: dJdU[1][1][0] = -2.*user->mu*u[0];
103: dJdU[1][0][1] = -2.*user->mu*u[0];
104: for (j=0;j<2;j++) {
105: vhv[j] = 0;
106: for (k=0;k<2;k++)
107: for (i=0;i<2;i++ )
108: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
109: }
110: VecRestoreArrayRead(U,&u);
111: VecRestoreArrayRead(Vl[0],&vl);
112: VecRestoreArrayRead(Vr,&vr);
113: VecRestoreArray(VHV[0],&vhv);
114: return(0);
115: }
117: /* ----------------------- Implicit form of the ODE -------------------- */
119: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
120: {
121: PetscErrorCode ierr;
122: User user = (User)ctx;
123: const PetscScalar *u,*udot;
124: PetscScalar *f;
127: VecGetArrayRead(U,&u);
128: VecGetArrayRead(Udot,&udot);
129: VecGetArray(F,&f);
130: f[0] = udot[0] - u[1];
131: f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
132: VecRestoreArrayRead(U,&u);
133: VecRestoreArrayRead(Udot,&udot);
134: VecRestoreArray(F,&f);
135: return(0);
136: }
138: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
139: {
140: PetscErrorCode ierr;
141: User user = (User)ctx;
142: PetscInt rowcol[] = {0,1};
143: PetscScalar J[2][2];
144: const PetscScalar *u;
147: VecGetArrayRead(U,&u);
148: J[0][0] = a; J[0][1] = -1.0;
149: J[1][0] = user->mu*(1.0 + 2.0*u[0]*u[1]); J[1][1] = a - user->mu*(1.0-u[0]*u[0]);
150: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
151: VecRestoreArrayRead(U,&u);
152: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
153: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
154: if (A != B) {
155: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
156: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
157: }
158: return(0);
159: }
161: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
162: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
163: {
164: PetscErrorCode ierr;
165: const PetscScalar *u;
166: PetscReal tfinal, dt;
167: User user = (User)ctx;
168: Vec interpolatedU;
171: TSGetTimeStep(ts,&dt);
172: TSGetMaxTime(ts,&tfinal);
174: while (user->next_output <= t && user->next_output <= tfinal) {
175: VecDuplicate(U,&interpolatedU);
176: TSInterpolate(ts,user->next_output,interpolatedU);
177: VecGetArrayRead(interpolatedU,&u);
178: PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
179: (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
180: (double)PetscRealPart(u[1]));
181: VecRestoreArrayRead(interpolatedU,&u);
182: VecDestroy(&interpolatedU);
183: user->next_output += 0.1;
184: }
185: return(0);
186: }
188: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
189: {
190: const PetscScalar *vl,*vr,*u;
191: PetscScalar *vhv;
192: PetscScalar dJdU[2][2][2]={{{0}}};
193: PetscInt i,j,k;
194: User user = (User)ctx;
195: PetscErrorCode ierr;
198: VecGetArrayRead(U,&u);
199: VecGetArrayRead(Vl[0],&vl);
200: VecGetArrayRead(Vr,&vr);
201: VecGetArray(VHV[0],&vhv);
202: dJdU[1][0][0] = 2.*user->mu*u[1];
203: dJdU[1][1][0] = 2.*user->mu*u[0];
204: dJdU[1][0][1] = 2.*user->mu*u[0];
205: for (j=0;j<2;j++) {
206: vhv[j] = 0;
207: for (k=0;k<2;k++)
208: for (i=0;i<2;i++ )
209: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
210: }
211: VecRestoreArrayRead(U,&u);
212: VecRestoreArrayRead(Vl[0],&vl);
213: VecRestoreArrayRead(Vr,&vr);
214: VecRestoreArray(VHV[0],&vhv);
215: return(0);
216: }
218: /* ------------------ User-defined routines for TAO -------------------------- */
220: static PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx)
221: {
222: User user_ptr = (User)ctx;
223: TS ts = user_ptr->ts;
224: const PetscScalar *x_ptr;
225: PetscScalar *y_ptr;
229: VecCopy(IC,user_ptr->U); /* set up the initial condition */
231: TSSetTime(ts,0.0);
232: TSSetStepNumber(ts,0);
233: TSSetTimeStep(ts,0.001); /* can be overwritten by command line options */
234: TSSetFromOptions(ts);
235: TSSolve(ts,user_ptr->U);
237: VecGetArrayRead(user_ptr->U,&x_ptr);
238: VecGetArray(user_ptr->Lambda[0],&y_ptr);
239: *f = (x_ptr[0]-user_ptr->ob[0])*(x_ptr[0]-user_ptr->ob[0])+(x_ptr[1]-user_ptr->ob[1])*(x_ptr[1]-user_ptr->ob[1]);
240: y_ptr[0] = 2.*(x_ptr[0]-user_ptr->ob[0]);
241: y_ptr[1] = 2.*(x_ptr[1]-user_ptr->ob[1]);
242: VecRestoreArray(user_ptr->Lambda[0],&y_ptr);
243: VecRestoreArrayRead(user_ptr->U,&x_ptr);
245: TSSetCostGradients(ts,1,user_ptr->Lambda,NULL);
246: TSAdjointSolve(ts);
247: VecCopy(user_ptr->Lambda[0],G);
248: return(0);
249: }
251: static PetscErrorCode FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
252: {
253: User user_ptr = (User)ctx;
254: PetscScalar harr[2];
255: PetscScalar *x_ptr;
256: const PetscInt rows[2] = {0,1};
257: PetscInt col;
261: VecCopy(U,user_ptr->U);
262: VecGetArray(user_ptr->Dir,&x_ptr);
263: x_ptr[0] = 1.;
264: x_ptr[1] = 0.;
265: VecRestoreArray(user_ptr->Dir,&x_ptr);
266: Adjoint2(user_ptr->U,harr,user_ptr);
267: col = 0;
268: MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);
270: VecCopy(U,user_ptr->U);
271: VecGetArray(user_ptr->Dir,&x_ptr);
272: x_ptr[0] = 0.;
273: x_ptr[1] = 1.;
274: VecRestoreArray(user_ptr->Dir,&x_ptr);
275: Adjoint2(user_ptr->U,harr,user_ptr);
276: col = 1;
277: MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);
279: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
280: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
281: if (H != Hpre) {
282: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
283: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
284: }
285: return(0);
286: }
288: static PetscErrorCode MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
289: {
290: User user_ptr = (User)ctx;
294: VecCopy(U,user_ptr->U);
295: return(0);
296: }
298: /* ------------ Routines calculating second-order derivatives -------------- */
300: /**
301: Compute the Hessian-vector product for the cost function using Second-order adjoint
302: */
303: PetscErrorCode Adjoint2(Vec U,PetscScalar arr[],User ctx)
304: {
305: TS ts = ctx->ts;
306: PetscScalar *x_ptr,*y_ptr;
307: Mat tlmsen;
311: TSAdjointReset(ts);
313: TSSetTime(ts,0.0);
314: TSSetStepNumber(ts,0);
315: TSSetTimeStep(ts,0.001);
316: TSSetFromOptions(ts);
317: TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir);
318: TSAdjointSetForward(ts,NULL);
319: TSSolve(ts,U);
321: /* Set terminal conditions for first- and second-order adjonts */
322: VecGetArray(U,&x_ptr);
323: VecGetArray(ctx->Lambda[0],&y_ptr);
324: y_ptr[0] = 2.*(x_ptr[0]-ctx->ob[0]);
325: y_ptr[1] = 2.*(x_ptr[1]-ctx->ob[1]);
326: VecRestoreArray(ctx->Lambda[0],&y_ptr);
327: VecRestoreArray(U,&x_ptr);
328: TSForwardGetSensitivities(ts,NULL,&tlmsen);
329: MatDenseGetColumn(tlmsen,0,&x_ptr);
330: VecGetArray(ctx->Lambda2[0],&y_ptr);
331: y_ptr[0] = 2.*x_ptr[0];
332: y_ptr[1] = 2.*x_ptr[1];
333: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
334: MatDenseRestoreColumn(tlmsen,&x_ptr);
336: TSSetCostGradients(ts,1,ctx->Lambda,NULL);
337: if (ctx->implicitform) {
338: TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
339: } else {
340: TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
341: }
342: TSAdjointSolve(ts);
344: VecGetArray(ctx->Lambda2[0],&x_ptr);
345: arr[0] = x_ptr[0];
346: arr[1] = x_ptr[1];
347: VecRestoreArray(ctx->Lambda2[0],&x_ptr);
349: TSAdjointReset(ts);
350: TSAdjointResetForward(ts);
351: return(0);
352: }
354: PetscErrorCode FiniteDiff(Vec U,PetscScalar arr[],User ctx)
355: {
356: Vec Up,G,Gp;
357: const PetscScalar eps = PetscRealConstant(1e-7);
358: PetscScalar *u;
359: Tao tao = NULL;
360: PetscReal f;
361: PetscErrorCode ierr;
364: VecDuplicate(U,&Up);
365: VecDuplicate(U,&G);
366: VecDuplicate(U,&Gp);
368: FormFunctionGradient(tao,U,&f,G,ctx);
370: VecCopy(U,Up);
371: VecGetArray(Up,&u);
372: u[0] += eps;
373: VecRestoreArray(Up,&u);
374: FormFunctionGradient(tao,Up,&f,Gp,ctx);
375: VecAXPY(Gp,-1,G);
376: VecScale(Gp,1./eps);
377: VecGetArray(Gp,&u);
378: arr[0] = u[0];
379: arr[1] = u[1];
380: VecRestoreArray(Gp,&u);
382: VecCopy(U,Up);
383: VecGetArray(Up,&u);
384: u[1] += eps;
385: VecRestoreArray(Up,&u);
386: FormFunctionGradient(tao,Up,&f,Gp,ctx);
387: VecAXPY(Gp,-1,G);
388: VecScale(Gp,1./eps);
389: VecGetArray(Gp,&u);
390: arr[2] = u[0];
391: arr[3] = u[1];
392: VecRestoreArray(Gp,&u);
394: VecDestroy(&G);
395: VecDestroy(&Gp);
396: VecDestroy(&Up);
397: return(0);
398: }
400: static PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
401: {
402: User user_ptr;
403: PetscScalar *y_ptr;
407: MatShellGetContext(mat,(void*)&user_ptr);
408: VecCopy(svec,user_ptr->Dir);
409: VecGetArray(y,&y_ptr);
410: Adjoint2(user_ptr->U,y_ptr,user_ptr);
411: VecRestoreArray(y,&y_ptr);
412: return(0);
413: }
415: int main(int argc,char **argv)
416: {
417: PetscBool monitor = PETSC_FALSE,mf = PETSC_TRUE;
418: PetscInt mode = 0;
419: PetscMPIInt size;
420: struct _n_User user;
421: Vec x; /* working vector for TAO */
422: PetscScalar *x_ptr,arr[4];
423: PetscScalar ic1 = 2.2,ic2 = -0.7; /* initial guess for TAO */
424: Tao tao;
425: KSP ksp;
426: PC pc;
429: /* Initialize program */
430: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
431: MPI_Comm_size(PETSC_COMM_WORLD,&size);
432: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
434: /* Set runtime options */
435: user.next_output = 0.0;
436: user.mu = 1.0e3;
437: user.steps = 0;
438: user.ftime = 0.5;
439: user.implicitform = PETSC_TRUE;
441: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
442: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
443: PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL);
444: PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL);
445: PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL);
446: PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL); /* matrix-free hessian for optimization */
447: PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);
449: /* Create necessary matrix and vectors, solve same ODE on every process */
450: MatCreate(PETSC_COMM_WORLD,&user.A);
451: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
452: MatSetFromOptions(user.A);
453: MatSetUp(user.A);
454: MatCreateVecs(user.A,&user.U,NULL);
455: MatCreateVecs(user.A,&user.Dir,NULL);
456: MatCreateVecs(user.A,&user.Lambda[0],NULL);
457: MatCreateVecs(user.A,&user.Lambda2[0],NULL);
458: MatCreateVecs(user.A,&user.Ihp1[0],NULL);
460: /* Create timestepping solver context */
461: TSCreate(PETSC_COMM_WORLD,&user.ts);
462: TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
463: if (user.implicitform) {
464: TSSetIFunction(user.ts,NULL,IFunction,&user);
465: TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
466: TSSetType(user.ts,TSCN);
467: } else {
468: TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
469: TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
470: TSSetType(user.ts,TSRK);
471: }
472: TSSetMaxTime(user.ts,user.ftime);
473: TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);
475: if (monitor) {
476: TSMonitorSet(user.ts,Monitor,&user,NULL);
477: }
479: /* Set ODE initial conditions */
480: VecGetArray(user.U,&x_ptr);
481: x_ptr[0] = 2.0;
482: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
483: VecRestoreArray(user.U,&x_ptr);
485: /* Set runtime options */
486: TSSetFromOptions(user.ts);
488: /* Obtain the observation */
489: TSSolve(user.ts,user.U);
490: VecGetArray(user.U,&x_ptr);
491: user.ob[0] = x_ptr[0];
492: user.ob[1] = x_ptr[1];
493: VecRestoreArray(user.U,&x_ptr);
495: VecDuplicate(user.U,&x);
496: VecGetArray(x,&x_ptr);
497: x_ptr[0] = ic1;
498: x_ptr[1] = ic2;
499: VecRestoreArray(x,&x_ptr);
501: /* Save trajectory of solution so that TSAdjointSolve() may be used */
502: TSSetSaveTrajectory(user.ts);
504: /* Compare finite difference and second-order adjoint. */
505: switch (mode) {
506: case 2 :
507: FiniteDiff(x,arr,&user);
508: PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n");
509: PetscPrintf(PETSC_COMM_WORLD,"%g %g\n%g %g\n",(double)arr[0],(double)arr[1],(double)arr[2],(double)arr[3]);
510: break;
511: case 1 : /* Compute the Hessian column by column */
512: VecCopy(x,user.U);
513: VecGetArray(user.Dir,&x_ptr);
514: x_ptr[0] = 1.;
515: x_ptr[1] = 0.;
516: VecRestoreArray(user.Dir,&x_ptr);
517: Adjoint2(user.U,arr,&user);
518: PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n");
519: PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
520: VecCopy(x,user.U);
521: VecGetArray(user.Dir,&x_ptr);
522: x_ptr[0] = 0.;
523: x_ptr[1] = 1.;
524: VecRestoreArray(user.Dir,&x_ptr);
525: Adjoint2(user.U,arr,&user);
526: PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n");
527: PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
528: break;
529: case 0 :
530: /* Create optimization context and set up */
531: TaoCreate(PETSC_COMM_WORLD,&tao);
532: TaoSetType(tao,TAOBLMVM);
533: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
535: if (mf) {
536: MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H);
537: MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat);
538: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
539: TaoSetHessianRoutine(tao,user.H,user.H,MatrixFreeHessian,(void *)&user);
540: } else { /* Create Hessian matrix */
541: MatCreate(PETSC_COMM_WORLD,&user.H);
542: MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2);
543: MatSetUp(user.H);
544: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
545: }
547: /* Not use any preconditioner */
548: TaoGetKSP(tao,&ksp);
549: if (ksp) {
550: KSPGetPC(ksp,&pc);
551: PCSetType(pc,PCNONE);
552: }
554: /* Set initial solution guess */
555: TaoSetInitialVector(tao,x);
556: TaoSetFromOptions(tao);
557: TaoSolve(tao);
558: TaoDestroy(&tao);
559: MatDestroy(&user.H);
560: break;
561: default:
562: break;
563: }
565: /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */
566: MatDestroy(&user.A);
567: VecDestroy(&user.U);
568: VecDestroy(&user.Lambda[0]);
569: VecDestroy(&user.Lambda2[0]);
570: VecDestroy(&user.Ihp1[0]);
571: VecDestroy(&user.Dir);
572: TSDestroy(&user.ts);
573: VecDestroy(&x);
574: PetscFinalize();
575: return(ierr);
576: }
578: /*TEST
579: build:
580: requires: !complex !single
582: test:
583: args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125
584: output_file: output/ex20opt_ic_1.out
586: test:
587: suffix: 2
588: args: -ts_type beuler -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none
589: output_file: output/ex20opt_ic_2.out
591: test:
592: suffix: 3
593: args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none
594: output_file: output/ex20opt_ic_3.out
596: test:
597: suffix: 4
598: args: -implicitform 0 -ts_dt 0.01 -ts_max_steps 2 -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view -mode 1 -my_tao_mf
599: TEST*/