Actual source code: ex20opt_ic.c

petsc-master 2019-05-20
Report Typos and Errors
  1: #define c11 1.0
  2: #define c12 0
  3: #define c21 2.0
  4: #define c22 1.0
  5: static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n";

  7: /**
  8:   Concepts: TS^time-dependent nonlinear problems
  9:   Concepts: TS^van der Pol equation DAE equivalent
 10:   Concepts: Optimization using adjoint sensitivities
 11:   Processors: 1
 12: */
 13: /**
 14:   Notes:
 15:   This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS.
 16:   The nonlinear problem is written in an ODE equivalent form.
 17:   The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions.
 18:   The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details.
 19: */

 21:  #include <petsctao.h>
 22:  #include <petscts.h>

 24: typedef struct _n_User *User;
 25: struct _n_User {
 26:   TS        ts;
 27:   PetscReal mu;
 28:   PetscReal next_output;

 30:   /* Sensitivity analysis support */
 31:   PetscInt  steps;
 32:   PetscReal ftime;
 33:   Mat       A;                       /* Jacobian matrix for ODE */
 34:   Mat       Jacp;                    /* JacobianP matrix for ODE*/
 35:   Mat       H;                       /* Hessian matrix for optimization */
 36:   Vec       U,Lambda[1],Mup[1];      /* first-order adjoint variables */
 37:   Vec       Lambda2[2];              /* second-order adjoint variables */
 38:   Vec       Ihp1[1];                 /* working space for Hessian evaluations */
 39:   Vec       Dir;                     /* direction vector */
 40:   PetscReal ob[2];                   /* observation used by the cost function */
 41:   PetscBool implicitform;            /* implicit ODE? */
 42: };
 43: PetscErrorCode Adjoint2(Vec,PetscScalar[],User);

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

 47: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
 48: {
 49:   PetscErrorCode    ierr;
 50:   User              user = (User)ctx;
 51:   PetscScalar       *f;
 52:   const PetscScalar *u;

 55:   VecGetArrayRead(U,&u);
 56:   VecGetArray(F,&f);
 57:   f[0] = u[1];
 58:   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
 59:   VecRestoreArrayRead(U,&u);
 60:   VecRestoreArray(F,&f);
 61:   return(0);
 62: }

 64: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
 65: {
 66:   PetscErrorCode    ierr;
 67:   User              user = (User)ctx;
 68:   PetscReal         mu   = user->mu;
 69:   PetscInt          rowcol[] = {0,1};
 70:   PetscScalar       J[2][2];
 71:   const PetscScalar *u;

 74:   VecGetArrayRead(U,&u);
 75:   J[0][0] = 0;
 76:   J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
 77:   J[0][1] = 1.0;
 78:   J[1][1] = mu*(1.0-u[0]*u[0]);
 79:   MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
 80:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 82:   if (A != B) {
 83:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 84:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 85:   }
 86:   VecRestoreArrayRead(U,&u);
 87:   return(0);
 88: }

 90: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
 91: {
 92:   const PetscScalar *vl,*vr,*u;
 93:   PetscScalar       *vhv;
 94:   PetscScalar       dJdU[2][2][2]={{{0}}};
 95:   PetscInt          i,j,k;
 96:   User              user = (User)ctx;
 97:   PetscErrorCode    ierr;

100:   VecGetArrayRead(U,&u);
101:   VecGetArrayRead(Vl[0],&vl);
102:   VecGetArrayRead(Vr,&vr);
103:   VecGetArray(VHV[0],&vhv);

105:   dJdU[1][0][0] = -2.*user->mu*u[1];
106:   dJdU[1][1][0] = -2.*user->mu*u[0];
107:   dJdU[1][0][1] = -2.*user->mu*u[0];
108:   for (j=0;j<2;j++) {
109:     vhv[j] = 0;
110:     for (k=0;k<2;k++)
111:       for (i=0;i<2;i++ )
112:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
113:   }
114:   VecRestoreArrayRead(U,&u);
115:   VecRestoreArrayRead(Vl[0],&vl);
116:   VecRestoreArrayRead(Vr,&vr);
117:   VecRestoreArray(VHV[0],&vhv);
118:   return(0);
119: }

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

