Actual source code: ex20opt_ic.c

  1: static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n";

  3: /*
  4:   Notes:
  5:   This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS.
  6:   The nonlinear problem is written in an ODE equivalent form.
  7:   The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions.
  8:   The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details.
  9: */

 11: #include <petsctao.h>
 12: #include <petscts.h>

 14: typedef struct _n_User *User;
 15: struct _n_User {
 16:   TS        ts;
 17:   PetscReal mu;
 18:   PetscReal next_output;

 20:   /* Sensitivity analysis support */
 21:   PetscInt  steps;
 22:   PetscReal ftime;
 23:   Mat       A;                       /* Jacobian matrix for ODE */
 24:   Mat       Jacp;                    /* JacobianP matrix for ODE*/
 25:   Mat       H;                       /* Hessian matrix for optimization */
 26:   Vec       U,Lambda[1],Mup[1];      /* first-order adjoint variables */
 27:   Vec       Lambda2[2];              /* second-order adjoint variables */
 28:   Vec       Ihp1[1];                 /* working space for Hessian evaluations */
 29:   Vec       Dir;                     /* direction vector */
 30:   PetscReal ob[2];                   /* observation used by the cost function */
 31:   PetscBool implicitform;            /* implicit ODE? */
 32: };
 33: PetscErrorCode Adjoint2(Vec,PetscScalar[],User);

 35: /* ----------------------- Explicit form of the ODE  -------------------- */

 37: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
 38: {
 39:   User              user = (User)ctx;
 40:   PetscScalar       *f;
 41:   const PetscScalar *u;

 44:   VecGetArrayRead(U,&u);
 45:   VecGetArray(F,&f);
 46:   f[0] = u[1];
 47:   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
 48:   VecRestoreArrayRead(U,&u);
 49:   VecRestoreArray(F,&f);
 50:   return 0;
 51: }

 53: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
 54: {
 55:   User              user = (User)ctx;
 56:   PetscReal         mu   = user->mu;
 57:   PetscInt          rowcol[] = {0,1};
 58:   PetscScalar       J[2][2];
 59:   const PetscScalar *u;

 62:   VecGetArrayRead(U,&u);
 63:   J[0][0] = 0;
 64:   J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
 65:   J[0][1] = 1.0;
 66:   J[1][1] = mu*(1.0-u[0]*u[0]);
 67:   MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
 68:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 70:   if (A != B) {
 71:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 72:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 73:   }
 74:   VecRestoreArrayRead(U,&u);
 75:   return 0;
 76: }

 78: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
 79: {
 80:   const PetscScalar *vl,*vr,*u;
 81:   PetscScalar       *vhv;
 82:   PetscScalar       dJdU[2][2][2]={{{0}}};
 83:   PetscInt          i,j,k;
 84:   User              user = (User)ctx;

 87:   VecGetArrayRead(U,&u);
 88:   VecGetArrayRead(Vl[0],&vl);
 89:   VecGetArrayRead(Vr,&vr);
 90:   VecGetArray(VHV[0],&vhv);

 92:   dJdU[1][0][0] = -2.*user->mu*u[1];
 93:   dJdU[1][1][0] = -2.*user->mu*u[0];
 94:   dJdU[1][0][1] = -2.*user->mu*u[0];
 95:   for (j=0;j<2;j++) {
 96:     vhv[j] = 0;
 97:     for (k=0;k<2;k++)
 98:       for (i=0;i<2;i++)
 99:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
100:   }
101:   VecRestoreArrayRead(U,&u);
102:   VecRestoreArrayRead(Vl[0],&vl);
103:   VecRestoreArrayRead(Vr,&vr);
104:   VecRestoreArray(VHV[0],&vhv);
105:   return 0;
106: }

108: /* ----------------------- Implicit form of the ODE  -------------------- */

110: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
111: {
112:   User              user = (User)ctx;
113:   const PetscScalar *u,*udot;
114:   PetscScalar       *f;

117:   VecGetArrayRead(U,&u);
118:   VecGetArrayRead(Udot,&udot);
119:   VecGetArray(F,&f);
120:   f[0] = udot[0] - u[1];
121:   f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
122:   VecRestoreArrayRead(U,&u);
123:   VecRestoreArrayRead(Udot,&udot);
124:   VecRestoreArray(F,&f);
125:   return 0;
126: }

