Actual source code: eptorsion1.c

  1: /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------

  5:   Elastic-plastic torsion problem.

  7:   The elastic plastic torsion problem arises from the determination
  8:   of the stress field on an infinitely long cylindrical bar, which is
  9:   equivalent to the solution of the following problem:

 11:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}

 13:   where C is the torsion angle per unit length.

 15:   The multiprocessor version of this code is eptorsion2.c.

 17: ---------------------------------------------------------------------- */

 19: /*
 20:   Include "petsctao.h" so that we can use TAO solvers.  Note that this
 21:   file automatically includes files for lower-level support, such as those
 22:   provided by the PETSc library:
 23:      petsc.h       - base PETSc routines   petscvec.h - vectors
 24:      petscsys.h    - system routines        petscmat.h - matrices
 25:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 26:      petscviewer.h - viewers               petscpc.h  - preconditioners
 27: */

 29: #include <petsctao.h>

 31: static  char help[]=
 32: "Demonstrates use of the TAO package to solve \n\
 33: unconstrained minimization problems on a single processor.  This example \n\
 34: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
 35: test suite.\n\
 36: The command line options are:\n\
 37:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 38:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 39:   -par <param>, where <param> = angle of twist per unit length\n\n";

 41: /*T
 42:    Concepts: TAO^Solving an unconstrained minimization problem
 43:    Routines: TaoCreate(); TaoSetType();
 44:    Routines: TaoSetInitialVector();
 45:    Routines: TaoSetObjectiveAndGradientRoutine();
 46:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 47:    Routines: TaoGetKSP(); TaoSolve();
 48:    Routines: TaoDestroy();
 49:    Processors: 1
 50: T*/

 52: /*
 53:    User-defined application context - contains data needed by the
 54:    application-provided call-back routines, FormFunction(),
 55:    FormGradient(), and FormHessian().
 56: */

 58: typedef struct {
 59:    PetscReal  param;      /* nonlinearity parameter */
 60:    PetscInt   mx, my;     /* discretization in x- and y-directions */
 61:    PetscInt   ndim;       /* problem dimension */
 62:    Vec        s, y, xvec; /* work space for computing Hessian */
 63:    PetscReal  hx, hy;     /* mesh spacing in x- and y-directions */
 64: } AppCtx;

 66: /* -------- User-defined Routines --------- */

 68: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
 69: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
 70: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
 71: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
 72: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
 73: PetscErrorCode HessianProduct(void*,Vec,Vec);
 74: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
 75: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);

 77: PetscErrorCode main(int argc,char **argv)
 78: {
 79:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 80:   PetscInt           mx=10;               /* discretization in x-direction */
 81:   PetscInt           my=10;               /* discretization in y-direction */
 82:   Vec                x;                   /* solution, gradient vectors */
 83:   PetscBool          flg;                 /* A return value when checking for use options */
 84:   Tao                tao;                 /* Tao solver context */
 85:   Mat                H;                   /* Hessian matrix */
 86:   AppCtx             user;                /* application context */
 87:   PetscMPIInt        size;                /* number of processes */
 88:   PetscReal          one=1.0;

 90:   PetscBool          test_lmvm = PETSC_FALSE;
 91:   KSP                ksp;
 92:   PC                 pc;
 93:   Mat                M;
 94:   Vec                in, out, out2;
 95:   PetscReal          mult_solve_dist;

 97:   /* Initialize TAO,PETSc */
 98:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
 99:   MPI_Comm_size(MPI_COMM_WORLD,&size);
100:   if (size >1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors");

102:   /* Specify default parameters for the problem, check for command-line overrides */
103:   user.param = 5.0;
104:   PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
105:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
106:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
107:   PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg);

109:   PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
110:   PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",mx,my);
111:   user.ndim = mx * my; user.mx = mx; user.my = my;
112:   user.hx = one/(mx+1); user.hy = one/(my+1);

114:   /* Allocate vectors */
115:   VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
116:   VecDuplicate(user.y,&user.s);
117:   VecDuplicate(user.y,&x);

119:   /* Create TAO solver and set desired solution method */
120:   TaoCreate(PETSC_COMM_SELF,&tao);
121:   TaoSetType(tao,TAOLMVM);

123:   /* Set solution vector with an initial guess */
124:   FormInitialGuess(&user,x);
125:   TaoSetInitialVector(tao,x);

127:   /* Set routine for function and gradient evaluation */
128:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

130:   /* From command line options, determine if using matrix-free hessian */
131:   PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
132:   if (flg) {
133:     MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
134:     MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
135:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);

137:     TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);
138:   } else {
139:     MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
140:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
141:     TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
142:   }

