Actual source code: ex18.c
2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
3: Uses 2-dimensional distributed arrays.\n\
4: A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
5: \n\
6: Solves the linear systems via multilevel methods \n\
7: \n\
8: The command line\n\
9: options are:\n\
10: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
11: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
12: -beta <beta>, where <beta> indicates the exponent in T \n\n";
14: /*T
15: Concepts: SNES^solving a system of nonlinear equations
16: Concepts: DA^using distributed arrays
17: Concepts: multigrid;
18: Processors: n
19: T*/
21: /*
22:
23: This example models the partial differential equation
24:
25: - Div(alpha* T^beta (GRAD T)) = 0.
26:
27: where beta = 2.5 and alpha = 1.0
28:
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
30:
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding. The degrees of freedom are cell centered.
33:
34: A finite volume approximation with the usual 5-point stencil
35: is used to discretize the boundary value problem to obtain a
36: nonlinear system of equations.
38: This code was contributed by David Keyes
39:
40: */
42: #include petscsnes.h
43: #include petscda.h
44: #include petscmg.h
45: #include petscdmmg.h
47: /* User-defined application context */
49: typedef struct {
50: PetscReal tleft,tright; /* Dirichlet boundary conditions */
51: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
52: } AppCtx;
54: #define POWFLOP 5 /* assume a pow() takes five flops */
62: int main(int argc,char **argv)
63: {
64: DMMG *dmmg;
65: SNES snes;
66: AppCtx user;
68: PetscInt its,lits;
69: PetscReal litspit;
70: DA da;
72: PetscInitialize(&argc,&argv,PETSC_NULL,help);
74: /* set problem parameters */
75: user.tleft = 1.0;
76: user.tright = 0.1;
77: user.beta = 2.5;
78: PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
79: PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
80: PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
81: user.bm1 = user.beta - 1.0;
82: user.coef = user.beta/2.0;
85: /*
86: Create the multilevel DA data structure
87: */
88: DMMGCreate(PETSC_COMM_WORLD,3,&user,&dmmg);
90: /*
91: Set the DA (grid structure) for the grids.
92: */
93: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
94: DMMGSetDM(dmmg,(DM)da);
95: DADestroy(da);
97: /*
98: Create the nonlinear solver, and tell the DMMG structure to use it
99: */
100: DMMGSetSNES(dmmg,FormFunction,FormJacobian);
102: /*
103: PreLoadBegin() means that the following section of code is run twice. The first time
104: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
105: The second time through (the actually timed code) the maximum iterations is set to 10
106: Preload of the executable is done to eliminate from the timing the time spent bring the
107: executable into memory from disk (paging in).
108: */
109: PreLoadBegin(PETSC_TRUE,"Solve");
110: DMMGSetInitialGuess(dmmg,FormInitialGuess);
111: DMMGSolve(dmmg);
112: PreLoadEnd();
113: snes = DMMGGetSNES(dmmg);
114: SNESGetIterationNumber(snes,&its);
115: SNESGetNumberLinearIterations(snes,&lits);
116: litspit = ((PetscReal)lits)/((PetscReal)its);
117: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
119: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);
121: DMMGDestroy(dmmg);
122: PetscFinalize();
124: return 0;
125: }
126: /* -------------------- Form initial approximation ----------------- */
129: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
130: {
131: AppCtx *user = (AppCtx*)dmmg->user;
132: PetscInt i,j,xs,ys,xm,ym;
134: PetscReal tleft = user->tleft;
135: PetscScalar **x;
139: /* Get ghost points */
140: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
141: DAVecGetArray((DA)dmmg->dm,X,&x);
143: /* Compute initial guess */
144: for (j=ys; j<ys+ym; j++) {
145: for (i=xs; i<xs+xm; i++) {
146: x[j][i] = tleft;
147: }
148: }
149: DAVecRestoreArray((DA)dmmg->dm,X,&x);
150: return(0);
151: }
152: /* -------------------- Evaluate Function F(x) --------------------- */
155: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
156: {
157: DMMG dmmg = (DMMG)ptr;
158: AppCtx *user = (AppCtx*)dmmg->user;
160: PetscInt i,j,mx,my,xs,ys,xm,ym;
161: PetscScalar zero = 0.0,one = 1.0;
162: PetscScalar hx,hy,hxdhy,hydhx;
163: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
164: PetscScalar tleft,tright,beta;
165: PetscScalar **x,**f;
166: Vec localX;
169: DAGetLocalVector((DA)dmmg->dm,&localX);
170: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
171: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
172: hxdhy = hx/hy; hydhx = hy/hx;
173: tleft = user->tleft; tright = user->tright;
174: beta = user->beta;
175:
176: /* Get ghost points */
177: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
178: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
179: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
180: DAVecGetArray((DA)dmmg->dm,localX,&x);
181: DAVecGetArray((DA)dmmg->dm,F,&f);
183: /* Evaluate function */
184: for (j=ys; j<ys+ym; j++) {
185: for (i=xs; i<xs+xm; i++) {
186: t0 = x[j][i];
188: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
190: /* general interior volume */
192: tw = x[j][i-1];
193: aw = 0.5*(t0 + tw);
194: dw = PetscPowScalar(aw,beta);
195: fw = dw*(t0 - tw);
197: te = x[j][i+1];
198: ae = 0.5*(t0 + te);
199: de = PetscPowScalar(ae,beta);
200: fe = de*(te - t0);
202: ts = x[j-1][i];
203: as = 0.5*(t0 + ts);
204: ds = PetscPowScalar(as,beta);
205: fs = ds*(t0 - ts);
206:
207: tn = x[j+1][i];
208: an = 0.5*(t0 + tn);
209: dn = PetscPowScalar(an,beta);
210: fn = dn*(tn - t0);
212: } else if (i == 0) {
214: /* left-hand boundary */
215: tw = tleft;
216: aw = 0.5*(t0 + tw);
217: dw = PetscPowScalar(aw,beta);
218: fw = dw*(t0 - tw);
220: te = x[j][i+1];
221: ae = 0.5*(t0 + te);
222: de = PetscPowScalar(ae,beta);
223: fe = de*(te - t0);
225: if (j > 0) {
226: ts = x[j-1][i];
227: as = 0.5*(t0 + ts);
228: ds = PetscPowScalar(as,beta);
229: fs = ds*(t0 - ts);
230: } else {
231: fs = zero;
232: }
234: if (j < my-1) {
235: tn = x[j+1][i];
236: an = 0.5*(t0 + tn);
237: dn = PetscPowScalar(an,beta);
238: fn = dn*(tn - t0);
239: } else {
240: fn = zero;
241: }
243: } else if (i == mx-1) {
245: /* right-hand boundary */
246: tw = x[j][i-1];
247: aw = 0.5*(t0 + tw);
248: dw = PetscPowScalar(aw,beta);
249: fw = dw*(t0 - tw);
250:
251: te = tright;
252: ae = 0.5*(t0 + te);
253: de = PetscPowScalar(ae,beta);
254: fe = de*(te - t0);
255:
256: if (j > 0) {
257: ts = x[j-1][i];
258: as = 0.5*(t0 + ts);
259: ds = PetscPowScalar(as,beta);
260: fs = ds*(t0 - ts);
261: } else {
262: fs = zero;
263: }
264:
265: if (j < my-1) {
266: tn = x[j+1][i];
267: an = 0.5*(t0 + tn);
268: dn = PetscPowScalar(an,beta);
269: fn = dn*(tn - t0);
270: } else {
271: fn = zero;
272: }
274: } else if (j == 0) {
276: /* bottom boundary,and i <> 0 or mx-1 */
277: tw = x[j][i-1];
278: aw = 0.5*(t0 + tw);
279: dw = PetscPowScalar(aw,beta);
280: fw = dw*(t0 - tw);
282: te = x[j][i+1];
283: ae = 0.5*(t0 + te);
284: de = PetscPowScalar(ae,beta);
285: fe = de*(te - t0);
287: fs = zero;
289: tn = x[j+1][i];
290: an = 0.5*(t0 + tn);
291: dn = PetscPowScalar(an,beta);
292: fn = dn*(tn - t0);
294: } else if (j == my-1) {
296: /* top boundary,and i <> 0 or mx-1 */
297: tw = x[j][i-1];
298: aw = 0.5*(t0 + tw);
299: dw = PetscPowScalar(aw,beta);
300: fw = dw*(t0 - tw);
302: te = x[j][i+1];
303: ae = 0.5*(t0 + te);
304: de = PetscPowScalar(ae,beta);
305: fe = de*(te - t0);
307: ts = x[j-1][i];
308: as = 0.5*(t0 + ts);
309: ds = PetscPowScalar(as,beta);
310: fs = ds*(t0 - ts);
312: fn = zero;
314: }
316: f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
318: }
319: }
320: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
321: DAVecRestoreArray((DA)dmmg->dm,F,&f);
322: DARestoreLocalVector((DA)dmmg->dm,&localX);
323: PetscLogFlops((22 + 4*POWFLOP)*ym*xm);
324: return(0);
325: }
326: /* -------------------- Evaluate Jacobian F(x) --------------------- */
329: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
330: {
331: DMMG dmmg = (DMMG)ptr;
332: AppCtx *user = (AppCtx*)dmmg->user;
333: Mat jac = *J;
335: PetscInt i,j,mx,my,xs,ys,xm,ym;
336: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
337: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
338: PetscScalar tleft,tright,beta,bm1,coef;
339: PetscScalar v[5],**x;
340: Vec localX;
341: MatStencil col[5],row;
344: DAGetLocalVector((DA)dmmg->dm,&localX);
345: *flg = SAME_NONZERO_PATTERN;
346: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
347: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
348: hxdhy = hx/hy; hydhx = hy/hx;
349: tleft = user->tleft; tright = user->tright;
350: beta = user->beta; bm1 = user->bm1; coef = user->coef;
352: /* Get ghost points */
353: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
354: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
355: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
356: DAVecGetArray((DA)dmmg->dm,localX,&x);
358: /* Evaluate Jacobian of function */
359: for (j=ys; j<ys+ym; j++) {
360: for (i=xs; i<xs+xm; i++) {
361: t0 = x[j][i];
363: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
365: /* general interior volume */
367: tw = x[j][i-1];
368: aw = 0.5*(t0 + tw);
369: bw = PetscPowScalar(aw,bm1);
370: /* dw = bw * aw */
371: dw = PetscPowScalar(aw,beta);
372: gw = coef*bw*(t0 - tw);
374: te = x[j][i+1];
375: ae = 0.5*(t0 + te);
376: be = PetscPowScalar(ae,bm1);
377: /* de = be * ae; */
378: de = PetscPowScalar(ae,beta);
379: ge = coef*be*(te - t0);
381: ts = x[j-1][i];
382: as = 0.5*(t0 + ts);
383: bs = PetscPowScalar(as,bm1);
384: /* ds = bs * as; */
385: ds = PetscPowScalar(as,beta);
386: gs = coef*bs*(t0 - ts);
387:
388: tn = x[j+1][i];
389: an = 0.5*(t0 + tn);
390: bn = PetscPowScalar(an,bm1);
391: /* dn = bn * an; */
392: dn = PetscPowScalar(an,beta);
393: gn = coef*bn*(tn - t0);
395: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
396: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
397: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
398: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
399: v[4] = - hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
400: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
402: } else if (i == 0) {
404: /* left-hand boundary */
405: tw = tleft;
406: aw = 0.5*(t0 + tw);
407: bw = PetscPowScalar(aw,bm1);
408: /* dw = bw * aw */
409: dw = PetscPowScalar(aw,beta);
410: gw = coef*bw*(t0 - tw);
411:
412: te = x[j][i + 1];
413: ae = 0.5*(t0 + te);
414: be = PetscPowScalar(ae,bm1);
415: /* de = be * ae; */
416: de = PetscPowScalar(ae,beta);
417: ge = coef*be*(te - t0);
418:
419: /* left-hand bottom boundary */
420: if (j == 0) {
422: tn = x[j+1][i];
423: an = 0.5*(t0 + tn);
424: bn = PetscPowScalar(an,bm1);
425: /* dn = bn * an; */
426: dn = PetscPowScalar(an,beta);
427: gn = coef*bn*(tn - t0);
428:
429: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
430: v[1] = - hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
431: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
432: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
433:
434: /* left-hand interior boundary */
435: } else if (j < my-1) {
437: ts = x[j-1][i];
438: as = 0.5*(t0 + ts);
439: bs = PetscPowScalar(as,bm1);
440: /* ds = bs * as; */
441: ds = PetscPowScalar(as,beta);
442: gs = coef*bs*(t0 - ts);
443:
444: tn = x[j+1][i];
445: an = 0.5*(t0 + tn);
446: bn = PetscPowScalar(an,bm1);
447: /* dn = bn * an; */
448: dn = PetscPowScalar(an,beta);
449: gn = coef*bn*(tn - t0);
450:
451: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
452: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
453: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
454: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
455: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
456: /* left-hand top boundary */
457: } else {
459: ts = x[j-1][i];
460: as = 0.5*(t0 + ts);
461: bs = PetscPowScalar(as,bm1);
462: /* ds = bs * as; */
463: ds = PetscPowScalar(as,beta);
464: gs = coef*bs*(t0 - ts);
465:
466: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
467: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
468: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
469: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
470: }
472: } else if (i == mx-1) {
473:
474: /* right-hand boundary */
475: tw = x[j][i-1];
476: aw = 0.5*(t0 + tw);
477: bw = PetscPowScalar(aw,bm1);
478: /* dw = bw * aw */
479: dw = PetscPowScalar(aw,beta);
480: gw = coef*bw*(t0 - tw);
481:
482: te = tright;
483: ae = 0.5*(t0 + te);
484: be = PetscPowScalar(ae,bm1);
485: /* de = be * ae; */
486: de = PetscPowScalar(ae,beta);
487: ge = coef*be*(te - t0);
488:
489: /* right-hand bottom boundary */
490: if (j == 0) {
492: tn = x[j+1][i];
493: an = 0.5*(t0 + tn);
494: bn = PetscPowScalar(an,bm1);
495: /* dn = bn * an; */
496: dn = PetscPowScalar(an,beta);
497: gn = coef*bn*(tn - t0);
498:
499: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
500: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
501: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
502: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
503:
504: /* right-hand interior boundary */
505: } else if (j < my-1) {
507: ts = x[j-1][i];
508: as = 0.5*(t0 + ts);
509: bs = PetscPowScalar(as,bm1);
510: /* ds = bs * as; */
511: ds = PetscPowScalar(as,beta);
512: gs = coef*bs*(t0 - ts);
513:
514: tn = x[j+1][i];
515: an = 0.5*(t0 + tn);
516: bn = PetscPowScalar(an,bm1);
517: /* dn = bn * an; */
518: dn = PetscPowScalar(an,beta);
519: gn = coef*bn*(tn - t0);
520:
521: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
522: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
523: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
524: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
525: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
526: /* right-hand top boundary */
527: } else {
529: ts = x[j-1][i];
530: as = 0.5*(t0 + ts);
531: bs = PetscPowScalar(as,bm1);
532: /* ds = bs * as; */
533: ds = PetscPowScalar(as,beta);
534: gs = coef*bs*(t0 - ts);
535:
536: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
537: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
538: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
539: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
540: }
542: /* bottom boundary,and i <> 0 or mx-1 */
543: } else if (j == 0) {
545: tw = x[j][i-1];
546: aw = 0.5*(t0 + tw);
547: bw = PetscPowScalar(aw,bm1);
548: /* dw = bw * aw */
549: dw = PetscPowScalar(aw,beta);
550: gw = coef*bw*(t0 - tw);
552: te = x[j][i+1];
553: ae = 0.5*(t0 + te);
554: be = PetscPowScalar(ae,bm1);
555: /* de = be * ae; */
556: de = PetscPowScalar(ae,beta);
557: ge = coef*be*(te - t0);
559: tn = x[j+1][i];
560: an = 0.5*(t0 + tn);
561: bn = PetscPowScalar(an,bm1);
562: /* dn = bn * an; */
563: dn = PetscPowScalar(an,beta);
564: gn = coef*bn*(tn - t0);
565:
566: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
567: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
568: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
569: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
570: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
571:
572: /* top boundary,and i <> 0 or mx-1 */
573: } else if (j == my-1) {
574:
575: tw = x[j][i-1];
576: aw = 0.5*(t0 + tw);
577: bw = PetscPowScalar(aw,bm1);
578: /* dw = bw * aw */
579: dw = PetscPowScalar(aw,beta);
580: gw = coef*bw*(t0 - tw);
582: te = x[j][i+1];
583: ae = 0.5*(t0 + te);
584: be = PetscPowScalar(ae,bm1);
585: /* de = be * ae; */
586: de = PetscPowScalar(ae,beta);
587: ge = coef*be*(te - t0);
589: ts = x[j-1][i];
590: as = 0.5*(t0 + ts);
591: bs = PetscPowScalar(as,bm1);
592: /* ds = bs * as; */
593: ds = PetscPowScalar(as,beta);
594: gs = coef*bs*(t0 - ts);
596: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
597: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
598: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
599: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
600: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
601:
602: }
603: }
604: }
605: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
606: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
607: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
608: DARestoreLocalVector((DA)dmmg->dm,&localX);
610: PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
611: return(0);
612: }