128: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
129: {
130:   User              user = (User)ctx;
131:   PetscInt          rowcol[] = {0,1};
132:   PetscScalar       J[2][2];
133:   const PetscScalar *u;

136:   VecGetArrayRead(U,&u);
137:   J[0][0] = a;     J[0][1] =  -1.0;
138:   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]);
139:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
140:   VecRestoreArrayRead(U,&u);
141:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
142:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
143:   if (A != B) {
144:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
145:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
146:   }
147:   return 0;
148: }

150: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
151: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
152: {
153:   const PetscScalar *u;
154:   PetscReal         tfinal, dt;
155:   User              user = (User)ctx;
156:   Vec               interpolatedU;

159:   TSGetTimeStep(ts,&dt);
160:   TSGetMaxTime(ts,&tfinal);

162:   while (user->next_output <= t && user->next_output <= tfinal) {
163:     VecDuplicate(U,&interpolatedU);
164:     TSInterpolate(ts,user->next_output,interpolatedU);
165:     VecGetArrayRead(interpolatedU,&u);
166:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
167:                         (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
168:                         (double)PetscRealPart(u[1])));
169:     VecRestoreArrayRead(interpolatedU,&u);
170:     VecDestroy(&interpolatedU);
171:     user->next_output += 0.1;
172:   }
173:   return 0;
174: }

176: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
177: {
178:   const PetscScalar *vl,*vr,*u;
179:   PetscScalar       *vhv;
180:   PetscScalar       dJdU[2][2][2]={{{0}}};
181:   PetscInt          i,j,k;
182:   User              user = (User)ctx;

185:   VecGetArrayRead(U,&u);
186:   VecGetArrayRead(Vl[0],&vl);
187:   VecGetArrayRead(Vr,&vr);
188:   VecGetArray(VHV[0],&vhv);
189:   dJdU[1][0][0] = 2.*user->mu*u[1];
190:   dJdU[1][1][0] = 2.*user->mu*u[0];
191:   dJdU[1][0][1] = 2.*user->mu*u[0];
192:   for (j=0;j<2;j++) {
193:     vhv[j] = 0;
194:     for (k=0;k<2;k++)
195:       for (i=0;i<2;i++)
196:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
197:   }
198:   VecRestoreArrayRead(U,&u);
199:   VecRestoreArrayRead(Vl[0],&vl);
200:   VecRestoreArrayRead(Vr,&vr);
201:   VecRestoreArray(VHV[0],&vhv);
202:   return 0;
203: }

205: /* ------------------ User-defined routines for TAO -------------------------- */

207: static PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx)
208: {
209:   User              user_ptr = (User)ctx;
210:   TS                ts = user_ptr->ts;
211:   const PetscScalar *x_ptr;
212:   PetscScalar       *y_ptr;

215:   VecCopy(IC,user_ptr->U); /* set up the initial condition */

217:   TSSetTime(ts,0.0);
218:   TSSetStepNumber(ts,0);
219:   TSSetTimeStep(ts,0.001); /* can be overwritten by command line options */
220:   TSSetFromOptions(ts);
221:   TSSolve(ts,user_ptr->U);

223:   VecGetArrayRead(user_ptr->U,&x_ptr);
224:   VecGetArray(user_ptr->Lambda[0],&y_ptr);
225:   *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]);
226:   y_ptr[0] = 2.*(x_ptr[0]-user_ptr->ob[0]);
227:   y_ptr[1] = 2.*(x_ptr[1]-user_ptr->ob[1]);
228:   VecRestoreArray(user_ptr->Lambda[0],&y_ptr);
229:   VecRestoreArrayRead(user_ptr->U,&x_ptr);

231:   TSSetCostGradients(ts,1,user_ptr->Lambda,NULL);
232:   TSAdjointSolve(ts);
233:   VecCopy(user_ptr->Lambda[0],G);
234:   return 0;
235: }