144:   /* Test the LMVM matrix */
145:   if (test_lmvm) {
146:     PetscOptionsSetValue(NULL, "-tao_type", "bntr");
147:     PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
148:   }

150:   /* Check for any TAO command line options */
151:   TaoSetFromOptions(tao);

153:   /* SOLVE THE APPLICATION */
154:   TaoSolve(tao);

156:   /* Test the LMVM matrix */
157:   if (test_lmvm) {
158:     TaoGetKSP(tao, &ksp);
159:     KSPGetPC(ksp, &pc);
160:     PCLMVMGetMatLMVM(pc, &M);
161:     VecDuplicate(x, &in);
162:     VecDuplicate(x, &out);
163:     VecDuplicate(x, &out2);
164:     VecSet(in, 5.0);
165:     MatMult(M, in, out);
166:     MatSolve(M, out, out2);
167:     VecAXPY(out2, -1.0, in);
168:     VecNorm(out2, NORM_2, &mult_solve_dist);
169:     PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist);
170:     VecDestroy(&in);
171:     VecDestroy(&out);
172:     VecDestroy(&out2);
173:   }

175:   TaoDestroy(&tao);
176:   VecDestroy(&user.s);
177:   VecDestroy(&user.y);
178:   VecDestroy(&x);
179:   MatDestroy(&H);

181:   PetscFinalize();
182:   return ierr;
183: }

185: /* ------------------------------------------------------------------- */
186: /*
187:     FormInitialGuess - Computes an initial approximation to the solution.

189:     Input Parameters:
190: .   user - user-defined application context
191: .   X    - vector

193:     Output Parameters:
194: .   X    - vector
195: */
196: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
197: {
198:   PetscReal      hx = user->hx, hy = user->hy, temp;
199:   PetscReal      val;
201:   PetscInt       i, j, k, nx = user->mx, ny = user->my;

203:   /* Compute initial guess */
205:   for (j=0; j<ny; j++) {
206:     temp = PetscMin(j+1,ny-j)*hy;
207:     for (i=0; i<nx; i++) {
208:       k   = nx*j + i;
209:       val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
210:       VecSetValues(X,1,&k,&val,ADD_VALUES);
211:     }
212:   }
213:   VecAssemblyBegin(X);
214:   VecAssemblyEnd(X);
215:   return(0);
216: }

218: /* ------------------------------------------------------------------- */
219: /*
220:    FormFunctionGradient - Evaluates the function and corresponding gradient.

222:    Input Parameters:
223:    tao - the Tao context
224:    X   - the input vector
225:    ptr - optional user-defined context, as set by TaoSetFunction()

227:    Output Parameters:
228:    f   - the newly evaluated function
229:    G   - the newly evaluated gradient
230: */
231: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
232: {

236:   FormFunction(tao,X,f,ptr);
237:   FormGradient(tao,X,G,ptr);
238:   return(0);
239: }

