Actual source code: ex5.c
petsc-3.5.4 2015-05-23
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 <petscdm.h>
51: #include <petscdmda.h>
52: #include <petscsnes.h>
53: #include <petscmatlab.h>
55: /*
56: User-defined application context - contains data needed by the
57: application-provided call-back routines, FormJacobianLocal() and
58: FormFunctionLocal().
59: */
60: typedef struct {
61: PassiveReal param; /* test problem parameter */
62: } AppCtx;
64: /*
65: User-defined routines
66: */
67: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
68: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
69: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
70: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
71: #if defined(PETSC_HAVE_MATLAB_ENGINE)
72: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
73: #endif
74: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);
78: int main(int argc,char **argv)
79: {
80: SNES snes; /* nonlinear solver */
81: Vec x; /* solution vector */
82: AppCtx user; /* user-defined work context */
83: PetscInt its; /* iterations for convergence */
85: PetscReal bratu_lambda_max = 6.81;
86: PetscReal bratu_lambda_min = 0.;
87: PetscBool flg = PETSC_FALSE;
88: DM da;
89: #if defined(PETSC_HAVE_MATLAB_ENGINE)
90: Vec r = NULL;
91: PetscBool matlab_function = PETSC_FALSE;
92: #endif
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Initialize program
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscInitialize(&argc,&argv,(char*)0,help);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Initialize problem parameters
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: user.param = 6.0;
104: PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
105: 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);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create nonlinear solver context
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: SNESCreate(PETSC_COMM_WORLD,&snes);
111: SNESSetNGS(snes, NonlinearGS, NULL);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create distributed array (DMDA) to manage parallel grid and vectors
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
117: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
118: DMSetApplicationContext(da,&user);
119: SNESSetDM(snes,da);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Extract global vectors from DMDA; then duplicate for remaining
122: vectors that are the same types
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: DMCreateGlobalVector(da,&x);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Set local function evaluation routine
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
130: PetscOptionsGetBool(NULL,"-fd",&flg,NULL);
131: if (!flg) {
132: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
133: }
135: PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
136: if (flg) {
137: DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
138: }
140: #if defined(PETSC_HAVE_MATLAB_ENGINE)
141: PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
142: if (matlab_function) {
143: VecDuplicate(x,&r);
144: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
145: }
146: #endif
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Customize nonlinear solver; set runtime options
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: SNESSetFromOptions(snes);
153: FormInitialGuess(da,&user,x);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Solve nonlinear system
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: SNESSolve(snes,NULL,x);
159: SNESGetIterationNumber(snes,&its);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Free work space. All PETSc objects should be destroyed when they
166: are no longer needed.
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: #if defined(PETSC_HAVE_MATLAB_ENGINE)
169: VecDestroy(&r);
170: #endif
171: VecDestroy(&x);
172: SNESDestroy(&snes);
173: DMDestroy(&da);
174: PetscFinalize();
175: return 0;
176: }
177: /* ------------------------------------------------------------------- */
180: /*
181: FormInitialGuess - Forms initial approximation.
183: Input Parameters:
184: user - user-defined application context
185: X - vector
187: Output Parameter:
188: X - vector
189: */
190: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
191: {
192: PetscInt i,j,Mx,My,xs,ys,xm,ym;
194: PetscReal lambda,temp1,temp,hx,hy;
195: PetscScalar **x;
198: 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);
200: lambda = user->param;
201: hx = 1.0/(PetscReal)(Mx-1);
202: hy = 1.0/(PetscReal)(My-1);
203: temp1 = lambda/(lambda + 1.0);
205: /*
206: Get a pointer to vector data.
207: - For default PETSc vectors, VecGetArray() returns a pointer to
208: the data array. Otherwise, the routine is implementation dependent.
209: - You MUST call VecRestoreArray() when you no longer need access to
210: the array.
211: */
212: DMDAVecGetArray(da,X,&x);
214: /*
215: Get local grid boundaries (for 2-dimensional DMDA):
216: xs, ys - starting grid indices (no ghost points)
217: xm, ym - widths of local grid (no ghost points)
219: */
220: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
222: /*
223: Compute initial guess over the locally owned part of the grid
224: */
225: for (j=ys; j<ys+ym; j++) {
226: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
227: for (i=xs; i<xs+xm; i++) {
228: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
229: /* boundary conditions are all zero Dirichlet */
230: x[j][i] = 0.0;
231: } else {
232: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
233: }
234: }
235: }
237: /*
238: Restore vector
239: */
240: DMDAVecRestoreArray(da,X,&x);
241: return(0);
242: }
243: /* ------------------------------------------------------------------- */
246: /*
247: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
250: */
251: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
252: {
254: PetscInt i,j;
255: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
256: PetscScalar u,ue,uw,un,us,uxx,uyy;
259: lambda = user->param;
260: hx = 1.0/(PetscReal)(info->mx-1);
261: hy = 1.0/(PetscReal)(info->my-1);
262: sc = hx*hy*lambda;
263: hxdhy = hx/hy;
264: hydhx = hy/hx;
265: /*
266: Compute function over the locally owned part of the grid
267: */
268: for (j=info->ys; j<info->ys+info->ym; j++) {
269: for (i=info->xs; i<info->xs+info->xm; i++) {
270: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
271: f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
272: } else {
273: u = x[j][i];
274: uw = x[j][i-1];
275: ue = x[j][i+1];
276: un = x[j-1][i];
277: us = x[j+1][i];
279: if (i-1 == 0) uw = 0.;
280: if (i+1 == info->mx-1) ue = 0.;
281: if (j-1 == 0) un = 0.;
282: if (j+1 == info->my-1) us = 0.;
284: uxx = (2.0*u - uw - ue)*hydhx;
285: uyy = (2.0*u - un - us)*hxdhy;
286: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
287: }
288: }
289: }
290: PetscLogFlops(11.0*info->ym*info->xm);
291: return(0);
292: }
297: /*
298: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
301: */
302: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
303: {
305: PetscInt i,j;
306: PetscReal lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
307: PetscScalar u,ue,uw,un,us,uxux,uyuy;
308: MPI_Comm comm;
311: *obj = 0;
312: PetscObjectGetComm((PetscObject)info,&comm);
313: lambda = user->param;
314: hx = 1.0/(PetscReal)(info->mx-1);
315: hy = 1.0/(PetscReal)(info->my-1);
316: sc = hx*hy*lambda;
317: hxdhy = hx/hy;
318: hydhx = hy/hx;
319: /*
320: Compute function over the locally owned part of the grid
321: */
322: for (j=info->ys; j<info->ys+info->ym; j++) {
323: for (i=info->xs; i<info->xs+info->xm; i++) {
324: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
325: lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
326: } else {
327: u = x[j][i];
328: uw = x[j][i-1];
329: ue = x[j][i+1];
330: un = x[j-1][i];
331: us = x[j+1][i];
333: if (i-1 == 0) uw = 0.;
334: if (i+1 == info->mx-1) ue = 0.;
335: if (j-1 == 0) un = 0.;
336: if (j+1 == info->my-1) us = 0.;
338: /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
340: uxux = u*(2.*u - ue - uw)*hydhx;
341: uyuy = u*(2.*u - un - us)*hxdhy;
343: lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
344: }
345: }
346: }
347: PetscLogFlops(12.0*info->ym*info->xm);
348: MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
349: return(0);
350: }
354: /*
355: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
356: */
357: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jacpre,Mat jac,AppCtx *user)
358: {
360: PetscInt i,j,k;
361: MatStencil col[5],row;
362: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
365: lambda = user->param;
366: hx = 1.0/(PetscReal)(info->mx-1);
367: hy = 1.0/(PetscReal)(info->my-1);
368: sc = hx*hy*lambda;
369: hxdhy = hx/hy;
370: hydhx = hy/hx;
373: /*
374: Compute entries for the locally owned part of the Jacobian.
375: - Currently, all PETSc parallel matrix formats are partitioned by
376: contiguous chunks of rows across the processors.
377: - Each processor needs to insert only elements that it owns
378: locally (but any non-local elements will be sent to the
379: appropriate processor during matrix assembly).
380: - Here, we set all entries for a particular row at once.
381: - We can set matrix entries either using either
382: MatSetValuesLocal() or MatSetValues(), as discussed above.
383: */
384: for (j=info->ys; j<info->ys+info->ym; j++) {
385: for (i=info->xs; i<info->xs+info->xm; i++) {
386: row.j = j; row.i = i;
387: /* boundary points */
388: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
389: v[0] = 2.0*(hydhx + hxdhy);
390: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
391: } else {
392: k = 0;
393: /* interior grid points */
394: if (j-1 != 0) {
395: v[k] = -hxdhy;
396: col[k].j = j - 1; col[k].i = i;
397: k++;
398: }
399: if (i-1 != 0) {
400: v[k] = -hydhx;
401: col[k].j = j; col[k].i = i-1;
402: k++;
403: }
405: v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;
407: if (i+1 != info->mx-1) {
408: v[k] = -hydhx;
409: col[k].j = j; col[k].i = i+1;
410: k++;
411: }
412: if (j+1 != info->mx-1) {
413: v[k] = -hxdhy;
414: col[k].j = j + 1; col[k].i = i;
415: k++;
416: }
417: MatSetValuesStencil(jac,1,&row,k,col,v,INSERT_VALUES);
418: }
419: }
420: }
422: /*
423: Assemble matrix, using the 2-step process:
424: MatAssemblyBegin(), MatAssemblyEnd().
425: */
426: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
427: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
428: /*
429: Tell the matrix we will never add a new nonzero location to the
430: matrix. If we do, it will generate an error.
431: */
432: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
433: return(0);
434: }
436: #if defined(PETSC_HAVE_MATLAB_ENGINE)
439: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
440: {
441: AppCtx *user = (AppCtx*)ptr;
443: PetscInt Mx,My;
444: PetscReal lambda,hx,hy;
445: Vec localX,localF;
446: MPI_Comm comm;
447: DM da;
450: SNESGetDM(snes,&da);
451: DMGetLocalVector(da,&localX);
452: DMGetLocalVector(da,&localF);
453: PetscObjectSetName((PetscObject)localX,"localX");
454: PetscObjectSetName((PetscObject)localF,"localF");
455: 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);
457: lambda = user->param;
458: hx = 1.0/(PetscReal)(Mx-1);
459: hy = 1.0/(PetscReal)(My-1);
461: PetscObjectGetComm((PetscObject)snes,&comm);
462: /*
463: Scatter ghost points to local vector,using the 2-step process
464: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
465: By placing code between these two statements, computations can be
466: done while messages are in transition.
467: */
468: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
469: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
470: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
471: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
472: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
474: /*
475: Insert values into global vector
476: */
477: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
478: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
479: DMRestoreLocalVector(da,&localX);
480: DMRestoreLocalVector(da,&localF);
481: return(0);
482: }
483: #endif
485: /* ------------------------------------------------------------------- */
488: /*
489: Applies some sweeps on nonlinear Gauss-Seidel on each process
491: */
492: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
493: {
494: PetscInt i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
496: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
497: PetscScalar **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
498: PetscReal atol,rtol,stol;
499: DM da;
500: AppCtx *user;
501: Vec localX,localB;
504: tot_its = 0;
505: SNESNGSGetSweeps(snes,&sweeps);
506: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
507: SNESGetDM(snes,&da);
508: DMGetApplicationContext(da,(void**)&user);
510: 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);
512: lambda = user->param;
513: hx = 1.0/(PetscReal)(Mx-1);
514: hy = 1.0/(PetscReal)(My-1);
515: sc = hx*hy*lambda;
516: hxdhy = hx/hy;
517: hydhx = hy/hx;
520: DMGetLocalVector(da,&localX);
521: if (B) {
522: DMGetLocalVector(da,&localB);
523: }
524: for (l=0; l<sweeps; l++) {
526: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
527: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
528: if (B) {
529: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
530: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
531: }
532: /*
533: Get a pointer to vector data.
534: - For default PETSc vectors, VecGetArray() returns a pointer to
535: the data array. Otherwise, the routine is implementation dependent.
536: - You MUST call VecRestoreArray() when you no longer need access to
537: the array.
538: */
539: DMDAVecGetArray(da,localX,&x);
540: if (B) DMDAVecGetArray(da,localB,&b);
541: /*
542: Get local grid boundaries (for 2-dimensional DMDA):
543: xs, ys - starting grid indices (no ghost points)
544: xm, ym - widths of local grid (no ghost points)
545: */
546: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
548: for (j=ys; j<ys+ym; j++) {
549: for (i=xs; i<xs+xm; i++) {
550: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
551: /* boundary conditions are all zero Dirichlet */
552: x[j][i] = 0.0;
553: } else {
554: if (B) bij = b[j][i];
555: else bij = 0.;
557: u = x[j][i];
558: un = x[j-1][i];
559: us = x[j+1][i];
560: ue = x[j][i-1];
561: uw = x[j][i+1];
563: for (k=0; k<its; k++) {
564: eu = PetscExpScalar(u);
565: uxx = (2.0*u - ue - uw)*hydhx;
566: uyy = (2.0*u - un - us)*hxdhy;
567: F = uxx + uyy - sc*eu - bij;
568: if (k == 0) F0 = F;
569: J = 2.0*(hydhx + hxdhy) - sc*eu;
570: y = F/J;
571: u -= y;
572: tot_its++;
574: if (atol > PetscAbsReal(PetscRealPart(F)) ||
575: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
576: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
577: break;
578: }
579: }
580: x[j][i] = u;
581: }
582: }
583: }
584: /*
585: Restore vector
586: */
587: DMDAVecRestoreArray(da,localX,&x);
588: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
589: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
590: }
591: PetscLogFlops(tot_its*(21.0));
592: DMRestoreLocalVector(da,&localX);
593: if (B) {
594: DMDAVecRestoreArray(da,localB,&b);
595: DMRestoreLocalVector(da,&localB);
596: }
597: return(0);
598: }