123: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
124: {
125:   PetscErrorCode    ierr;
126:   User              user = (User)ctx;
127:   const PetscScalar *u,*udot;
128:   PetscScalar       *f;

131:   VecGetArrayRead(U,&u);
132:   VecGetArrayRead(Udot,&udot);
133:   VecGetArray(F,&f);
134:   f[0] = udot[0] - u[1];
135:   f[1] = c21*(udot[0]-u[1]) + udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
136:   VecRestoreArrayRead(U,&u);
137:   VecRestoreArrayRead(Udot,&udot);
138:   VecRestoreArray(F,&f);
139:   return(0);
140: }

142: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
143: {
144:   PetscErrorCode    ierr;
145:   User              user = (User)ctx;
146:   PetscInt          rowcol[] = {0,1};
147:   PetscScalar       J[2][2];
148:   const PetscScalar *u;

151:   VecGetArrayRead(U,&u);
152:   J[0][0] = a;     J[0][1] =  -1.0;
153:   J[1][0] = c21*a + user->mu*(1.0 + 2.0*u[0]*u[1]);   J[1][1] = -c21 + a - user->mu*(1.0-u[0]*u[0]);
154:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
155:   VecRestoreArrayRead(U,&u);
156:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
157:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
158:   if (A != B) {
159:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
160:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
161:   }
162:   return(0);
163: }

165: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
166: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
167: {
168:   PetscErrorCode    ierr;
169:   const PetscScalar *u;
170:   PetscReal         tfinal, dt;
171:   User              user = (User)ctx;
172:   Vec               interpolatedU;

175:   TSGetTimeStep(ts,&dt);
176:   TSGetMaxTime(ts,&tfinal);

178:   while (user->next_output <= t && user->next_output <= tfinal) {
179:     VecDuplicate(U,&interpolatedU);
180:     TSInterpolate(ts,user->next_output,interpolatedU);
181:     VecGetArrayRead(interpolatedU,&u);
182:     PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
183:                        (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
184:                        (double)PetscRealPart(u[1]));
185:     VecRestoreArrayRead(interpolatedU,&u);
186:     VecDestroy(&interpolatedU);
187:     user->next_output += 0.1;
188:   }
189:   return(0);
190: }

192: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
193: {
194:   const PetscScalar *vl,*vr,*u;
195:   PetscScalar       *vhv;
196:   PetscScalar       dJdU[2][2][2]={{{0}}};
197:   PetscInt          i,j,k;
198:   User              user = (User)ctx;
199:   PetscErrorCode    ierr;

202:   VecGetArrayRead(U,&u);
203:   VecGetArrayRead(Vl[0],&vl);
204:   VecGetArrayRead(Vr,&vr);
205:   VecGetArray(VHV[0],&vhv);
206:   dJdU[1][0][0] = 2.*user->mu*u[1];
207:   dJdU[1][1][0] = 2.*user->mu*u[0];
208:   dJdU[1][0][1] = 2.*user->mu*u[0];
209:   for (j=0;j<2;j++) {
210:     vhv[j] = 0;
211:     for (k=0;k<2;k++)
212:       for (i=0;i<2;i++ )
213:         vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
214:   }
215:   VecRestoreArrayRead(U,&u);
216:   VecRestoreArrayRead(Vl[0],&vl);
217:   VecRestoreArrayRead(Vr,&vr);
218:   VecRestoreArray(VHV[0],&vhv);
219:   return(0);
220: }

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

224: static PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx)
225: {
226:   User              user_ptr = (User)ctx;
227:   TS                ts = user_ptr->ts;
228:   const PetscScalar *x_ptr;
229:   PetscScalar       *y_ptr;

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

235:   TSSetTime(ts,0.0);
236:   TSSetStepNumber(ts,0);
237:   TSSetTimeStep(ts,0.001); /* can be overwritten by command line options */
238:   TSSetFromOptions(ts);
239:   TSSolve(ts,user_ptr->U);

241:   VecGetArrayRead(user_ptr->U,&x_ptr);
242:   VecGetArray(user_ptr->Lambda[0],&y_ptr);
243:   *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]);;
244:   y_ptr[0] = 2.*(x_ptr[0]-user_ptr->ob[0]);
245:   y_ptr[1] = 2.*(x_ptr[1]-user_ptr->ob[1]);
246:   VecRestoreArray(user_ptr->Lambda[0],&y_ptr);
247:   VecRestoreArrayRead(user_ptr->U,&x_ptr);