241: /* ------------------------------------------------------------------- */
242: /*
243:    FormFunction - Evaluates the function, f(X).

245:    Input Parameters:
246: .  tao - the Tao context
247: .  X   - the input vector
248: .  ptr - optional user-defined context, as set by TaoSetFunction()

250:    Output Parameters:
251: .  f    - the newly evaluated function
252: */
253: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
254: {
255:   AppCtx            *user = (AppCtx *) ptr;
256:   PetscReal         hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
257:   PetscReal         zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
258:   PetscReal         v, cdiv3 = user->param/three;
259:   const PetscScalar *x;
260:   PetscErrorCode    ierr;
261:   PetscInt          nx = user->mx, ny = user->my, i, j, k;

264:   /* Get pointer to vector data */
265:   VecGetArrayRead(X,&x);

267:   /* Compute function contributions over the lower triangular elements */
268:   for (j=-1; j<ny; j++) {
269:     for (i=-1; i<nx; i++) {
270:       k = nx*j + i;
271:       v = zero;
272:       vr = zero;
273:       vt = zero;
274:       if (i >= 0 && j >= 0) v = x[k];
275:       if (i < nx-1 && j > -1) vr = x[k+1];
276:       if (i > -1 && j < ny-1) vt = x[k+nx];
277:       dvdx = (vr-v)/hx;
278:       dvdy = (vt-v)/hy;
279:       fquad += dvdx*dvdx + dvdy*dvdy;
280:       flin -= cdiv3*(v+vr+vt);
281:     }
282:   }

284:   /* Compute function contributions over the upper triangular elements */
285:   for (j=0; j<=ny; j++) {
286:     for (i=0; i<=nx; i++) {
287:       k = nx*j + i;
288:       vb = zero;
289:       vl = zero;
290:       v = zero;
291:       if (i < nx && j > 0) vb = x[k-nx];
292:       if (i > 0 && j < ny) vl = x[k-1];
293:       if (i < nx && j < ny) v = x[k];
294:       dvdx = (v-vl)/hx;
295:       dvdy = (v-vb)/hy;
296:       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
297:       flin = flin - cdiv3*(vb+vl+v);
298:     }
299:   }

301:   /* Restore vector */
302:   VecRestoreArrayRead(X,&x);

304:   /* Scale the function */
305:   area = p5*hx*hy;
306:   *f = area*(p5*fquad+flin);

308:   PetscLogFlops(24.0*nx*ny);
309:   return(0);
310: }

312: /* ------------------------------------------------------------------- */
313: /*
314:     FormGradient - Evaluates the gradient, G(X).

316:     Input Parameters:
317: .   tao  - the Tao context
318: .   X    - input vector
319: .   ptr  - optional user-defined context

321:     Output Parameters:
322: .   G - vector containing the newly evaluated gradient
323: */
324: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
325: {
326:   AppCtx            *user = (AppCtx *) ptr;
327:   PetscReal         zero=0.0, p5=0.5, three = 3.0, area, val;
328:   PetscErrorCode    ierr;
329:   PetscInt          nx = user->mx, ny = user->my, ind, i, j, k;
330:   PetscReal         hx = user->hx, hy = user->hy;
331:   PetscReal         vb, vl, vr, vt, dvdx, dvdy;
332:   PetscReal         v, cdiv3 = user->param/three;
333:   const PetscScalar *x;

336:   /* Initialize gradient to zero */
337:   VecSet(G, zero);

339:   /* Get pointer to vector data */
340:   VecGetArrayRead(X,&x);

342:   /* Compute gradient contributions over the lower triangular elements */
343:   for (j=-1; j<ny; j++) {
344:     for (i=-1; i<nx; i++) {
345:       k  = nx*j + i;
346:       v  = zero;
347:       vr = zero;
348:       vt = zero;
349:       if (i >= 0 && j >= 0)    v = x[k];
350:       if (i < nx-1 && j > -1) vr = x[k+1];
351:       if (i > -1 && j < ny-1) vt = x[k+nx];
352:       dvdx = (vr-v)/hx;
353:       dvdy = (vt-v)/hy;
354:       if (i != -1 && j != -1) {
355:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
356:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
357:       }
358:       if (i != nx-1 && j != -1) {
359:         ind = k+1; val =  dvdx/hx - cdiv3;
360:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
361:       }
362:       if (i != -1 && j != ny-1) {
363:         ind = k+nx; val = dvdy/hy - cdiv3;
364:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
365:       }
366:     }
367:   }

369:   /* Compute gradient contributions over the upper triangular elements */
370:   for (j=0; j<=ny; j++) {
371:     for (i=0; i<=nx; i++) {
372:       k = nx*j + i;
373:       vb = zero;
374:       vl = zero;
375:       v  = zero;
376:       if (i < nx && j > 0) vb = x[k-nx];
377:       if (i > 0 && j < ny) vl = x[k-1];
378:       if (i < nx && j < ny) v = x[k];
379:       dvdx = (v-vl)/hx;
380:       dvdy = (v-vb)/hy;
381:       if (i != nx && j != 0) {
382:         ind = k-nx; val = - dvdy/hy - cdiv3;
383:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
384:       }
385:       if (i != 0 && j != ny) {
386:         ind = k-1; val =  - dvdx/hx - cdiv3;
387:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
388:       }
389:       if (i != nx && j != ny) {
390:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
391:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
392:       }
393:     }
394:   }
395:   VecRestoreArrayRead(X,&x);

397:   /* Assemble gradient vector */
398:   VecAssemblyBegin(G);
399:   VecAssemblyEnd(G);

401:   /* Scale the gradient */
402:   area = p5*hx*hy;
403:   VecScale(G, area);
404:   PetscLogFlops(24.0*nx*ny);
405:   return(0);
406: }

