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