Actual source code: ex20opt_ic.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  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*/