Actual source code: eptorsion1.c

  1: /*$Id$*/

  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(); TaoAppGetKSP();
 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:    PetscInt   mx, my;     /* discretization in x- and y-directions */
 67:    PetscInt   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:   PetscInt      mx=10;               /* discretization in x-direction */
 89:   PetscInt      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:   TaoInt      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);


124:   PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
125:   PetscPrintf(PETSC_COMM_SELF,"mx: %d     my: %d   \n\n",mx,my);  
126:   user.ndim = mx * my; user.mx = mx; user.my = my;

128:   user.hx = one/(mx+1); user.hy = one/(my+1);


131:   /* Allocate vectors */
132:   info = VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y); CHKERRQ(info);
133:   info = VecDuplicate(user.y,&user.s); CHKERRQ(info);
134:   info = VecDuplicate(user.y,&x); CHKERRQ(info);

136:   /* The TAO code begins here */

138:   /* Create TAO solver and set desired solution method */
139:   info = TaoCreate(PETSC_COMM_SELF,"tao_lmvm",&tao); CHKERRQ(info);
140:   info = TaoApplicationCreate(PETSC_COMM_SELF,&eptorsionapp); CHKERRQ(info);

142:   /* Set solution vector with an initial guess */
143:   info = FormInitialGuess(&user,x); CHKERRQ(info);
144:   info = TaoAppSetInitialSolutionVec(eptorsionapp,x); CHKERRQ(info);

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

149:   /* From command line options, determine if using matrix-free hessian */
150:   info = PetscOptionsHasName(TAO_NULL,"-my_tao_mf",&flg); CHKERRQ(info);
151:   if (flg) {
152:     info = MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,
153:                           user.ndim,(void*)&user,&H); CHKERRQ(info);
154:     info = MatShellSetOperation(H,MATOP_MULT,(void(*)())HessianProductMat); CHKERRQ
155: (info);
156:     info = MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);

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

161:     /* Set null preconditioner.  Alternatively, set user-provided 
162:        preconditioner or explicitly form preconditioning matrix */
163:     info = TaoAppGetKSP(eptorsionapp,&ksp); CHKERRQ(info);
164:     if (ksp){
165:       info = KSPGetPC(ksp,&pc); CHKERRQ(info);
166:       info = PCSetType(pc,PCNONE); CHKERRQ(info);
167:     }
168:   } else {

170:     info = MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,TAO_NULL,&H); CHKERRQ(info);
171:     info = MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);

173:     info = TaoAppSetHessianRoutine(eptorsionapp,FormHessian,(void *)&user); CHKERRQ(info);
174:     info = TaoAppSetHessianMat(eptorsionapp,H,H); CHKERRQ(info); 
175:   }

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


181:   /* Modify the PETSc KSP structure */
182:   info = TaoAppGetKSP(eptorsionapp,&ksp); CHKERRQ(info);
183:   if (ksp) {                                              
184:     info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
185:   }

187:   /* SOLVE THE APPLICATION */
188:   info = TaoSolveApplication(eptorsionapp,tao);  CHKERRQ(info);
189:   if (ksp) {
190:     KSPView(ksp,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(info);
191:   }

193:   /* 
194:      To View TAO solver information use
195:       info = TaoView(tao); CHKERRQ(info);
196:   */

198:   /* Get information on termination */
199:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
200:   if (reason <= 0){
201:     PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
202:     PetscPrintf(PETSC_COMM_WORLD,"Iter: %d,   f: %4.2e,  residual: %4.2e\n",iter,ff,gnorm); CHKERRQ(info);
203:   }

205:   /* Free TAO data structures */
206:   info = TaoDestroy(tao); CHKERRQ(info);
207:   info = TaoAppDestroy(eptorsionapp); CHKERRQ(info);

209:   /* Free PETSc data structures */
210:   info = VecDestroy(user.s); CHKERRQ(info);
211:   info = VecDestroy(user.y); CHKERRQ(info);
212:   info = VecDestroy(x); CHKERRQ(info);
213:   info = MatDestroy(H); CHKERRQ(info);

215:   /* Finalize TAO, PETSc */
216:   TaoFinalize();
217:   PetscFinalize();

219:   return 0;
220: }