408: /* ------------------------------------------------------------------- */
409: /*
410:    FormHessian - Forms the Hessian matrix.

412:    Input Parameters:
413: .  tao - the Tao context
414: .  X    - the input vector
415: .  ptr  - optional user-defined context, as set by TaoSetHessian()

417:    Output Parameters:
418: .  H     - Hessian matrix
419: .  PrecH - optionally different preconditioning Hessian
420: .  flag  - flag indicating matrix structure

422:    Notes:
423:    This routine is intended simply as an example of the interface
424:    to a Hessian evaluation routine.  Since this example compute the
425:    Hessian a column at a time, it is not particularly efficient and
426:    is not recommended.
427: */
428: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
429: {
430:   AppCtx         *user = (AppCtx *) ptr;
432:   PetscInt       i,j, ndim = user->ndim;
433:   PetscReal      *y, zero = 0.0, one = 1.0;
434:   PetscBool      assembled;

437:   user->xvec = X;

439:   /* Initialize Hessian entries and work vector to zero */
440:   MatAssembled(H,&assembled);
441:   if (assembled){MatZeroEntries(H);}

443:   VecSet(user->s, zero);

445:   /* Loop over matrix columns to compute entries of the Hessian */
446:   for (j=0; j<ndim; j++) {
447:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
448:     VecAssemblyBegin(user->s);
449:     VecAssemblyEnd(user->s);

451:     HessianProduct(ptr,user->s,user->y);

453:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
454:     VecAssemblyBegin(user->s);
455:     VecAssemblyEnd(user->s);

457:     VecGetArray(user->y,&y);
458:     for (i=0; i<ndim; i++) {
459:       if (y[i] != zero) {
460:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
461:       }
462:     }
463:     VecRestoreArray(user->y,&y);
464:   }
465:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
466:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
467:   return(0);
468: }

470: /* ------------------------------------------------------------------- */
471: /*
472:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
473:    products.

475:    Input Parameters:
476: .  tao - the Tao context
477: .  X    - the input vector
478: .  ptr  - optional user-defined context, as set by TaoSetHessian()

480:    Output Parameters:
481: .  H     - Hessian matrix
482: .  PrecH - optionally different preconditioning Hessian
483: .  flag  - flag indicating matrix structure
484: */
485: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
486: {
487:   AppCtx *user = (AppCtx *) ptr;

489:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
491:   user->xvec = X;
492:   return(0);
493: }

