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: }