237: static PetscErrorCode FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
238: {
239:   User           user_ptr = (User)ctx;
240:   PetscScalar    harr[2];
241:   PetscScalar    *x_ptr;
242:   const PetscInt rows[2] = {0,1};
243:   PetscInt       col;

246:   VecCopy(U,user_ptr->U);
247:   VecGetArray(user_ptr->Dir,&x_ptr);
248:   x_ptr[0] = 1.;
249:   x_ptr[1] = 0.;
250:   VecRestoreArray(user_ptr->Dir,&x_ptr);
251:   Adjoint2(user_ptr->U,harr,user_ptr);
252:   col      = 0;
253:   MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);

255:   VecCopy(U,user_ptr->U);
256:   VecGetArray(user_ptr->Dir,&x_ptr);
257:   x_ptr[0] = 0.;
258:   x_ptr[1] = 1.;
259:   VecRestoreArray(user_ptr->Dir,&x_ptr);
260:   Adjoint2(user_ptr->U,harr,user_ptr);
261:   col      = 1;
262:   MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);

264:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
265:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
266:   if (H != Hpre) {
267:     MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
268:     MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
269:   }
270:   return 0;
271: }

273: static PetscErrorCode MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
274: {
275:   User           user_ptr = (User)ctx;

278:   VecCopy(U,user_ptr->U);
279:   return 0;
280: }

282: /* ------------ Routines calculating second-order derivatives -------------- */

284: /*
285:   Compute the Hessian-vector product for the cost function using Second-order adjoint
286: */
287: PetscErrorCode Adjoint2(Vec U,PetscScalar arr[],User ctx)
288: {
289:   TS             ts = ctx->ts;
290:   PetscScalar    *x_ptr,*y_ptr;
291:   Mat            tlmsen;

294:   TSAdjointReset(ts);

296:   TSSetTime(ts,0.0);
297:   TSSetStepNumber(ts,0);
298:   TSSetTimeStep(ts,0.001);
299:   TSSetFromOptions(ts);
300:   TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir);
301:   TSAdjointSetForward(ts,NULL);
302:   TSSolve(ts,U);

304:   /* Set terminal conditions for first- and second-order adjonts */
305:   VecGetArray(U,&x_ptr);
306:   VecGetArray(ctx->Lambda[0],&y_ptr);
307:   y_ptr[0] = 2.*(x_ptr[0]-ctx->ob[0]);
308:   y_ptr[1] = 2.*(x_ptr[1]-ctx->ob[1]);
309:   VecRestoreArray(ctx->Lambda[0],&y_ptr);
310:   VecRestoreArray(U,&x_ptr);
311:   TSForwardGetSensitivities(ts,NULL,&tlmsen);
312:   MatDenseGetColumn(tlmsen,0,&x_ptr);
313:   VecGetArray(ctx->Lambda2[0],&y_ptr);
314:   y_ptr[0] = 2.*x_ptr[0];
315:   y_ptr[1] = 2.*x_ptr[1];
316:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
317:   MatDenseRestoreColumn(tlmsen,&x_ptr);

319:   TSSetCostGradients(ts,1,ctx->Lambda,NULL);
320:   if (ctx->implicitform) {
321:     TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
322:   } else {
323:     TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
324:   }
325:   TSAdjointSolve(ts);

327:   VecGetArray(ctx->Lambda2[0],&x_ptr);
328:   arr[0] = x_ptr[0];
329:   arr[1] = x_ptr[1];
330:   VecRestoreArray(ctx->Lambda2[0],&x_ptr);

332:   TSAdjointReset(ts);
333:   TSAdjointResetForward(ts);
334:   return 0;
335: }