495: /* ------------------------------------------------------------------- */
496: /*
497:    HessianProductMat - Computes the matrix-vector product
498:    y = mat*svec.

500:    Input Parameters:
501: .  mat  - input matrix
502: .  svec - input vector

504:    Output Parameters:
505: .  y    - solution vector
506: */
507: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
508: {
509:   void           *ptr;

513:   MatShellGetContext(mat,&ptr);
514:   HessianProduct(ptr,svec,y);
515:   return(0);
516: }

518: /* ------------------------------------------------------------------- */
519: /*
520:    Hessian Product - Computes the matrix-vector product:
521:    y = f''(x)*svec.

523:    Input Parameters:
524: .  ptr  - pointer to the user-defined context
525: .  svec - input vector

527:    Output Parameters:
528: .  y    - product vector
529: */
530: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
531: {
532:   AppCtx            *user = (AppCtx *)ptr;
533:   PetscReal         p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
534:   const PetscScalar *x, *s;
535:   PetscReal         v, vb, vl, vr, vt, hxhx, hyhy;
536:   PetscErrorCode    ierr;
537:   PetscInt          nx, ny, i, j, k, ind;

540:   nx   = user->mx;
541:   ny   = user->my;
542:   hx   = user->hx;
543:   hy   = user->hy;
544:   hxhx = one/(hx*hx);
545:   hyhy = one/(hy*hy);

547:   /* Get pointers to vector data */
548:   VecGetArrayRead(user->xvec,&x);
549:   VecGetArrayRead(svec,&s);

551:   /* Initialize product vector to zero */
552:   VecSet(y, zero);

554:   /* Compute f''(x)*s over the lower triangular elements */
555:   for (j=-1; j<ny; j++) {
556:     for (i=-1; i<nx; i++) {
557:        k = nx*j + i;
558:        v = zero;
559:        vr = zero;
560:        vt = zero;
561:        if (i != -1 && j != -1) v = s[k];
562:        if (i != nx-1 && j != -1) {
563:          vr = s[k+1];
564:          ind = k+1; val = hxhx*(vr-v);
565:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
566:        }
567:        if (i != -1 && j != ny-1) {
568:          vt = s[k+nx];
569:          ind = k+nx; val = hyhy*(vt-v);
570:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
571:        }
572:        if (i != -1 && j != -1) {
573:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
574:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
575:        }
576:      }
577:    }

579:   /* Compute f''(x)*s over the upper triangular elements */
580:   for (j=0; j<=ny; j++) {
581:     for (i=0; i<=nx; i++) {
582:        k = nx*j + i;
583:        v = zero;
584:        vl = zero;
585:        vb = zero;
586:        if (i != nx && j != ny) v = s[k];
587:        if (i != nx && j != 0) {
588:          vb = s[k-nx];
589:          ind = k-nx; val = hyhy*(vb-v);
590:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
591:        }
592:        if (i != 0 && j != ny) {
593:          vl = s[k-1];
594:          ind = k-1; val = hxhx*(vl-v);
595:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
596:        }
597:        if (i != nx && j != ny) {
598:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
599:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
600:        }
601:     }
602:   }
603:   /* Restore vector data */
604:   VecRestoreArrayRead(svec,&s);
605:   VecRestoreArrayRead(user->xvec,&x);

607:   /* Assemble vector */
608:   VecAssemblyBegin(y);
609:   VecAssemblyEnd(y);

611:   /* Scale resulting vector by area */
612:   area = p5*hx*hy;
613:   VecScale(y, area);
614:   PetscLogFlops(18.0*nx*ny);
615:   return(0);
616: }

618: /*TEST

620:    build:
621:       requires: !complex

623:    test:
624:       suffix: 1
625:       args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4

627:    test:
628:       suffix: 2
629:       args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4

631:    test:
632:       suffix: 3
633:       args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian

635:    test:
636:      suffix: 4
637:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls

639:    test:
640:      suffix: 5
641:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm

643:    test:
644:      suffix: 6
645:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1

647: TEST*/