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