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++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
 99:   }
100:   VecRestoreArrayRead(U, &u);
101:   VecRestoreArrayRead(Vl[0], &vl);
102:   VecRestoreArrayRead(Vr, &vr);
103:   VecRestoreArray(VHV[0], &vhv);
104:   return 0;
105: }

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

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

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

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

135:   VecGetArrayRead(U, &u);
136:   J[0][0] = a;
137:   J[0][1] = -1.0;
138:   J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
139:   J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
140:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
141:   VecRestoreArrayRead(U, &u);
142:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
143:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
144:   if (A != B) {
145:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
146:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
147:   }
148:   return 0;
149: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

276:   VecCopy(U, user_ptr->U);
277:   return 0;
278: }

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

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

292:   TSAdjointReset(ts);

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

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

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

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

330:   TSAdjointReset(ts);
331:   TSAdjointResetForward(ts);
332:   return 0;
333: }

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

344:   VecDuplicate(U, &Up);
345:   VecDuplicate(U, &G);
346:   VecDuplicate(U, &Gp);

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

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

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

374:   VecDestroy(&G);
375:   VecDestroy(&Gp);
376:   VecDestroy(&Up);
377:   return 0;
378: }

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

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

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

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

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

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

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

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

454:   if (monitor) TSMonitorSet(user.ts, Monitor, &user, NULL);

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

462:   /* Set runtime options */
463:   TSSetFromOptions(user.ts);

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

472:   VecDuplicate(user.U, &x);
473:   VecGetArray(x, &x_ptr);
474:   x_ptr[0] = ic1;
475:   x_ptr[1] = ic2;
476:   VecRestoreArray(x, &x_ptr);

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

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

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

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

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

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

555: /*TEST
556:     build:
557:       requires: !complex !single

559:     test:
560:       args:  -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125
561:       output_file: output/ex20opt_ic_1.out

563:     test:
564:       suffix: 2
565:       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
566:       output_file: output/ex20opt_ic_2.out

568:     test:
569:       suffix: 3
570:       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
571:       output_file: output/ex20opt_ic_3.out

573:     test:
574:       suffix: 4
575:       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
576: TEST*/