Actual source code: ex5.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
3: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";
9: /*T
10: Concepts: SNES^parallel Bratu example
11: Concepts: DMDA^using distributed arrays;
12: Concepts: IS coloirng types;
13: Processors: n
14: T*/
16: /* ------------------------------------------------------------------------
18: Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: the partial differential equation
21: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
23: with boundary conditions
25: u = 0 for x = 0, x = 1, y = 0, y = 1.
27: A finite difference approximation with the usual 5-point stencil
28: is used to discretize the boundary value problem to obtain a nonlinear
29: system of equations.
32: This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
33: as SNESSetDM() is provided. Example usage
35: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_levels 3 -pc_mg_galerkin -da_grid_x 17 -da_grid_y 17
36: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
38: or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
39: multigrid levels, it will be determined automatically based on the number of refinements done)
41: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin -snes_grid_sequence 3
42: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
44: ------------------------------------------------------------------------- */
46: /*
47: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
48: Include "petscsnes.h" so that we can use SNES solvers. Note that this
49: */
50: #include <petscdmda.h>
51: #include <petscsnes.h>
52: #include <petscmatlab.h>
54: /*
55: User-defined application context - contains data needed by the
56: application-provided call-back routines, FormJacobianLocal() and
57: FormFunctionLocal().
58: */
59: typedef struct {
60: PassiveReal param; /* test problem parameter */
61: } AppCtx;
63: /*
64: User-defined routines
65: */
66: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
67: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
68: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,MatStructure*,AppCtx*);
69: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
70: #if defined(PETSC_HAVE_MATLAB_ENGINE)
71: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
72: #endif
73: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);
77: int main(int argc,char **argv)
78: {
79: SNES snes; /* nonlinear solver */
80: Vec x; /* solution vector */
81: AppCtx user; /* user-defined work context */
82: PetscInt its; /* iterations for convergence */
84: PetscReal bratu_lambda_max = 6.81;
85: PetscReal bratu_lambda_min = 0.;
86: PetscBool flg = PETSC_FALSE;
87: DM da;
88: #if defined(PETSC_HAVE_MATLAB_ENGINE)
89: Vec r = NULL;
90: PetscBool matlab_function = PETSC_FALSE;
91: #endif
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Initialize program
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: PetscInitialize(&argc,&argv,(char*)0,help);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Initialize problem parameters
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: user.param = 6.0;
103: PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
104: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda, %g, is out of range, [%g, %g]", user.param, bratu_lambda_min, bratu_lambda_max);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create nonlinear solver context
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: SNESCreate(PETSC_COMM_WORLD,&snes);
110: SNESSetGS(snes, NonlinearGS, NULL);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create distributed array (DMDA) to manage parallel grid and vectors
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
116: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
117: DMSetApplicationContext(da,&user);
118: SNESSetDM(snes,da);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Extract global vectors from DMDA; then duplicate for remaining
121: vectors that are the same types
122: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: DMCreateGlobalVector(da,&x);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Set local function evaluation routine
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
129: PetscOptionsGetBool(NULL,"-fd",&flg,NULL);
130: if (!flg) {
131: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
132: }
134: PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
135: if (flg) {
136: DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
137: }
139: #if defined(PETSC_HAVE_MATLAB_ENGINE)
140: PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
141: if (matlab_function) {
142: VecDuplicate(x,&r);
143: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
144: }
145: #endif
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Customize nonlinear solver; set runtime options
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: SNESSetFromOptions(snes);
152: FormInitialGuess(da,&user,x);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Solve nonlinear system
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: SNESSolve(snes,NULL,x);
158: SNESGetIterationNumber(snes,&its);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Free work space. All PETSc objects should be destroyed when they
165: are no longer needed.
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: #if defined(PETSC_HAVE_MATLAB_ENGINE)
168: VecDestroy(&r);
169: #endif
170: VecDestroy(&x);
171: SNESDestroy(&snes);
172: DMDestroy(&da);
173: PetscFinalize();
174: return 0;
175: }
176: /* ------------------------------------------------------------------- */
179: /*
180: FormInitialGuess - Forms initial approximation.
182: Input Parameters:
183: user - user-defined application context
184: X - vector
186: Output Parameter:
187: X - vector
188: */
189: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
190: {
191: PetscInt i,j,Mx,My,xs,ys,xm,ym;
193: PetscReal lambda,temp1,temp,hx,hy;
194: PetscScalar **x;
197: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
199: lambda = user->param;
200: hx = 1.0/(PetscReal)(Mx-1);
201: hy = 1.0/(PetscReal)(My-1);
202: temp1 = lambda/(lambda + 1.0);
204: /*
205: Get a pointer to vector data.
206: - For default PETSc vectors, VecGetArray() returns a pointer to
207: the data array. Otherwise, the routine is implementation dependent.
208: - You MUST call VecRestoreArray() when you no longer need access to
209: the array.
210: */
211: DMDAVecGetArray(da,X,&x);
213: /*
214: Get local grid boundaries (for 2-dimensional DMDA):
215: xs, ys - starting grid indices (no ghost points)
216: xm, ym - widths of local grid (no ghost points)
218: */
219: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
221: /*
222: Compute initial guess over the locally owned part of the grid
223: */
224: for (j=ys; j<ys+ym; j++) {
225: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
226: for (i=xs; i<xs+xm; i++) {
227: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
228: /* boundary conditions are all zero Dirichlet */
229: x[j][i] = 0.0;
230: } else {
231: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
232: }
233: }
234: }
236: /*
237: Restore vector
238: */
239: DMDAVecRestoreArray(da,X,&x);
240: return(0);
241: }
242: /* ------------------------------------------------------------------- */
245: /*
246: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
249: */
250: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
251: {
253: PetscInt i,j;
254: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
255: PetscScalar u,ue,uw,un,us,uxx,uyy;
258: lambda = user->param;
259: hx = 1.0/(PetscReal)(info->mx-1);
260: hy = 1.0/(PetscReal)(info->my-1);
261: sc = hx*hy*lambda;
262: hxdhy = hx/hy;
263: hydhx = hy/hx;
264: /*
265: Compute function over the locally owned part of the grid
266: */
267: for (j=info->ys; j<info->ys+info->ym; j++) {
268: for (i=info->xs; i<info->xs+info->xm; i++) {
269: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
270: f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
271: } else {
272: u = x[j][i];
273: uw = x[j][i-1];
274: ue = x[j][i+1];
275: un = x[j-1][i];
276: us = x[j+1][i];
278: if (i-1 == 0) uw = 0.;
279: if (i+1 == info->mx-1) ue = 0.;
280: if (j-1 == 0) un = 0.;
281: if (j+1 == info->my-1) us = 0.;
283: uxx = (2.0*u - uw - ue)*hydhx;
284: uyy = (2.0*u - un - us)*hxdhy;
285: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
286: }
287: }
288: }
289: PetscLogFlops(11.0*info->ym*info->xm);
290: return(0);
291: }
296: /*
297: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
300: */
301: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
302: {
304: PetscInt i,j;
305: PetscReal lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
306: PetscScalar u,ue,uw,un,us,uxux,uyuy;
307: MPI_Comm comm;
310: *obj = 0;
311: PetscObjectGetComm((PetscObject)info,&comm);
312: lambda = user->param;
313: hx = 1.0/(PetscReal)(info->mx-1);
314: hy = 1.0/(PetscReal)(info->my-1);
315: sc = hx*hy*lambda;
316: hxdhy = hx/hy;
317: hydhx = hy/hx;
318: /*
319: Compute function over the locally owned part of the grid
320: */
321: for (j=info->ys; j<info->ys+info->ym; j++) {
322: for (i=info->xs; i<info->xs+info->xm; i++) {
323: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
324: lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
325: } else {
326: u = x[j][i];
327: uw = x[j][i-1];
328: ue = x[j][i+1];
329: un = x[j-1][i];
330: us = x[j+1][i];
332: if (i-1 == 0) uw = 0.;
333: if (i+1 == info->mx-1) ue = 0.;
334: if (j-1 == 0) un = 0.;
335: if (j+1 == info->my-1) us = 0.;
337: /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
339: uxux = u*(2.*u - ue - uw)*hydhx;
340: uyuy = u*(2.*u - un - us)*hxdhy;
342: lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
343: }
344: }
345: }
346: PetscLogFlops(12.0*info->ym*info->xm);
347: MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
348: return(0);
349: }
353: /*
354: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
355: */
356: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jacpre,Mat jac,MatStructure *flg,AppCtx *user)
357: {
359: PetscInt i,j,k;
360: MatStencil col[5],row;
361: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
364: lambda = user->param;
365: hx = 1.0/(PetscReal)(info->mx-1);
366: hy = 1.0/(PetscReal)(info->my-1);
367: sc = hx*hy*lambda;
368: hxdhy = hx/hy;
369: hydhx = hy/hx;
372: /*
373: Compute entries for the locally owned part of the Jacobian.
374: - Currently, all PETSc parallel matrix formats are partitioned by
375: contiguous chunks of rows across the processors.
376: - Each processor needs to insert only elements that it owns
377: locally (but any non-local elements will be sent to the
378: appropriate processor during matrix assembly).
379: - Here, we set all entries for a particular row at once.
380: - We can set matrix entries either using either
381: MatSetValuesLocal() or MatSetValues(), as discussed above.
382: */
383: for (j=info->ys; j<info->ys+info->ym; j++) {
384: for (i=info->xs; i<info->xs+info->xm; i++) {
385: row.j = j; row.i = i;
386: /* boundary points */
387: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
388: v[0] = 2.0*(hydhx + hxdhy);
389: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
390: } else {
391: k = 0;
392: /* interior grid points */
393: if (j-1 != 0) {
394: v[k] = -hxdhy;
395: col[k].j = j - 1; col[k].i = i;
396: k++;
397: }
398: if (i-1 != 0) {
399: v[k] = -hydhx;
400: col[k].j = j; col[k].i = i-1;
401: k++;
402: }
404: v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;
406: if (i+1 != info->mx-1) {
407: v[k] = -hydhx;
408: col[k].j = j; col[k].i = i+1;
409: k++;
410: }
411: if (j+1 != info->mx-1) {
412: v[k] = -hxdhy;
413: col[k].j = j + 1; col[k].i = i;
414: k++;
415: }
416: MatSetValuesStencil(jac,1,&row,k,col,v,INSERT_VALUES);
417: }
418: }
419: }
421: /*
422: Assemble matrix, using the 2-step process:
423: MatAssemblyBegin(), MatAssemblyEnd().
424: */
425: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
426: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
427: /*
428: Tell the matrix we will never add a new nonzero location to the
429: matrix. If we do, it will generate an error.
430: */
431: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
432: return(0);
433: }
435: #if defined(PETSC_HAVE_MATLAB_ENGINE)
438: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
439: {
440: AppCtx *user = (AppCtx*)ptr;
442: PetscInt Mx,My;
443: PetscReal lambda,hx,hy;
444: Vec localX,localF;
445: MPI_Comm comm;
446: DM da;
449: SNESGetDM(snes,&da);
450: DMGetLocalVector(da,&localX);
451: DMGetLocalVector(da,&localF);
452: PetscObjectSetName((PetscObject)localX,"localX");
453: PetscObjectSetName((PetscObject)localF,"localF");
454: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
456: lambda = user->param;
457: hx = 1.0/(PetscReal)(Mx-1);
458: hy = 1.0/(PetscReal)(My-1);
460: PetscObjectGetComm((PetscObject)snes,&comm);
461: /*
462: Scatter ghost points to local vector,using the 2-step process
463: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
464: By placing code between these two statements, computations can be
465: done while messages are in transition.
466: */
467: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
468: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
469: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
470: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
471: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
473: /*
474: Insert values into global vector
475: */
476: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
477: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
478: DMRestoreLocalVector(da,&localX);
479: DMRestoreLocalVector(da,&localF);
480: return(0);
481: }
482: #endif
484: /* ------------------------------------------------------------------- */
487: /*
488: Applies some sweeps on nonlinear Gauss-Seidel on each process
490: */
491: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
492: {
493: PetscInt i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
495: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
496: PetscScalar **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
497: PetscReal atol,rtol,stol;
498: DM da;
499: AppCtx *user;
500: Vec localX,localB;
503: tot_its = 0;
504: SNESGSGetSweeps(snes,&sweeps);
505: SNESGSGetTolerances(snes,&atol,&rtol,&stol,&its);
506: SNESGetDM(snes,&da);
507: DMGetApplicationContext(da,(void**)&user);
509: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
511: lambda = user->param;
512: hx = 1.0/(PetscReal)(Mx-1);
513: hy = 1.0/(PetscReal)(My-1);
514: sc = hx*hy*lambda;
515: hxdhy = hx/hy;
516: hydhx = hy/hx;
519: DMGetLocalVector(da,&localX);
520: if (B) {
521: DMGetLocalVector(da,&localB);
522: }
523: for (l=0; l<sweeps; l++) {
525: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
526: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
527: if (B) {
528: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
529: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
530: }
531: /*
532: Get a pointer to vector data.
533: - For default PETSc vectors, VecGetArray() returns a pointer to
534: the data array. Otherwise, the routine is implementation dependent.
535: - You MUST call VecRestoreArray() when you no longer need access to
536: the array.
537: */
538: DMDAVecGetArray(da,localX,&x);
539: if (B) DMDAVecGetArray(da,localB,&b);
540: /*
541: Get local grid boundaries (for 2-dimensional DMDA):
542: xs, ys - starting grid indices (no ghost points)
543: xm, ym - widths of local grid (no ghost points)
544: */
545: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
547: for (j=ys; j<ys+ym; j++) {
548: for (i=xs; i<xs+xm; i++) {
549: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
550: /* boundary conditions are all zero Dirichlet */
551: x[j][i] = 0.0;
552: } else {
553: if (B) bij = b[j][i];
554: else bij = 0.;
556: u = x[j][i];
557: un = x[j-1][i];
558: us = x[j+1][i];
559: ue = x[j][i-1];
560: uw = x[j][i+1];
562: for (k=0; k<its; k++) {
563: eu = PetscExpScalar(u);
564: uxx = (2.0*u - ue - uw)*hydhx;
565: uyy = (2.0*u - un - us)*hxdhy;
566: F = uxx + uyy - sc*eu - bij;
567: if (k == 0) F0 = F;
568: J = 2.0*(hydhx + hxdhy) - sc*eu;
569: y = F/J;
570: u -= y;
571: tot_its++;
573: if (atol > PetscAbsReal(PetscRealPart(F)) ||
574: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
575: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
576: break;
577: }
578: }
579: x[j][i] = u;
580: }
581: }
582: }
583: /*
584: Restore vector
585: */
586: DMDAVecRestoreArray(da,localX,&x);
587: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
588: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
589: }
590: PetscLogFlops(tot_its*(21.0));
591: DMRestoreLocalVector(da,&localX);
592: if (B) {
593: DMDAVecRestoreArray(da,localB,&b);
594: DMRestoreLocalVector(da,&localB);
595: }
596: return(0);
597: }