337: PetscErrorCode FiniteDiff(Vec U,PetscScalar arr[],User ctx)
338: {
339:   Vec               Up,G,Gp;
340:   const PetscScalar eps = PetscRealConstant(1e-7);
341:   PetscScalar       *u;
342:   Tao               tao = NULL;
343:   PetscReal         f;

346:   VecDuplicate(U,&Up);
347:   VecDuplicate(U,&G);
348:   VecDuplicate(U,&Gp);

350:   FormFunctionGradient(tao,U,&f,G,ctx);

352:   VecCopy(U,Up);
353:   VecGetArray(Up,&u);
354:   u[0] += eps;
355:   VecRestoreArray(Up,&u);
356:   FormFunctionGradient(tao,Up,&f,Gp,ctx);
357:   VecAXPY(Gp,-1,G);
358:   VecScale(Gp,1./eps);
359:   VecGetArray(Gp,&u);
360:   arr[0] = u[0];
361:   arr[1] = u[1];
362:   VecRestoreArray(Gp,&u);

364:   VecCopy(U,Up);
365:   VecGetArray(Up,&u);
366:   u[1] += eps;
367:   VecRestoreArray(Up,&u);
368:   FormFunctionGradient(tao,Up,&f,Gp,ctx);
369:   VecAXPY(Gp,-1,G);
370:   VecScale(Gp,1./eps);
371:   VecGetArray(Gp,&u);
372:   arr[2] = u[0];
373:   arr[3] = u[1];
374:   VecRestoreArray(Gp,&u);

376:   VecDestroy(&G);
377:   VecDestroy(&Gp);
378:   VecDestroy(&Up);
379:   return 0;
380: }

382: static PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
383: {
384:   User           user_ptr;
385:   PetscScalar    *y_ptr;

388:   MatShellGetContext(mat,&user_ptr);
389:   VecCopy(svec,user_ptr->Dir);
390:   VecGetArray(y,&y_ptr);
391:   Adjoint2(user_ptr->U,y_ptr,user_ptr);
392:   VecRestoreArray(y,&y_ptr);
393:   return 0;
394: }

396: int main(int argc,char **argv)
397: {
398:   PetscBool      monitor = PETSC_FALSE,mf = PETSC_TRUE;
399:   PetscInt       mode = 0;
400:   PetscMPIInt    size;
401:   struct _n_User user;
402:   Vec            x; /* working vector for TAO */
403:   PetscScalar    *x_ptr,arr[4];
404:   PetscScalar    ic1 = 2.2,ic2 = -0.7; /* initial guess for TAO */
405:   Tao            tao;
406:   KSP            ksp;
407:   PC             pc;

409:   /* Initialize program */
410:   PetscInitialize(&argc,&argv,NULL,help);
411:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

414:   /* Set runtime options */
415:   user.next_output  = 0.0;
416:   user.mu           = 1.0e3;
417:   user.steps        = 0;
418:   user.ftime        = 0.5;
419:   user.implicitform = PETSC_TRUE;

421:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
422:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
423:   PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL);
424:   PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL);
425:   PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL);
426:   PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL); /* matrix-free hessian for optimization */
427:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);

429:   /* Create necessary matrix and vectors, solve same ODE on every process */
430:   MatCreate(PETSC_COMM_WORLD,&user.A);
431:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
432:   MatSetFromOptions(user.A);
433:   MatSetUp(user.A);
434:   MatCreateVecs(user.A,&user.U,NULL);
435:   MatCreateVecs(user.A,&user.Dir,NULL);
436:   MatCreateVecs(user.A,&user.Lambda[0],NULL);
437:   MatCreateVecs(user.A,&user.Lambda2[0],NULL);
438:   MatCreateVecs(user.A,&user.Ihp1[0],NULL);

440:   /* Create timestepping solver context */
441:   TSCreate(PETSC_COMM_WORLD,&user.ts);
442:   TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
443:   if (user.implicitform) {
444:     TSSetIFunction(user.ts,NULL,IFunction,&user);
445:     TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
446:     TSSetType(user.ts,TSCN);
447:   } else {
448:     TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
449:     TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
450:     TSSetType(user.ts,TSRK);
451:   }
452:   TSSetMaxTime(user.ts,user.ftime);
453:   TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);

455:   if (monitor) {
456:     TSMonitorSet(user.ts,Monitor,&user,NULL);
457:   }

459:   /* Set ODE initial conditions */
460:   VecGetArray(user.U,&x_ptr);
461:   x_ptr[0] = 2.0;
462:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
463:   VecRestoreArray(user.U,&x_ptr);