222: /* ------------------------------------------------------------------- */
225: /* 
226:     FormInitialGuess - Computes an initial approximation to the solution.

228:     Input Parameters:
229: .   user - user-defined application context
230: .   X    - vector

232:     Output Parameters:
233: .   X    - vector
234: */
235: int FormInitialGuess(AppCtx *user,Vec X)
236: {
237:   double hx = user->hx, hy = user->hy, temp;
238:   PetscScalar val;
239:   int    info;
240:   PetscInt i, j, k, nx = user->mx, ny = user->my;

242:   /* Compute initial guess */
243:   for (j=0; j<ny; j++) {
244:     temp = TaoMin(j+1,ny-j)*hy;
245:     for (i=0; i<nx; i++) {
246:       k   = nx*j + i;
247:       val = TaoMin((TaoMin(i+1,nx-i))*hx,temp);
248:       info = VecSetValues(X,1,&k,&val,ADD_VALUES); CHKERRQ(info);
249:     }
250:   }
251:   info = VecAssemblyBegin(X); CHKERRQ(info);
252:   info = VecAssemblyEnd(X); CHKERRQ(info);
253:   return 0;
254: }
255: /* ------------------------------------------------------------------- */
258: /* 
259:    FormFunctionGradient - Evaluates the function and corresponding gradient.
260:     
261:    Input Parameters:
262:    tao - the TAO_APPLICATION context
263:    X   - the input vector 
264:    ptr - optional user-defined context, as set by TaoSetFunction()

266:    Output Parameters:
267:    f   - the newly evaluated function
268:    G   - the newly evaluated gradient
269: */
270: int FormFunctionGradient(TAO_APPLICATION tao,Vec X,double *f,Vec G,void *ptr)
271: {
272:   int info;
273:   info = FormFunction(tao,X,f,ptr);CHKERRQ(info);
274:   info = FormGradient(tao,X,G,ptr);CHKERRQ(info);
275:   return 0;
276: }
277: /* ------------------------------------------------------------------- */
280: /* 
281:    FormFunction - Evaluates the function, f(X).

283:    Input Parameters:
284: .  taoapp - the TAO_APPLICATION context
285: .  X   - the input vector 
286: .  ptr - optional user-defined context, as set by TaoSetFunction()

288:    Output Parameters:
289: .  f    - the newly evaluated function
290: */
291: int FormFunction(TAO_APPLICATION taoapp,Vec X,double *f,void *ptr)
292: {
293:   AppCtx *user = (AppCtx *) ptr;
294:   double hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
295:   double zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
296:   double v, cdiv3 = user->param/three;
297:   PetscScalar *x;
298:   int    info;
299:   PetscInt  nx = user->mx, ny = user->my, i, j, k;

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

304:   /* Compute function contributions over the lower triangular elements */
305:   for (j=-1; j<ny; j++) {
306:     for (i=-1; i<nx; i++) {
307:       k = nx*j + i;
308:       v = zero;
309:       vr = zero;
310:       vt = zero;
311:       if (i >= 0 && j >= 0) v = x[k];
312:       if (i < nx-1 && j > -1) vr = x[k+1];
313:       if (i > -1 && j < ny-1) vt = x[k+nx];
314:       dvdx = (vr-v)/hx;
315:       dvdy = (vt-v)/hy;
316:       fquad += dvdx*dvdx + dvdy*dvdy;
317:       flin -= cdiv3*(v+vr+vt);
318:     }
319:   }

321:   /* Compute function contributions over the upper triangular elements */
322:   for (j=0; j<=ny; j++) {
323:     for (i=0; i<=nx; i++) {
324:       k = nx*j + i;
325:       vb = zero;
326:       vl = zero;
327:       v = zero;
328:       if (i < nx && j > 0) vb = x[k-nx];
329:       if (i > 0 && j < ny) vl = x[k-1];
330:       if (i < nx && j < ny) v = x[k];
331:       dvdx = (v-vl)/hx;
332:       dvdy = (v-vb)/hy;
333:       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
334:       flin = flin - cdiv3*(vb+vl+v);
335:     }
336:   }

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

341:   /* Scale the function */
342:   area = p5*hx*hy;
343:   *f = area*(p5*fquad+flin);

345:   info = PetscLogFlops(nx*ny*24); CHKERRQ(info);
346:   return 0;
347: }
348: /* ------------------------------------------------------------------- */
351: /*  
352:     FormGradient - Evaluates the gradient, G(X).              

354:     Input Parameters:
355: .   taoapp  - the TAO_APPLICATION context
356: .   X    - input vector
357: .   ptr  - optional user-defined context
358:     
359:     Output Parameters:
360: .   G - vector containing the newly evaluated gradient
361: */
362: int FormGradient(TAO_APPLICATION taoapp,Vec X,Vec G,void *ptr)
363: {
364:   AppCtx *user = (AppCtx *) ptr;
365:   PetscScalar zero=0.0, p5=0.5, three = 3.0, area, val, *x;
366:   int    info;
367:   PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
368:   double hx = user->hx, hy = user->hy;
369:   double vb, vl, vr, vt, dvdx, dvdy;
370:   double v, cdiv3 = user->param/three;

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

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

378:   /* Compute gradient contributions over the lower triangular elements */
379:   for (j=-1; j<ny; j++) {
380:     for (i=-1; i<nx; i++) {
381:       k  = nx*j + i;
382:       v  = zero;
383:       vr = zero;
384:       vt = zero;
385:       if (i >= 0 && j >= 0)    v = x[k];
386:       if (i < nx-1 && j > -1) vr = x[k+1];
387:       if (i > -1 && j < ny-1) vt = x[k+nx];
388:       dvdx = (vr-v)/hx;
389:       dvdy = (vt-v)/hy;
390:       if (i != -1 && j != -1) {
391:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
392:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
393:       }
394:       if (i != nx-1 && j != -1) {
395:         ind = k+1; val =  dvdx/hx - cdiv3;
396:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
397:       }
398:       if (i != -1 && j != ny-1) {
399:         ind = k+nx; val = dvdy/hy - cdiv3;
400:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
401:       }
402:     }
403:   }

405:   /* Compute gradient contributions over the upper triangular elements */
406:   for (j=0; j<=ny; j++) {
407:     for (i=0; i<=nx; i++) {
408:       k = nx*j + i;
409:       vb = zero;
410:       vl = zero;
411:       v  = zero;
412:       if (i < nx && j > 0) vb = x[k-nx];
413:       if (i > 0 && j < ny) vl = x[k-1];
414:       if (i < nx && j < ny) v = x[k];
415:       dvdx = (v-vl)/hx;
416:       dvdy = (v-vb)/hy;
417:       if (i != nx && j != 0) {
418:         ind = k-nx; val = - dvdy/hy - cdiv3;
419:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
420:       }
421:       if (i != 0 && j != ny) {
422:         ind = k-1; val =  - dvdx/hx - cdiv3;
423:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
424:       }
425:       if (i != nx && j != ny) {
426:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
427:         info = VecSetValues(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
428:       }
429:     }
430:   }

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

435:   /* Assemble gradient vector */
436:   info = VecAssemblyBegin(G); CHKERRQ(info);
437:   info = VecAssemblyEnd(G); CHKERRQ(info);

439:   /* Scale the gradient */
440:   area = p5*hx*hy;
441:   info = VecScale(G, area); CHKERRQ(info);
442:   
443:   info = PetscLogFlops(nx*ny*24); CHKERRQ(info);
444:   return 0;
445: }

447: /* ------------------------------------------------------------------- */
450: /* 
451:    FormHessian - Forms the Hessian matrix.

453:    Input Parameters:
454: .  taoapp - the TAO_APPLICATION context
455: .  X    - the input vector
456: .  ptr  - optional user-defined context, as set by TaoSetHessian()
457:    
458:    Output Parameters:
459: .  H     - Hessian matrix
460: .  PrecH - optionally different preconditioning Hessian
461: .  flag  - flag indicating matrix structure

463:    Notes:
464:    This routine is intended simply as an example of the interface
465:    to a Hessian evaluation routine.  Since this example compute the
466:    Hessian a column at a time, it is not particularly efficient and
467:    is not recommended.
468: */
469: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *HH,Mat *Hpre, MatStructure *flg, void *ptr)
470: {
471:   AppCtx     *user = (AppCtx *) ptr;
472:   int info;
473:   PetscInt   i,j, ndim = user->ndim;
474:   PetscScalar  *y, zero = 0.0, one = 1.0;
475:   Mat H=*HH;
476:   *Hpre = H;
477:   PetscTruth assembled;

479:   /* Set location of vector */
480:   user->xvec = X;

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

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

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

491:     info = VecSetValues(user->s,1,&j,&one,INSERT_VALUES); CHKERRQ(info);
492:     info = VecAssemblyBegin(user->s); CHKERRQ(info);
493:     info = VecAssemblyEnd(user->s); CHKERRQ(info);

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

497:     info = VecSetValues(user->s,1,&j,&zero,INSERT_VALUES); CHKERRQ(info);
498:     info = VecAssemblyBegin(user->s); CHKERRQ(info);
499:     info = VecAssemblyEnd(user->s); CHKERRQ(info);

501:     info = VecGetArray(user->y,&y); CHKERRQ(info);
502:     for (i=0; i<ndim; i++) {
503:       if (y[i] != zero) {
504:         info = MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES); CHKERRQ(info);
505:       }
506:     }
507:     info = VecRestoreArray(user->y,&y); CHKERRQ(info);