249:   TSSetCostGradients(ts,1,user_ptr->Lambda,NULL);
250:   TSAdjointSolve(ts);
251:   VecCopy(user_ptr->Lambda[0],G);
252:   return(0);
253: }

255: static PetscErrorCode FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
256: {
257:   User           user_ptr = (User)ctx;
258:   PetscScalar    harr[2];
259:   PetscScalar    *x_ptr;
260:   const PetscInt rows[2] = {0,1};
261:   PetscInt       col;

265:   VecCopy(U,user_ptr->U);
266:   VecGetArray(user_ptr->Dir,&x_ptr);
267:   x_ptr[0] = 1.;
268:   x_ptr[1] = 0.;
269:   VecRestoreArray(user_ptr->Dir,&x_ptr);
270:   Adjoint2(user_ptr->U,harr,user_ptr);
271:   col      = 0;
272:   MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);

274:   VecCopy(U,user_ptr->U);
275:   VecGetArray(user_ptr->Dir,&x_ptr);
276:   x_ptr[0] = 0.;
277:   x_ptr[1] = 1.;
278:   VecRestoreArray(user_ptr->Dir,&x_ptr);
279:   Adjoint2(user_ptr->U,harr,user_ptr);
280:   col      = 1;
281:   MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES);

283:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
284:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
285:   if (H != Hpre) {
286:     MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
287:     MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
288:   }
289:   return(0);
290: }

292: static PetscErrorCode MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx)
293: {
294:   User           user_ptr = (User)ctx;

298:   VecCopy(U,user_ptr->U);
299:   return(0);
300: }

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

304: /**
305:   Compute the Hessian-vector product for the cost function using Second-order adjoint
306: */
307: PetscErrorCode Adjoint2(Vec U,PetscScalar arr[],User ctx)
308: {
309:   TS             ts = ctx->ts;
310:   PetscScalar    *x_ptr,*y_ptr;
311:   Mat            tlmsen;

315:   TSAdjointReset(ts);

317:   TSSetTime(ts,0.0);
318:   TSSetStepNumber(ts,0);
319:   TSSetTimeStep(ts,0.001);
320:   TSSetFromOptions(ts);
321:   TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir);
322:   TSAdjointSetForward(ts,NULL);
323:   TSSolve(ts,U);

325:   /* Set terminal conditions for first- and second-order adjonts */
326:   VecGetArray(U,&x_ptr);
327:   VecGetArray(ctx->Lambda[0],&y_ptr);
328:   y_ptr[0] = 2.*(x_ptr[0]-ctx->ob[0]);
329:   y_ptr[1] = 2.*(x_ptr[1]-ctx->ob[1]);
330:   VecRestoreArray(ctx->Lambda[0],&y_ptr);
331:   VecRestoreArray(U,&x_ptr);
332:   TSForwardGetSensitivities(ts,NULL,&tlmsen);
333:   MatDenseGetColumn(tlmsen,0,&x_ptr);
334:   VecGetArray(ctx->Lambda2[0],&y_ptr);
335:   y_ptr[0] = 2.*x_ptr[0];
336:   y_ptr[1] = 2.*x_ptr[1];
337:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
338:   MatDenseRestoreColumn(tlmsen,&x_ptr);

340:   TSSetCostGradients(ts,1,ctx->Lambda,NULL);
341:   if (ctx->implicitform) {
342:     TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
343:   } else {
344:     TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx);
345:   }
346:   TSAdjointSolve(ts);

348:   VecGetArray(ctx->Lambda2[0],&x_ptr);
349:   arr[0] = x_ptr[0];
350:   arr[1] = x_ptr[1];
351:   VecRestoreArray(ctx->Lambda2[0],&x_ptr);

353:   TSAdjointReset(ts);
354:   TSAdjointResetForward(ts);
355:   return(0);
356: }