465:   /* Set runtime options */
466:   TSSetFromOptions(user.ts);

468:   /* Obtain the observation */
469:   TSSolve(user.ts,user.U);
470:   VecGetArray(user.U,&x_ptr);
471:   user.ob[0] = x_ptr[0];
472:   user.ob[1] = x_ptr[1];
473:   VecRestoreArray(user.U,&x_ptr);

475:   VecDuplicate(user.U,&x);
476:   VecGetArray(x,&x_ptr);
477:   x_ptr[0] = ic1;
478:   x_ptr[1] = ic2;
479:   VecRestoreArray(x,&x_ptr);

481:   /* Save trajectory of solution so that TSAdjointSolve() may be used */
482:   TSSetSaveTrajectory(user.ts);

484:   /* Compare finite difference and second-order adjoint. */
485:   switch (mode) {
486:     case 2 :
487:       FiniteDiff(x,arr,&user);
488:       PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n");
489:       PetscPrintf(PETSC_COMM_WORLD,"%g %g\n%g %g\n",(double)arr[0],(double)arr[1],(double)arr[2],(double)arr[3]);
490:       break;
491:     case 1 : /* Compute the Hessian column by column */
492:       VecCopy(x,user.U);
493:       VecGetArray(user.Dir,&x_ptr);
494:       x_ptr[0] = 1.;
495:       x_ptr[1] = 0.;
496:       VecRestoreArray(user.Dir,&x_ptr);
497:       Adjoint2(user.U,arr,&user);
498:       PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n");
499:       PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
500:       VecCopy(x,user.U);
501:       VecGetArray(user.Dir,&x_ptr);
502:       x_ptr[0] = 0.;
503:       x_ptr[1] = 1.;
504:       VecRestoreArray(user.Dir,&x_ptr);
505:       Adjoint2(user.U,arr,&user);
506:       PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n");
507:       PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
508:       break;
509:     case 0 :
510:       /* Create optimization context and set up */
511:       TaoCreate(PETSC_COMM_WORLD,&tao);
512:       TaoSetType(tao,TAOBLMVM);
513:       TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);

515:       if (mf) {
516:         MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H);
517:         MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat);
518:         MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
519:         TaoSetHessian(tao,user.H,user.H,MatrixFreeHessian,(void *)&user);
520:       } else { /* Create Hessian matrix */
521:         MatCreate(PETSC_COMM_WORLD,&user.H);
522:         MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2);
523:         MatSetUp(user.H);
524:         TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
525:       }

527:       /* Not use any preconditioner */
528:       TaoGetKSP(tao,&ksp);
529:       if (ksp) {
530:         KSPGetPC(ksp,&pc);
531:         PCSetType(pc,PCNONE);
532:       }

534:       /* Set initial solution guess */
535:       TaoSetSolution(tao,x);
536:       TaoSetFromOptions(tao);
537:       TaoSolve(tao);
538:       TaoDestroy(&tao);
539:       MatDestroy(&user.H);
540:       break;
541:     default:
542:       break;
543:   }

545:   /* Free work space.  All PETSc objects should be destroyed when they are no longer needed. */
546:   MatDestroy(&user.A);
547:   VecDestroy(&user.U);
548:   VecDestroy(&user.Lambda[0]);
549:   VecDestroy(&user.Lambda2[0]);
550:   VecDestroy(&user.Ihp1[0]);
551:   VecDestroy(&user.Dir);
552:   TSDestroy(&user.ts);
553:   VecDestroy(&x);
554:   PetscFinalize();
555:   return 0;
556: }

558: /*TEST
559:     build:
560:       requires: !complex !single

562:     test:
563:       args:  -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125
564:       output_file: output/ex20opt_ic_1.out

566:     test:
567:       suffix: 2
568:       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
569:       output_file: output/ex20opt_ic_2.out

571:     test:
572:       suffix: 3
573:       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
574:       output_file: output/ex20opt_ic_3.out

576:     test:
577:       suffix: 4
578:       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
579: TEST*/