Actual source code: ex4.c
petsc-3.3-p6 2013-02-11
2: /* NOTE: THIS PROGRAM HAS NOT YET BEEN SET UP IN TUTORIAL STYLE. */
4: static char help[] ="This program demonstrates use of the SNES package to solve systems of nonlinear equations on a single processor.\n\
5: Both of these examples employ\n\
6: sparse storage of the Jacobian matrices and are taken from the MINPACK-2\n\
7: test suite. By default the Bratu (SFI - solid fuel ignition) test problem\n\
8: is solved. The command line options are:\n\
9: -cavity : Solve FDC (flow in a driven cavity) problem\n\
10: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
11: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
12: problem FDC: <parameter> = Reynolds number (par > 0)\n\
13: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
14: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
16: /*
17: 1) Solid Fuel Ignition (SFI) problem. This problem is modeled by
18: the partial differential equation
19:
20: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
21:
22: with boundary conditions
23:
24: u = 0 for x = 0, x = 1, y = 0, y = 1.
25:
26: A finite difference approximation with the usual 5-point stencil
27: is used to discretize the boundary value problem to obtain a nonlinear
28: system of equations.
29:
30: 2) Flow in a Driven Cavity (FDC) problem. The problem is
31: formulated as a boundary value problem, which is discretized by
32: standard finite difference approximations to obtain a system of
33: nonlinear equations.
34: */
36: #include <petscsnes.h>
38: typedef struct {
39: PetscReal param; /* test problem parameter */
40: PetscInt mx; /* Discretization in x-direction */
41: PetscInt my; /* Discretization in y-direction */
42: } AppCtx;
44: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
45: FormFunction1(SNES,Vec,Vec,void*),
46: FormInitialGuess1(AppCtx*,Vec);
47: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
48: FormFunction2(SNES,Vec,Vec,void*),
49: FormInitialGuess2(AppCtx*,Vec);
50:
53: int main(int argc, char **argv)
54: {
55: SNES snes; /* SNES context */
56: const SNESType method = SNESLS; /* default nonlinear solution method */
57: Vec x, r; /* solution, residual vectors */
58: Mat J; /* Jacobian matrix */
59: AppCtx user; /* user-defined application context */
60: PetscDraw draw; /* drawing context */
61: PetscInt its, N, nfails;
63: PetscBool cavity;
64: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
65: PetscScalar *xvalues;
67: PetscInitialize(&argc, &argv,(char *)0,help);
68: /* PetscDrawOpenX(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw); */
69: PetscDrawCreate(PETSC_COMM_WORLD,0,"Solution",300,0,300,300,&draw);
70: PetscDrawSetType(draw,PETSC_DRAW_X);
72: user.mx = 4;
73: user.my = 4;
74: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
75: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
76: PetscOptionsHasName(PETSC_NULL,"-cavity",&cavity);
77: if (cavity) user.param = 100.0;
78: else user.param = 6.0;
79: PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
80: if (!cavity && (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min)) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
81: N = user.mx*user.my;
82:
83: /* Set up data structures */
84: VecCreateSeq(PETSC_COMM_SELF,N,&x);
85: VecDuplicate(x,&r);
86: MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,PETSC_NULL,&J);
88: /* Create nonlinear solver */
89: SNESCreate(PETSC_COMM_WORLD,&snes);
90: SNESSetType(snes,method);
92: /* Set various routines and compute initial guess for nonlinear solver */
93: if (cavity){
94: FormInitialGuess2(&user,x);
95: SNESSetFunction(snes,r,FormFunction2,(void*)&user);
96: SNESSetJacobian(snes,J,J,FormJacobian2,(void*)&user);
97: } else {
98: FormInitialGuess1(&user,x);
99: SNESSetFunction(snes,r,FormFunction1,(void*)&user);
100: SNESSetJacobian(snes,J,J,FormJacobian1,(void*)&user);
101: }
103: /* Set options and solve nonlinear system */
104: SNESSetFromOptions(snes);
105: SNESSolve(snes,PETSC_NULL,x);
106: SNESGetIterationNumber(snes,&its);
107: SNESGetNonlinearStepFailures(snes,&nfails);
109: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D, ",its);
110: PetscPrintf(PETSC_COMM_SELF,"number of unsuccessful steps = %D\n\n",nfails);
111: VecGetArray(x,&xvalues);
112: PetscDrawTensorContour(draw,user.mx,user.my,0,0,xvalues);
113: VecRestoreArray(x,&xvalues);
115: /* Free data structures */
116: VecDestroy(&x); VecDestroy(&r);
117: MatDestroy(&J); SNESDestroy(&snes);
118: PetscDrawDestroy(&draw);
119: PetscFinalize();
121: return 0;
122: }
123: /* ------------------------------------------------------------------ */
124: /* Bratu (Solid Fuel Ignition) Test Problem */
125: /* ------------------------------------------------------------------ */
127: /* -------------------- Form initial approximation ----------------- */
131: PetscErrorCode FormInitialGuess1(AppCtx *user,Vec X)
132: {
133: PetscInt i, j, row, mx, my;
135: PetscReal lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
136: PetscScalar *x;
138: mx = user->mx;
139: my = user->my;
140: lambda = user->param;
142: hx = 1.0 / (PetscReal)(mx-1);
143: hy = 1.0 / (PetscReal)(my-1);
144: sc = hx*hy;
145: hxdhy = hx/hy;
146: hydhx = hy/hx;
148: VecGetArray(X,&x);
149: temp1 = lambda/(lambda + 1.0);
150: for (j=0; j<my; j++) {
151: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
152: for (i=0; i<mx; i++) {
153: row = i + j*mx;
154: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
155: x[row] = 0.0;
156: continue;
157: }
158: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
159: }
160: }
161: VecRestoreArray(X,&x);
162: return 0;
163: }
164: /* -------------------- Evaluate Function F(x) --------------------- */
165:
168: PetscErrorCode FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
169: {
170: AppCtx *user = (AppCtx*)ptr;
171: PetscInt i, j, row, mx, my;
173: PetscReal two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx;
174: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc,*x,*f;
176: mx = user->mx;
177: my = user->my;
178: lambda = user->param;
180: hx = one / (PetscReal)(mx-1);
181: hy = one / (PetscReal)(my-1);
182: sc = hx*hy;
183: hxdhy = hx/hy;
184: hydhx = hy/hx;
186: VecGetArray(X,&x);
187: VecGetArray(F,&f);
188: for (j=0; j<my; j++) {
189: for (i=0; i<mx; i++) {
190: row = i + j*mx;
191: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
192: f[row] = x[row];
193: continue;
194: }
195: u = x[row];
196: ub = x[row - mx];
197: ul = x[row - 1];
198: ut = x[row + mx];
199: ur = x[row + 1];
200: uxx = (-ur + two*u - ul)*hydhx;
201: uyy = (-ut + two*u - ub)*hxdhy;
202: f[row] = uxx + uyy - sc*lambda*exp(u);
203: }
204: }
205: VecRestoreArray(X,&x);
206: VecRestoreArray(F,&f);
207: return 0;
208: }
209: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
213: PetscErrorCode FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
214: {
215: AppCtx *user = (AppCtx*)ptr;
216: Mat jac = *J;
217: PetscInt i, j, row, mx, my, col[5];
219: PetscScalar two = 2.0, one = 1.0, lambda, v[5],sc, *x;
220: PetscReal hx, hy, hxdhy, hydhx;
223: mx = user->mx;
224: my = user->my;
225: lambda = user->param;
227: hx = 1.0 / (PetscReal)(mx-1);
228: hy = 1.0 / (PetscReal)(my-1);
229: sc = hx*hy;
230: hxdhy = hx/hy;
231: hydhx = hy/hx;
233: VecGetArray(X,&x);
234: for (j=0; j<my; j++) {
235: for (i=0; i<mx; i++) {
236: row = i + j*mx;
237: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
238: MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
239: continue;
240: }
241: v[0] = -hxdhy; col[0] = row - mx;
242: v[1] = -hydhx; col[1] = row - 1;
243: v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
244: v[3] = -hydhx; col[3] = row + 1;
245: v[4] = -hxdhy; col[4] = row + mx;
246: MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
247: }
248: }
249: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
250: VecRestoreArray(X,&x);
251: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
252: *flag = SAME_NONZERO_PATTERN;
253: return 0;
254: }
255: /* ------------------------------------------------------------------ */
256: /* Driven Cavity Test Problem */
257: /* ------------------------------------------------------------------ */
259: /* -------------------- Form initial approximation ----------------- */
263: PetscErrorCode FormInitialGuess2(AppCtx *user,Vec X)
264: {
266: PetscInt i, j, row, mx, my;
267: PetscScalar xx,yy,*x;
268: PetscReal hx, hy;
270: mx = user->mx;
271: my = user->my;
273: /* Test for invalid input parameters */
274: if ((mx <= 0) || (my <= 0)) SETERRQ(PETSC_COMM_SELF,1,0);
276: hx = 1.0 / (PetscReal)(mx-1);
277: hy = 1.0 / (PetscReal)(my-1);
279: VecGetArray(X,&x);
280: yy = 0.0;
281: for (j=0; j<my; j++) {
282: xx = 0.0;
283: for (i=0; i<mx; i++) {
284: row = i + j*mx;
285: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
286: x[row] = 0.0;
287: }
288: else {
289: x[row] = - xx*(1.0 - xx)*yy*(1.0 - yy);
290: }
291: xx = xx + hx;
292: }
293: yy = yy + hy;
294: }
295: VecRestoreArray(X,&x);
296: return 0;
297: }
298: /* -------------------- Evaluate Function F(x) --------------------- */
302: PetscErrorCode FormFunction2(SNES snes,Vec X,Vec F,void *pptr)
303: {
304: AppCtx *user = (AppCtx*)pptr;
305: PetscInt i, j, row, mx, my;
307: PetscScalar two = 2.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
308: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
309: PetscScalar *x,*f, hx2, hy2, hxhy2;
310: PetscReal hx, hy;
312: mx = user->mx;
313: my = user->my;
314: hx = 1.0 / (PetscReal)(mx-1);
315: hy = 1.0 / (PetscReal)(my-1);
316: hx2 = hx*hx;
317: hy2 = hy*hy;
318: hxhy2 = hx2*hy2;
319: rey = user->param;
321: VecGetArray(X,&x);
322: VecGetArray(F,&f);
323: for (j=0; j<my; j++) {
324: for (i=0; i<mx; i++) {
325: row = i + j*mx;
326: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
327: f[row] = x[row];
328: continue;
329: }
330: if (i == 1 || j == 1) {
331: pbl = zero;
332: }
333: else {
334: pbl = x[row-mx-1];
335: }
336: if (j == 1) {
337: pb = zero;
338: pbb = x[row];
339: }
340: else if (j == 2) {
341: pb = x[row-mx];
342: pbb = zero;
343: }
344: else {
345: pb = x[row-mx];
346: pbb = x[row-2*mx];
347: }
348: if (j == 1 || i == mx-2) {
349: pbr = zero;
350: }
351: else {
352: pbr = x[row-mx+1];
353: }
354: if (i == 1) {
355: pl = zero;
356: pll = x[row];
357: }
358: else if (i == 2) {
359: pl = x[row-1];
360: pll = zero;
361: }
362: else {
363: pl = x[row-1];
364: pll = x[row-2];
365: }
366: p = x[row];
367: if (i == mx-3) {
368: pr = x[row+1];
369: prr = zero;
370: }
371: else if (i == mx-2) {
372: pr = zero;
373: prr = x[row];
374: }
375: else {
376: pr = x[row+1];
377: prr = x[row+2];
378: }
379: if (j == my-2 || i == 1) {
380: ptl = zero;
381: }
382: else {
383: ptl = x[row+mx-1];
384: }
385: if (j == my-3) {
386: pt = x[row+mx];
387: ptt = zero;
388: }
389: else if (j == my-2) {
390: pt = zero;
391: ptt = x[row] + two*hy;
392: }
393: else {
394: pt = x[row+mx];
395: ptt = x[row+2*mx];
396: }
397: if (j == my-2 || i == mx-2) {
398: ptr = zero;
399: }
400: else {
401: ptr = x[row+mx+1];
402: }
404: dpdy = (pt - pb)/(two*hy);
405: dpdx = (pr - pl)/(two*hx);
407: /* Laplacians at each point in the 5 point stencil */
408: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
409: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
410: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
411: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
412: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
414: f[row] = hxhy2*((prlap - two*plap + pllap)/hx2
415: + (ptlap - two*plap + pblap)/hy2
416: - rey*(dpdy*(prlap - pllap)/(two*hx)
417: - dpdx*(ptlap - pblap)/(two*hy)));
418: }
419: }
420: VecRestoreArray(X,&x);
421: VecRestoreArray(F,&f);
422: return 0;
423: }
424: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
428: PetscErrorCode FormJacobian2(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *pptr)
429: {
430: AppCtx *user = (AppCtx*)pptr;
431: PetscInt i, j, row, mx, my, col;
433: PetscScalar two = 2.0, one = 1.0, zero = 0.0, pb, pbb,pbr, pl,pll,p,pr,prr;
434: PetscScalar ptl,pt,ptt,dpdy,dpdx,pblap,ptlap,rey,pbl,ptr,pllap,plap,prlap;
435: PetscScalar val,four = 4.0, three = 3.0,*x;
436: PetscReal hx, hy,hx2, hy2, hxhy2;
437: PetscBool assembled;
439: mx = user->mx;
440: my = user->my;
441: hx = 1.0 / (PetscReal)(mx-1);
442: hy = 1.0 / (PetscReal)(my-1);
443: hx2 = hx*hx;
444: hy2 = hy*hy;
445: hxhy2 = hx2*hy2;
446: rey = user->param;
448: MatAssembled(*J,&assembled);
449: if (assembled) {
450: MatZeroEntries(*J);
451: }
452: VecGetArray(X,&x);
453: for (j=0; j<my; j++) {
454: for (i=0; i<mx; i++) {
455: row = i + j*mx;
456: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
457: MatSetValues(*J,1,&row,1,&row,&one,ADD_VALUES);
458: continue;
459: }
460: if (i == 1 || j == 1) {
461: pbl = zero;
462: }
463: else {
464: pbl = x[row-mx-1];
465: }
466: if (j == 1) {
467: pb = zero;
468: pbb = x[row];
469: }
470: else if (j == 2) {
471: pb = x[row-mx];
472: pbb = zero;
473: }
474: else {
475: pb = x[row-mx];
476: pbb = x[row-2*mx];
477: }
478: if (j == 1 || i == mx-2) {
479: pbr = zero;
480: }
481: else {
482: pbr = x[row-mx+1];
483: }
484: if (i == 1) {
485: pl = zero;
486: pll = x[row];
487: }
488: else if (i == 2) {
489: pl = x[row-1];
490: pll = zero;
491: }
492: else {
493: pl = x[row-1];
494: pll = x[row-2];
495: }
496: p = x[row];
497: if (i == mx-3) {
498: pr = x[row+1];
499: prr = zero;
500: }
501: else if (i == mx-2) {
502: pr = zero;
503: prr = x[row];
504: }
505: else {
506: pr = x[row+1];
507: prr = x[row+2];
508: }
509: if (j == my-2 || i == 1) {
510: ptl = zero;
511: }
512: else {
513: ptl = x[row+mx-1];
514: }
515: if (j == my-3) {
516: pt = x[row+mx];
517: ptt = zero;
518: }
519: else if (j == my-2) {
520: pt = zero;
521: ptt = x[row] + two*hy;
522: }
523: else {
524: pt = x[row+mx];
525: ptt = x[row+2*mx];
526: }
527: if (j == my-2 || i == mx-2) {
528: ptr = zero;
529: }
530: else {
531: ptr = x[row+mx+1];
532: }
534: dpdy = (pt - pb)/(two*hy);
535: dpdx = (pr - pl)/(two*hx);
537: /* Laplacians at each point in the 5 point stencil */
538: pblap = (pbr - two*pb + pbl)/hx2 + (p - two*pb + pbb)/hy2;
539: pllap = (p - two*pl + pll)/hx2 + (ptl - two*pl + pbl)/hy2;
540: plap = (pr - two*p + pl)/hx2 + (pt - two*p + pb)/hy2;
541: prlap = (prr - two*pr + p)/hx2 + (ptr - two*pr + pbr)/hy2;
542: ptlap = (ptr - two*pt + ptl)/hx2 + (ptt - two*pt + p)/hy2;
544: if (j > 2) {
545: val = hxhy2*(one/hy2/hy2 - rey*dpdx/hy2/(two*hy));
546: col = row - 2*mx;
547: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
548: }
549: if (i > 2) {
550: val = hxhy2*(one/hx2/hx2 + rey*dpdy/hx2/(two*hx));
551: col = row - 2;
552: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
553: }
554: if (i < mx-3) {
555: val = hxhy2*(one/hx2/hx2 - rey*dpdy/hx2/(two*hx));
556: col = row + 2;
557: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
558: }
559: if (j < my-3) {
560: val = hxhy2*(one/hy2/hy2 + rey*dpdx/hy2/(two*hy));
561: col = row + 2*mx;
562: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
563: }
564: if (i != 1 && j != 1) {
565: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
566: col = row - mx - 1;
567: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
568: }
569: if (j != 1 && i != mx-2) {
570: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
571: col = row - mx + 1;
572: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
573: }
574: if (j != my-2 && i != 1) {
575: val = hxhy2*(two/hy2/hx2 + rey*(dpdy/hy2/(two*hx) + dpdx/hx2/(two*hy)));
576: col = row + mx - 1;
577: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
578: }
579: if (j != my-2 && i != mx-2) {
580: val = hxhy2*(two/hy2/hx2 - rey*(dpdy/hy2/(two*hx) - dpdx/hx2/(two*hy)));
581: col = row + mx + 1;
582: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
583: }
584: if (j != 1) {
585: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
586: + rey*((prlap - pllap)/(two*hx)/(two*hy)
587: + dpdx*(one/hx2 + one/hy2)/hy));
588: col = row - mx;
589: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
590: }
591: if (i != 1) {
592: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
593: - rey*((ptlap - pblap)/(two*hx)/(two*hy)
594: + dpdy*(one/hx2 + one/hy2)/hx));
595: col = row - 1;
596: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
597: }
598: if (i != mx-2) {
599: val = hxhy2*(-four*(one/hy2/hx2 + one/hx2/hx2)
600: + rey*((ptlap - pblap)/(two*hx)/(two*hy)
601: + dpdy*(one/hx2 + one/hy2)/hx));
602: col = row + 1;
603: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
604: }
605: if (j != my-2) {
606: val = hxhy2*(-four*(one/hy2/hx2 + one/hy2/hy2)
607: - rey*((prlap - pllap)/(two*hx)/(two*hy)
608: + dpdx*(one/hx2 + one/hy2)/hy));
609: col = row + mx;
610: MatSetValues(*J,1,&row,1,&col,&val,ADD_VALUES);
611: }
612: val = hxhy2*(two*(four/hx2/hy2 + three/hx2/hx2 + three/hy2/hy2));
613: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
614: if (j == 1) {
615: val = hxhy2*(one/hy2/hy2 - rey*(dpdx/hy2/(two*hy)));
616: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
617: }
618: if (i == 1) {
619: val = hxhy2*(one/hx2/hx2 + rey*(dpdy/hx2/(two*hx)));
620: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
621: }
622: if (i == mx-2) {
623: val = hxhy2*(one/hx2/hx2 - rey*(dpdy/hx2/(two*hx)));
624: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
625: }
626: if (j == my-2) {
627: val = hxhy2*(one/hy2/hy2 + rey*(dpdx/hy2/(two*hy)));
628: MatSetValues(*J,1,&row,1,&row,&val,ADD_VALUES);
629: }
630: }
631: }
632: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
633: VecRestoreArray(X,&x);
634: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
635: *flag = SAME_NONZERO_PATTERN;
636: return 0;
637: }