358: PetscErrorCode FiniteDiff(Vec U,PetscScalar arr[],User ctx)
359: {
360:   Vec               Up,G,Gp;
361:   const PetscScalar eps = PetscRealConstant(1e-7);
362:   PetscScalar       *u;
363:   Tao               tao = NULL;
364:   PetscReal         f;
365:   PetscErrorCode    ierr;

368:   VecDuplicate(U,&Up);
369:   VecDuplicate(U,&G);
370:   VecDuplicate(U,&Gp);

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

374:   VecCopy(U,Up);
375:   VecGetArray(Up,&u);
376:   u[0] += eps;
377:   VecRestoreArray(Up,&u);
378:   FormFunctionGradient(tao,Up,&f,Gp,ctx);
379:   VecAXPY(Gp,-1,G);
380:   VecScale(Gp,1./eps);
381:   VecGetArray(Gp,&u);
382:   arr[0] = u[0];
383:   arr[1] = u[1];
384:   VecRestoreArray(Gp,&u);

386:   VecCopy(U,Up);
387:   VecGetArray(Up,&u);
388:   u[1] += eps;
389:   VecRestoreArray(Up,&u);
390:   FormFunctionGradient(tao,Up,&f,Gp,ctx);
391:   VecAXPY(Gp,-1,G);
392:   VecScale(Gp,1./eps);
393:   VecGetArray(Gp,&u);
394:   arr[2] = u[0];
395:   arr[3] = u[1];
396:   VecRestoreArray(Gp,&u);

398:   VecDestroy(&G);
399:   VecDestroy(&Gp);
400:   VecDestroy(&Up);
401:   return(0);
402: }

404: static PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
405: {
406:   User           user_ptr;
407:   PetscScalar    *y_ptr;

411:   MatShellGetContext(mat,(void*)&user_ptr);
412:   VecCopy(svec,user_ptr->Dir);
413:   VecGetArray(y,&y_ptr);
414:   Adjoint2(user_ptr->U,y_ptr,user_ptr);
415:   VecRestoreArray(y,&y_ptr);
416:   return(0);
417: }

419: int main(int argc,char **argv)
420: {
421:   PetscBool      monitor = PETSC_FALSE,mf = PETSC_TRUE;
422:   PetscInt       mode = 0;
423:   PetscMPIInt    size;
424:   struct _n_User user;
425:   Vec            x; /* working vector for TAO */
426:   PetscScalar    *x_ptr,arr[4];
427:   PetscScalar    ic1 = 2.2,ic2 = -0.7; /* initial guess for TAO */
428:   Tao            tao;
429:   KSP            ksp;
430:   PC             pc;

433:   /* Initialize program */
434:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
435:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
436:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

438:   /* Set runtime options */
439:   user.next_output  = 0.0;
440:   user.mu           = 1.0e3;
441:   user.steps        = 0;
442:   user.ftime        = 0.5;
443:   user.implicitform = PETSC_TRUE;

445:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
446:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
447:   PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL);
448:   PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL);
449:   PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL);
450:   PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL); /* matrix-free hessian for optimization */
451:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);

453:   /* Create necessary matrix and vectors, solve same ODE on every process */
454:   MatCreate(PETSC_COMM_WORLD,&user.A);
455:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
456:   MatSetFromOptions(user.A);
457:   MatSetUp(user.A);
458:   MatCreateVecs(user.A,&user.U,NULL);
459:   MatCreateVecs(user.A,&user.Dir,NULL);
460:   MatCreateVecs(user.A,&user.Lambda[0],NULL);
461:   MatCreateVecs(user.A,&user.Lambda2[0],NULL);
462:   MatCreateVecs(user.A,&user.Ihp1[0],NULL);

464:   /* Create timestepping solver context */
465:   TSCreate(PETSC_COMM_WORLD,&user.ts);
466:   if (user.implicitform) {
467:     TSSetIFunction(user.ts,NULL,IFunction,&user);
468:     TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
469:     TSSetType(user.ts,TSCN);
470:   } else {
471:     TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
472:     TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
473:     TSSetType(user.ts,TSRK);
474:   }
475:   TSSetMaxTime(user.ts,user.ftime);
476:   TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);

