Actual source code: eptorsion1.c

      1: /*$Id: eptorsion1.c 1.49 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/

      3: /* Program usage: mpirun -np 1 eptorsion1 [-help] [all TAO options] */

      5: /* ----------------------------------------------------------------------

      7:   Elastic-plastic torsion problem.

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

     13:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
     14:
     15:   where C is the torsion angle per unit length.

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

     19: ---------------------------------------------------------------------- */

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

     31:  #include tao.h

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

     44: /*T
     45:    Concepts: TAO - Solving an unconstrained minimization problem
     46:    Routines: TaoInitialize(); TaoFinalize();
     47:    Routines: TaoApplicationCreate(); TaoAppDestroy();
     48:    Routines: TaoCreate(); TaoDestroy();
     49:    Routines: TaoAppSetObjectiveAndGradientRoutine();
     50:    Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     51:    Routines: TaoSetOptions();
     52:    Routines: TaoAppSetInitialSolutionVec();
     53:    Routines: TaoSolveApplication();
     54:    Routines: TaoGetSolutionStatus(); TaoGetKSP();
     55:    Processors: 1
     56: T*/

     58: /*
     59:    User-defined application context - contains data needed by the
     60:    application-provided call-back routines, FormFunction(),
     61:    FormGradient(), and FormHessian().
     62: */

     64: typedef struct {
     65:    PetscReal  param;      /* nonlinearity parameter */
     66:    int     mx, my;     /* discretization in x- and y-directions */
     67:    int     ndim;       /* problem dimension */
     68:    Vec     s, y, xvec; /* work space for computing Hessian */
     69:    PetscReal  hx, hy;     /* mesh spacing in x- and y-directions */
     70: } AppCtx;

     72: /* -------- User-defined Routines --------- */

     74: int FormInitialGuess(AppCtx*,Vec);
     75: int FormFunction(TAO_APPLICATION,Vec,double*,void*);
     76: int FormGradient(TAO_APPLICATION,Vec,Vec,void*);
     77: int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*, MatStructure *,void*);
     78: int HessianProductMat(Mat,Vec,Vec);
     79: int HessianProduct(void*,Vec,Vec);
     80: int MatrixFreeHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*);
     81: int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void *);

     85: int main(int argc,char **argv)
     86: {
     87:   int         info;                /* used to check for functions returning nonzeros */
     88:   int         mx=10;               /* discretization in x-direction */
     89:   int         my=10;               /* discretization in y-direction */
     90:   Vec         x;                   /* solution, gradient vectors */
     91:   PetscTruth  flg;                 /* A return value when checking for use options */
     92:   TAO_SOLVER  tao;                 /* TAO_SOLVER solver context */
     93:   TAO_APPLICATION eptorsionapp;    /* The PETSc application */
     94:   Mat         H;                   /* Hessian matrix */
     95:   int         iter;                /* iteration information */
     96:   double      ff,gnorm;
     97:   TaoTerminateReason reason;
     98:   KSP         ksp;                 /* PETSc Krylov subspace solver */
     99:   PC          pc;                  /* PETSc preconditioner */
    100:   AppCtx      user;                /* application context */
    101:   int         size;                /* number of processes */
    102:   double      one=1.0;

    105:   /* Initialize TAO,PETSc */
    106:   PetscInitialize(&argc,&argv,(char *)0,help);
    107:   TaoInitialize(&argc,&argv,(char *)0,help);

    109:   /* Optional:  Check  that only one processor is being used. */
    110:   info = MPI_Comm_size(MPI_COMM_WORLD,&size); CHKERRQ(info);
    111:   if (size >1) {
    112:     PetscPrintf(PETSC_COMM_SELF,"This example is intended for single processor use!\n");
    113:     PetscPrintf(PETSC_COMM_SELF,"Try the example eptorsion2!\n");
    114:     SETERRQ(1,"Incorrect number of processors");
    115:   }

    117:   /* Specify default parameters for the problem, check for command-line overrides */
    118:   user.param = 5.0;
    119:   info = PetscOptionsGetInt(TAO_NULL,"-my",&my,&flg); CHKERRQ(info);
    120:   info = PetscOptionsGetInt(TAO_NULL,"-mx",&mx,&flg); CHKERRQ(info);
    121:   info = PetscOptionsGetReal(TAO_NULL,"-par",&user.param,&flg); CHKERRQ(info);

    123:   user.ndim = mx * my; user.mx = mx; user.my = my;
    124:   user.hx = one/(mx+1); user.hy = one/(my+1);

    127:   /* Allocate vectors */
    128:   info = VecCreateSeq(MPI_COMM_SELF,user.ndim,&user.y); CHKERRQ(info);
    129:   info = VecDuplicate(user.y,&user.s); CHKERRQ(info);
    130:   info = VecDuplicate(user.y,&x); CHKERRQ(info);

    132:   /* The TAO code begins here */

    134:   /* Create TAO solver and set desired solution method */
    135:   info = TaoCreate(MPI_COMM_SELF,"tao_lmvm",&tao); CHKERRQ(info);
    136:   info = TaoApplicationCreate(PETSC_COMM_SELF,&eptorsionapp); CHKERRQ(info);

    138:   /* Set solution vector with an initial guess */
    139:   info = FormInitialGuess(&user,x); CHKERRQ(info);
    140:   info = TaoAppSetInitialSolutionVec(eptorsionapp,x); CHKERRQ(info);

    142:   /* Set routine for function and gradient evaluation */
    143:   info = TaoAppSetObjectiveAndGradientRoutine(eptorsionapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);

    145:   /* From command line options, determine if using matrix-free hessian */
    146:   info = PetscOptionsHasName(TAO_NULL,"-my_tao_mf",&flg); CHKERRQ(info);
    147:   if (flg) {
    148:     info = MatCreateShell(MPI_COMM_SELF,user.ndim,user.ndim,user.ndim,
    149:                           user.ndim,(void*)&user,&H); CHKERRQ(info);
    150:     info = MatShellSetOperation(H,MATOP_MULT,(void(*)())HessianProductMat); CHKERRQ
    151: (info);
    152:     info = MatSetOption(H,MAT_SYMMETRIC); CHKERRQ(info);

    154:     info = TaoAppSetHessianRoutine(eptorsionapp,MatrixFreeHessian,(void *)&user); CHKERRQ(info);
    155:     info = TaoAppSetHessianMat(eptorsionapp,H,H); CHKERRQ(info);

    157:     /* Set null preconditioner.  Alternatively, set user-provided
    158:        preconditioner or explicitly form preconditioning matrix */
    159:     info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
    160:     if (ksp){
    161:       info = KSPGetPC(ksp,&pc); CHKERRQ(info);
    162:       info = PCSetType(pc,PCNONE); CHKERRQ(info);
    163:     }
    164:   } else {
    165:     info = MatCreateSeqAIJ(MPI_COMM_SELF,user.ndim,user.ndim,5,TAO_NULL,&H); CHKERRQ(info);
    166:     info = MatSetOption(H,MAT_SYMMETRIC); CHKERRQ(info);
    167:     info = TaoAppSetHessianRoutine(eptorsionapp,FormHessian,(void *)&user); CHKERRQ(info);
    168:     info = TaoAppSetHessianMat(eptorsionapp,H,H); CHKERRQ(info);
    169:   }

    171:   /* Check for any TAO command line options */
    172:   info = TaoSetOptions(eptorsionapp,tao); CHKERRQ(info);

    175:   /* Modify the PETSc KSP structure */
    176:   info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
    177:   if (ksp) {
    178:     info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
    179:   }

    181:   /* SOLVE THE APPLICATION */
    182:   info = TaoSolveApplication(eptorsionapp,tao);  CHKERRQ(info);

    184:   /*
    185:      To View TAO solver information use
    186:       info = TaoView(tao); CHKERRQ(info);
    187:   */

    189:   /* Get information on termination */
    190:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
    191:   if (reason <= 0){
    192:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
    193:     PetscPrintf(MPI_COMM_WORLD,"Iter: %d,   f: %4.2e,  residual: %4.2e\n",iter,ff,gnorm); CHKERRQ(info);
    194:   }

    196:   /* Free TAO data structures */
    197:   info = TaoDestroy(tao); CHKERRQ(info);
    198:   info = TaoAppDestroy(eptorsionapp); CHKERRQ(info);

    200:   /* Free PETSc data structures */
    201:   info = VecDestroy(user.s); CHKERRQ(info);
    202:   info = VecDestroy(user.y); CHKERRQ(info);
    203:   info = VecDestroy(x); CHKERRQ(info);
    204:   info = MatDestroy(H); CHKERRQ(info);

    206:   /* Finalize TAO, PETSc */
    207:   TaoFinalize();
    208:   PetscFinalize();

    210:   return 0;
    211: }

    213: /* ------------------------------------------------------------------- */
    216: /*
    217:     FormInitialGuess - Computes an initial approximation to the solution.

    219:     Input Parameters:
    220: .   user - user-defined application context
    221: .   X    - vector

    223:     Output Parameters:
    224: .   X    - vector
    225: */
    226: int FormInitialGuess(AppCtx *user,Vec X)
    227: {
    228:   double hx = user->hx, hy = user->hy, temp;
    229:   PetscScalar val;
    230:   int    info, i, j, k, nx = user->mx, ny = user->my;

    232:   /* Compute initial guess */
    233:   for (j=0; j<ny; j++) {
    234:     temp = TaoMin(j+1,ny-j)*hy;
    235:     for (i=0; i<nx; i++) {
    236:       k   = nx*j + i;
    237:       val = TaoMin((TaoMin(i+1,nx-i))*hx,temp);
    238:       info = VecSetValues(X,1,&k,&val,ADD_VALUES); CHKERRQ(info);
    239:     }
    240:   }

    242:   return 0;
    243: }
    244: /* ------------------------------------------------------------------- */
    247: /*
    248:    FormFunctionGradient - Evaluates the function and corresponding gradient.
    249:
    250:    Input Parameters:
    251:    tao - the TAO_APPLICATION context
    252:    X   - the input vector
    253:    ptr - optional user-defined context, as set by TaoSetFunction()

    255:    Output Parameters:
    256:    f   - the newly evaluated function
    257:    G   - the newly evaluated gradient
    258: */
    259: int FormFunctionGradient(TAO_APPLICATION tao,Vec X,double *f,Vec G,void *ptr)
    260: {
    261:   int info;
    262:   info = FormFunction(tao,X,f,ptr);CHKERRQ(info);
    263:   info = FormGradient(tao,X,G,ptr);CHKERRQ(info);
    264:   return 0;
    265: }
    266: /* ------------------------------------------------------------------- */
    269: /*
    270:    FormFunction - Evaluates the function, f(X).

    272:    Input Parameters:
    273: .  taoapp - the TAO_APPLICATION context
    274: .  X   - the input vector
    275: .  ptr - optional user-defined context, as set by TaoSetFunction()

    277:    Output Parameters:
    278: .  f    - the newly evaluated function
    279: */
    280: int FormFunction(TAO_APPLICATION taoapp,Vec X,double *f,void *ptr)
    281: {
    282:   AppCtx *user = (AppCtx *) ptr;
    283:   double hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
    284:   double zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
    285:   double v, cdiv3 = user->param/three;
    286:   PetscScalar *x;
    287:   int    info, nx = user->mx, ny = user->my, i, j, k;

    289:   /* Get pointer to vector data */
    290:   info = VecGetArray(X,&x); CHKERRQ(info);

    292:   /* Compute function contributions over the lower triangular elements */
    293:   for (j=-1; j<ny; j++) {
    294:     for (i=-1; i<nx; i++) {
    295:       k = nx*j + i;
    296:       v = zero;
    297:       vr = zero;
    298:       vt = zero;
    299:       if (i >= 0 && j >= 0) v = x[k];
    300:       if (i < nx-1 && j > -1) vr = x[k+1];
    301:       if (i > -1 && j < ny-1) vt = x[k+nx];
    302:       dvdx = (vr-v)/hx;
    303:       dvdy = (vt-v)/hy;
    304:       fquad += dvdx*dvdx + dvdy*dvdy;
    305:       flin -= cdiv3*(v+vr+vt);
    306:     }
    307:   }

    309:   /* Compute function contributions over the upper triangular elements */
    310:   for (j=0; j<=ny; j++) {
    311:     for (i=0; i<=nx; i++) {
    312:       k = nx*j + i;
    313:       vb = zero;
    314:       vl = zero;
    315:       v = zero;
    316:       if (i < nx && j > 0) vb = x[k-nx];
    317:       if (i > 0 && j < ny) vl = x[k-1];
    318:       if (i < nx && j < ny) v = x[k];
    319:       dvdx = (v-vl)/hx;
    320:       dvdy = (v-vb)/hy;
    321:       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
    322:       flin = flin - cdiv3*(vb+vl+v);
    323:     }
    324:   }

    326:   /* Restore vector */
    327:   info = VecRestoreArray(X,&x); CHKERRQ(info);

    329:   /* Scale the function */
    330:   area = p5*hx*hy;
    331:   *f = area*(p5*fquad+flin);

    333:   info = PetscLogFlops(nx*ny*24); CHKERRQ(info);
    334:   return 0;
    335: }
    336: /* ------------------------------------------------------------------- */
    339: /*
    340:     FormGradient - Evaluates the gradient, G(X).

    342:     Input Parameters:
    343: .   taoapp  - the TAO_APPLICATION context
    344: .   X    - input vector
    345: .   ptr  - optional user-defined context
    346:
    347:     Output Parameters:
    348: .   G - vector containing the newly evaluated gradient
    349: */
    350: int FormGradient(TAO_APPLICATION taoapp,Vec X,Vec G,void *ptr)
    351: {
    352:   AppCtx *user = (AppCtx *) ptr;
    353:   PetscScalar zero=0.0, p5=0.5, three = 3.0, area, val, *x;
    354:   int    info, nx = user->mx, ny = user->my, ind, i, j, k;
    355:   double hx = user->hx, hy = user->hy;
    356:   double vb, vl, vr, vt, dvdx, dvdy;
    357:   double v, cdiv3 = user->param/three;

    359:   /* Initialize gradient to zero */
    360:   info = VecSet(G, zero); CHKERRQ(info);

    362:   /* Get pointer to vector data */
    363:   info = VecGetArray(X,&x); CHKERRQ(info);

    365:   /* Compute gradient contributions over the lower triangular elements */
    366:   for (j=-1; j<ny; j++) {
    367:     for (i=-1; i<nx; i++) {
    368:       k  = nx*j + i;
    369:       v  = zero;
    370:       vr = zero;
    371:       vt = zero;
    372:       if (i >= 0 && j >= 0)    v = x[k];
    373:       if (i < nx-1 && j > -1) vr = x[k+1];
    374:       if (i > -1 && j < ny-1) vt = x[k+nx];
    375:       dvdx = (vr-v)/hx;
    376:       dvdy = (vt-v)/hy;
    377:       if (i != -1 && j != -1) {
    378:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
    379:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    380:       }
    381:       if (i != nx-1 && j != -1) {
    382:         ind = k+1; val =  dvdx/hx - cdiv3;
    383:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    384:       }
    385:       if (i != -1 && j != ny-1) {
    386:         ind = k+nx; val = dvdy/hy - cdiv3;
    387:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    388:       }
    389:     }
    390:   }

    392:   /* Compute gradient contributions over the upper triangular elements */
    393:   for (j=0; j<=ny; j++) {
    394:     for (i=0; i<=nx; i++) {
    395:       k = nx*j + i;
    396:       vb = zero;
    397:       vl = zero;
    398:       v  = zero;
    399:       if (i < nx && j > 0) vb = x[k-nx];
    400:       if (i > 0 && j < ny) vl = x[k-1];
    401:       if (i < nx && j < ny) v = x[k];
    402:       dvdx = (v-vl)/hx;
    403:       dvdy = (v-vb)/hy;
    404:       if (i != nx && j != 0) {
    405:         ind = k-nx; val = - dvdy/hy - cdiv3;
    406:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    407:       }
    408:       if (i != 0 && j != ny) {
    409:         ind = k-1; val =  - dvdx/hx - cdiv3;
    410:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    411:       }
    412:       if (i != nx && j != ny) {
    413:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
    414:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    415:       }
    416:     }
    417:   }

    419:   /* Restore vector */
    420:   info = VecRestoreArray(X,&x); CHKERRQ(info);

    422:   /* Assemble gradient vector */
    423:   info = VecAssemblyBegin(G); CHKERRQ(info);
    424:   info = VecAssemblyEnd(G); CHKERRQ(info);

    426:   /* Scale the gradient */
    427:   area = p5*hx*hy;
    428:   info = VecScale(G, area); CHKERRQ(info);
    429:
    430:   info = PetscLogFlops(nx*ny*24); CHKERRQ(info);
    431:   return 0;
    432: }

    434: /* ------------------------------------------------------------------- */
    437: /*
    438:    FormHessian - Forms the Hessian matrix.

    440:    Input Parameters:
    441: .  taoapp - the TAO_APPLICATION context
    442: .  X    - the input vector
    443: .  ptr  - optional user-defined context, as set by TaoSetHessian()
    444:
    445:    Output Parameters:
    446: .  H     - Hessian matrix
    447: .  PrecH - optionally different preconditioning Hessian
    448: .  flag  - flag indicating matrix structure

    450:    Notes:
    451:    This routine is intended simply as an example of the interface
    452:    to a Hessian evaluation routine.  Since this example compute the
    453:    Hessian a column at a time, it is not particularly efficient and
    454:    is not recommended.
    455: */
    456: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *HH,Mat *Hpre, MatStructure *flg, void *ptr)
    457: {
    458:   AppCtx     *user = (AppCtx *) ptr;
    459:   int        i, j, info, ndim = user->ndim;
    460:   PetscScalar  *y, zero = 0.0, one = 1.0;
    461:   Mat H=*HH;
    462:   *Hpre = H;
    463:   PetscTruth assembled;

    465:   /* Set location of vector */
    466:   user->xvec = X;

    468:   /* Initialize Hessian entries and work vector to zero */
    469:   info = MatAssembled(H,&assembled); CHKERRQ(info);
    470:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

    472:   info = VecSet(user->s, zero); CHKERRQ(info);

    474:   /* Loop over matrix columns to compute entries of the Hessian */
    475:   for (j=0; j<ndim; j++) {

    477:     info = VecSetValues(user->s,1,&j,&one,INSERT_VALUES); CHKERRQ(info);
    478:     info = VecAssemblyBegin(user->s); CHKERRQ(info);
    479:     info = VecAssemblyEnd(user->s); CHKERRQ(info);

    481:     info = HessianProduct(ptr,user->s,user->y); CHKERRQ(info);

    483:     info = VecSetValues(user->s,1,&j,&zero,INSERT_VALUES); CHKERRQ(info);
    484:     info = VecAssemblyBegin(user->s); CHKERRQ(info);
    485:     info = VecAssemblyEnd(user->s); CHKERRQ(info);

    487:     info = VecGetArray(user->y,&y); CHKERRQ(info);
    488:     for (i=0; i<ndim; i++) {
    489:       if (y[i] != zero) {
    490:         info = MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES); CHKERRQ(info);
    491:       }
    492:     }
    493:     info = VecRestoreArray(user->y,&y); CHKERRQ(info);

    495:   }

    497:   *flg=SAME_NONZERO_PATTERN;

    499:   /* Assemble matrix  */
    500:   info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    501:   info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

    503:   return 0;
    504: }

    507: /* ------------------------------------------------------------------- */
    510: /*
    511:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
    512:    products.
    513:
    514:    Input Parameters:
    515: .  taoapp - the TAO_APPLICATION context
    516: .  X    - the input vector
    517: .  ptr  - optional user-defined context, as set by TaoSetHessian()
    518:
    519:    Output Parameters:
    520: .  H     - Hessian matrix
    521: .  PrecH - optionally different preconditioning Hessian
    522: .  flag  - flag indicating matrix structure
    523: */
    524: int MatrixFreeHessian(TAO_APPLICATION taoapp,Vec X,Mat *H,Mat *PrecH,
    525:                       MatStructure *flag,void *ptr)
    526: {
    527:   AppCtx     *user = (AppCtx *) ptr;

    529:   /* Sets location of vector for use in computing matrix-vector products
    530:      of the form H(X)*y  */

    532:   user->xvec = X;
    533:   return 0;
    534: }

    536: /* ------------------------------------------------------------------- */
    539: /*
    540:    HessianProductMat - Computes the matrix-vector product
    541:    y = mat*svec.

    543:    Input Parameters:
    544: .  mat  - input matrix
    545: .  svec - input vector

    547:    Output Parameters:
    548: .  y    - solution vector
    549: */
    550: int HessianProductMat(Mat mat,Vec svec,Vec y)
    551: {
    552:   void *ptr;
    553:   MatShellGetContext(mat,&ptr);
    554:   HessianProduct(ptr,svec,y);

    556:   return 0;
    557: }

    559: /* ------------------------------------------------------------------- */
    562: /*
    563:    Hessian Product - Computes the matrix-vector product:
    564:    y = f''(x)*svec.

    566:    Input Parameters
    567: .  ptr  - pointer to the user-defined context
    568: .  svec - input vector

    570:    Output Parameters:
    571: .  y    - product vector
    572: */
    573: int HessianProduct(void *ptr,Vec svec,Vec y)
    574: {
    575:   AppCtx *user = (AppCtx *)ptr;
    576:   PetscScalar p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
    577:   double v, vb, vl, vr, vt, hxhx, hyhy;
    578:   int    nx, ny, i, j, k, info, ind;

    580:   nx   = user->mx;
    581:   ny   = user->my;
    582:   hx   = user->hx;
    583:   hy   = user->hy;
    584:   hxhx = one/(hx*hx);
    585:   hyhy = one/(hy*hy);

    587:   /* Get pointers to vector data */
    588:   info = VecGetArray(user->xvec,&x); CHKERRQ(info);
    589:   info = VecGetArray(svec,&s); CHKERRQ(info);

    591:   /* Initialize product vector to zero */
    592:   info = VecSet(y, zero); CHKERRQ(info);

    594:   /* Compute f''(x)*s over the lower triangular elements */
    595:   for (j=-1; j<ny; j++) {
    596:     for (i=-1; i<nx; i++) {
    597:        k = nx*j + i;
    598:        v = zero;
    599:        vr = zero;
    600:        vt = zero;
    601:        if (i != -1 && j != -1) v = s[k];
    602:        if (i != nx-1 && j != -1) {
    603:          vr = s[k+1];
    604:          ind = k+1; val = hxhx*(vr-v);
    605:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    606:        }
    607:        if (i != -1 && j != ny-1) {
    608:          vt = s[k+nx];
    609:          ind = k+nx; val = hyhy*(vt-v);
    610:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    611:        }
    612:        if (i != -1 && j != -1) {
    613:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
    614:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    615:        }
    616:      }
    617:    }
    618:
    619:   /* Compute f''(x)*s over the upper triangular elements */
    620:   for (j=0; j<=ny; j++) {
    621:     for (i=0; i<=nx; i++) {
    622:        k = nx*j + i;
    623:        v = zero;
    624:        vl = zero;
    625:        vb = zero;
    626:        if (i != nx && j != ny) v = s[k];
    627:        if (i != nx && j != 0) {
    628:          vb = s[k-nx];
    629:          ind = k-nx; val = hyhy*(vb-v);
    630:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    631:        }
    632:        if (i != 0 && j != ny) {
    633:          vl = s[k-1];
    634:          ind = k-1; val = hxhx*(vl-v);
    635:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    636:        }
    637:        if (i != nx && j != ny) {
    638:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
    639:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    640:        }
    641:     }
    642:   }
    643:   /* Restore vector data */
    644:   info = VecRestoreArray(svec,&s); CHKERRQ(info);
    645:   info = VecRestoreArray(user->xvec,&x); CHKERRQ(info);

    647:   /* Assemble vector */
    648:   info = VecAssemblyBegin(y); CHKERRQ(info);
    649:   info = VecAssemblyEnd(y); CHKERRQ(info);
    650:
    651:   /* Scale resulting vector by area */
    652:   area = p5*hx*hy;
    653:   info = VecScale(y, area); CHKERRQ(info);

    655:   PetscLogFlops(nx*ny*18);
    656:
    657:   return 0;
    658: }