509:   }

511:   *flg=SAME_NONZERO_PATTERN;

513:   /* Assemble matrix  */
514:   info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
515:   info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

517:   return 0;
518: }


521: /* ------------------------------------------------------------------- */
524: /* 
525:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
526:    products.
527:     
528:    Input Parameters:
529: .  taoapp - the TAO_APPLICATION context
530: .  X    - the input vector
531: .  ptr  - optional user-defined context, as set by TaoSetHessian()
532:    
533:    Output Parameters:
534: .  H     - Hessian matrix
535: .  PrecH - optionally different preconditioning Hessian
536: .  flag  - flag indicating matrix structure
537: */
538: int MatrixFreeHessian(TAO_APPLICATION taoapp,Vec X,Mat *H,Mat *PrecH,
539:                       MatStructure *flag,void *ptr)
540: {
541:   AppCtx     *user = (AppCtx *) ptr;

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

546:   user->xvec = X;   
547:   return 0;
548: }

550: /* ------------------------------------------------------------------- */
553: /* 
554:    HessianProductMat - Computes the matrix-vector product
555:    y = mat*svec.

557:    Input Parameters:
558: .  mat  - input matrix
559: .  svec - input vector

561:    Output Parameters:
562: .  y    - solution vector
563: */
564: int HessianProductMat(Mat mat,Vec svec,Vec y)
565: {
566:   void *ptr;
567:   MatShellGetContext(mat,&ptr);
568:   HessianProduct(ptr,svec,y);

570:   return 0;
571: }