478:   if (monitor) {
479:     TSMonitorSet(user.ts,Monitor,&user,NULL);
480:   }

482:   /* Set ODE initial conditions */
483:   VecGetArray(user.U,&x_ptr);
484:   x_ptr[0] = 2.0;
485:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
486:   VecRestoreArray(user.U,&x_ptr);

488:   /* Set runtime options */
489:   TSSetFromOptions(user.ts);

491:   /* Obtain the observation */
492:   TSSolve(user.ts,user.U);
493:   VecGetArray(user.U,&x_ptr);
494:   user.ob[0] = x_ptr[0];
495:   user.ob[1] = x_ptr[1];
496:   VecRestoreArray(user.U,&x_ptr);

498:   VecDuplicate(user.U,&x);
499:   VecGetArray(x,&x_ptr);
500:   x_ptr[0] = ic1;
501:   x_ptr[1] = ic2;
502:   VecRestoreArray(x,&x_ptr);

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

507:   /* Compare finite difference and second-order adjoint. */
508:   switch (mode) {
509:     case 2 :
510:       FiniteDiff(x,arr,&user);
511:       PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n");
512:       PetscPrintf(PETSC_COMM_WORLD,"%g %g\n%g %g\n",(double)arr[0],(double)arr[1],(double)arr[2],(double)arr[3]);
513:       break;
514:     case 1 : /* Compute the Hessian column by column */
515:       VecCopy(x,user.U);
516:       VecGetArray(user.Dir,&x_ptr);
517:       x_ptr[0] = 1.;
518:       x_ptr[1] = 0.;
519:       VecRestoreArray(user.Dir,&x_ptr);
520:       Adjoint2(user.U,arr,&user);
521:       PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n");
522:       PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
523:       VecCopy(x,user.U);
524:       VecGetArray(user.Dir,&x_ptr);
525:       x_ptr[0] = 0.;
526:       x_ptr[1] = 1.;
527:       VecRestoreArray(user.Dir,&x_ptr);
528:       Adjoint2(user.U,arr,&user);
529:       PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n");
530:       PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1]);
531:       break;
532:     case 0 :
533:       /* Create optimization context and set up */
534:       TaoCreate(PETSC_COMM_WORLD,&tao);
535:       TaoSetType(tao,TAOBLMVM);
536:       TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

538:       if (mf) {
539:         MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H);
540:         MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat);
541:         MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
542:         TaoSetHessianRoutine(tao,user.H,user.H,MatrixFreeHessian,(void *)&user);
543:       } else { /* Create Hessian matrix */
544:         MatCreate(PETSC_COMM_WORLD,&user.H);
545:         MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2);
546:         MatSetUp(user.H);
547:         TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
548:       }

550:       /* Not use any preconditioner */
551:       TaoGetKSP(tao,&ksp);
552:       if (ksp) {
553:         KSPGetPC(ksp,&pc);
554:         PCSetType(pc,PCNONE);
555:       }

557:       /* Set initial solution guess */
558:       TaoSetInitialVector(tao,x);
559:       TaoSetFromOptions(tao);
560:       TaoSolve(tao);
561:       TaoDestroy(&tao);
562:       MatDestroy(&user.H);
563:       break;
564:     default:
565:       break;
566:   }

568:   /* Free work space.  All PETSc objects should be destroyed when they are no longer needed. */
569:   MatDestroy(&user.A);
570:   VecDestroy(&user.U);
571:   VecDestroy(&user.Lambda[0]);
572:   VecDestroy(&user.Lambda2[0]);
573:   VecDestroy(&user.Ihp1[0]);
574:   VecDestroy(&user.Dir);
575:   TSDestroy(&user.ts);
576:   VecDestroy(&x);
577:   PetscFinalize();
578:   return(ierr);
579: }

581: /*TEST
582:     build:
583:       requires: !complex !single

585:     test:
586:       args:  -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125
587:       output_file: output/ex20opt_ic_1.out

589:     test:
590:       suffix: 2
591:       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
592:       output_file: output/ex20opt_ic_2.out

594:     test:
595:       suffix: 3
596:       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
597:       output_file: output/ex20opt_ic_3.out

599:     test:
600:       suffix: 4
601:       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
602: TEST*/