Actual source code: eptorsion1.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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    - sysem 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>


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

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

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

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

 67: /* -------- User-defined Routines --------- */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

309:   PetscLogFlops(nx*ny*24);
310:   return(0);
311: }

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

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

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

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

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

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

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

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

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

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

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

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

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

438:   user->xvec = X;

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

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

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

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

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

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

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

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

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

490:   /* 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(nx*ny*18);
615:   return(0);
616: }


619: /*TEST

621:    build:
622:       requires: !complex

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

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

632:    test:
633:       suffix: 3
634:       args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
635:       
636:    test:
637:      suffix: 4
638:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
639:      
640:    test:
641:      suffix: 5
642:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
643:      
644:    test:
645:      suffix: 6
646:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1

648: TEST*/