573: /* ------------------------------------------------------------------- */
576: /* 
577:    Hessian Product - Computes the matrix-vector product: 
578:    y = f''(x)*svec.

580:    Input Parameters
581: .  ptr  - pointer to the user-defined context
582: .  svec - input vector

584:    Output Parameters:
585: .  y    - product vector
586: */
587: int HessianProduct(void *ptr,Vec svec,Vec y)
588: {
589:   AppCtx *user = (AppCtx *)ptr;
590:   PetscScalar p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
591:   double v, vb, vl, vr, vt, hxhx, hyhy;
592:   int    info;
593:   PetscInt nx, ny, i, j, k, ind;

595:   nx   = user->mx;
596:   ny   = user->my;
597:   hx   = user->hx;
598:   hy   = user->hy;
599:   hxhx = one/(hx*hx);
600:   hyhy = one/(hy*hy);

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

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

609:   /* Compute f''(x)*s over the lower triangular elements */
610:   for (j=-1; j<ny; j++) {
611:     for (i=-1; i<nx; i++) {
612:        k = nx*j + i;
613:        v = zero;
614:        vr = zero;
615:        vt = zero;
616:        if (i != -1 && j != -1) v = s[k];
617:        if (i != nx-1 && j != -1) {
618:          vr = s[k+1];
619:          ind = k+1; val = hxhx*(vr-v);
620:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
621:        }
622:        if (i != -1 && j != ny-1) {
623:          vt = s[k+nx];
624:          ind = k+nx; val = hyhy*(vt-v);
625:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
626:        }
627:        if (i != -1 && j != -1) {
628:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
629:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
630:        }
631:      }
632:    }
633:   
634:   /* Compute f''(x)*s over the upper triangular elements */
635:   for (j=0; j<=ny; j++) {
636:     for (i=0; i<=nx; i++) {
637:        k = nx*j + i;
638:        v = zero;
639:        vl = zero;
640:        vb = zero;
641:        if (i != nx && j != ny) v = s[k];
642:        if (i != nx && j != 0) {
643:          vb = s[k-nx];
644:          ind = k-nx; val = hyhy*(vb-v);
645:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
646:        }
647:        if (i != 0 && j != ny) {
648:          vl = s[k-1];
649:          ind = k-1; val = hxhx*(vl-v);
650:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
651:        }
652:        if (i != nx && j != ny) {
653:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
654:          info = VecSetValues(y,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
655:        }
656:     }
657:   }
658:   /* Restore vector data */
659:   info = VecRestoreArray(svec,&s); CHKERRQ(info);
660:   info = VecRestoreArray(user->xvec,&x); CHKERRQ(info);

662:   /* Assemble vector */
663:   info = VecAssemblyBegin(y); CHKERRQ(info);
664:   info = VecAssemblyEnd(y); CHKERRQ(info);
665:  
666:   /* Scale resulting vector by area */
667:   area = p5*hx*hy;
668:   info = VecScale(y, area); CHKERRQ(info);

670:   info = PetscLogFlops(nx*ny*18); CHKERRQ(info);
671:   
672:   return 0;
673: }