Actual source code: ex5.c
petsc-3.6.4 2016-04-12
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
93: KSP ksp;
94: PetscInt lits,slits;
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Initialize program
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: PetscInitialize(&argc,&argv,(char*)0,help);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Initialize problem parameters
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: user.param = 6.0;
106: PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
107: 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);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create nonlinear solver context
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: SNESCreate(PETSC_COMM_WORLD,&snes);
113: SNESSetCountersReset(snes,PETSC_FALSE);
114: SNESSetNGS(snes, NonlinearGS, NULL);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create distributed array (DMDA) to manage parallel grid and vectors
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
120: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
121: DMSetApplicationContext(da,&user);
122: SNESSetDM(snes,da);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Extract global vectors from DMDA; then duplicate for remaining
125: vectors that are the same types
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: DMCreateGlobalVector(da,&x);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Set local function evaluation routine
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
133: PetscOptionsGetBool(NULL,"-fd",&flg,NULL);
134: if (!flg) {
135: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
136: }
138: PetscOptionsGetBool(NULL,"-obj",&flg,NULL);
139: if (flg) {
140: DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
141: }
143: #if defined(PETSC_HAVE_MATLAB_ENGINE)
144: PetscOptionsGetBool(NULL,"-matlab_function",&matlab_function,0);
145: if (matlab_function) {
146: VecDuplicate(x,&r);
147: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
148: }
149: #endif
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Customize nonlinear solver; set runtime options
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: SNESSetFromOptions(snes);
156: FormInitialGuess(da,&user,x);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Solve nonlinear system
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: SNESSolve(snes,NULL,x);
162: SNESGetIterationNumber(snes,&its);
164: SNESGetLinearSolveIterations(snes,&slits);
165: SNESGetKSP(snes,&ksp);
166: KSPGetTotalIterations(ksp,&lits);
167: if (lits != slits) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Number of total linear iterations reported by SNES %D does not match reported by KSP %D",slits,lits);
168: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Free work space. All PETSc objects should be destroyed when they
173: are no longer needed.
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: #if defined(PETSC_HAVE_MATLAB_ENGINE)
176: VecDestroy(&r);
177: #endif
178: VecDestroy(&x);
179: SNESDestroy(&snes);
180: DMDestroy(&da);
181: PetscFinalize();
182: return 0;
183: }
184: /* ------------------------------------------------------------------- */
187: /*
188: FormInitialGuess - Forms initial approximation.
190: Input Parameters:
191: user - user-defined application context
192: X - vector
194: Output Parameter:
195: X - vector
196: */
197: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
198: {
199: PetscInt i,j,Mx,My,xs,ys,xm,ym;
201: PetscReal lambda,temp1,temp,hx,hy;
202: PetscScalar **x;
205: 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);
207: lambda = user->param;
208: hx = 1.0/(PetscReal)(Mx-1);
209: hy = 1.0/(PetscReal)(My-1);
210: temp1 = lambda/(lambda + 1.0);
212: /*
213: Get a pointer to vector data.
214: - For default PETSc vectors, VecGetArray() returns a pointer to
215: the data array. Otherwise, the routine is implementation dependent.
216: - You MUST call VecRestoreArray() when you no longer need access to
217: the array.
218: */
219: DMDAVecGetArray(da,X,&x);
221: /*
222: Get local grid boundaries (for 2-dimensional DMDA):
223: xs, ys - starting grid indices (no ghost points)
224: xm, ym - widths of local grid (no ghost points)
226: */
227: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
229: /*
230: Compute initial guess over the locally owned part of the grid
231: */
232: for (j=ys; j<ys+ym; j++) {
233: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
234: for (i=xs; i<xs+xm; i++) {
235: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
236: /* boundary conditions are all zero Dirichlet */
237: x[j][i] = 0.0;
238: } else {
239: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
240: }
241: }
242: }
244: /*
245: Restore vector
246: */
247: DMDAVecRestoreArray(da,X,&x);
248: return(0);
249: }
250: /* ------------------------------------------------------------------- */
253: /*
254: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
257: */
258: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
259: {
261: PetscInt i,j;
262: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
263: PetscScalar u,ue,uw,un,us,uxx,uyy;
266: lambda = user->param;
267: hx = 1.0/(PetscReal)(info->mx-1);
268: hy = 1.0/(PetscReal)(info->my-1);
269: sc = hx*hy*lambda;
270: hxdhy = hx/hy;
271: hydhx = hy/hx;
272: /*
273: Compute function over the locally owned part of the grid
274: */
275: for (j=info->ys; j<info->ys+info->ym; j++) {
276: for (i=info->xs; i<info->xs+info->xm; i++) {
277: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
278: f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
279: } else {
280: u = x[j][i];
281: uw = x[j][i-1];
282: ue = x[j][i+1];
283: un = x[j-1][i];
284: us = x[j+1][i];
286: if (i-1 == 0) uw = 0.;
287: if (i+1 == info->mx-1) ue = 0.;
288: if (j-1 == 0) un = 0.;
289: if (j+1 == info->my-1) us = 0.;
291: uxx = (2.0*u - uw - ue)*hydhx;
292: uyy = (2.0*u - un - us)*hxdhy;
293: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
294: }
295: }
296: }
297: PetscLogFlops(11.0*info->ym*info->xm);
298: return(0);
299: }
304: /*
305: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
308: */
309: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
310: {
312: PetscInt i,j;
313: PetscReal lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
314: PetscScalar u,ue,uw,un,us,uxux,uyuy;
315: MPI_Comm comm;
318: *obj = 0;
319: PetscObjectGetComm((PetscObject)info,&comm);
320: lambda = user->param;
321: hx = 1.0/(PetscReal)(info->mx-1);
322: hy = 1.0/(PetscReal)(info->my-1);
323: sc = hx*hy*lambda;
324: hxdhy = hx/hy;
325: hydhx = hy/hx;
326: /*
327: Compute function over the locally owned part of the grid
328: */
329: for (j=info->ys; j<info->ys+info->ym; j++) {
330: for (i=info->xs; i<info->xs+info->xm; i++) {
331: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
332: lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
333: } else {
334: u = x[j][i];
335: uw = x[j][i-1];
336: ue = x[j][i+1];
337: un = x[j-1][i];
338: us = x[j+1][i];
340: if (i-1 == 0) uw = 0.;
341: if (i+1 == info->mx-1) ue = 0.;
342: if (j-1 == 0) un = 0.;
343: if (j+1 == info->my-1) us = 0.;
345: /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
347: uxux = u*(2.*u - ue - uw)*hydhx;
348: uyuy = u*(2.*u - un - us)*hxdhy;
350: lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
351: }
352: }
353: }
354: PetscLogFlops(12.0*info->ym*info->xm);
355: MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
356: return(0);
357: }
361: /*
362: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
363: */
364: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
365: {
367: PetscInt i,j,k;
368: MatStencil col[5],row;
369: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
372: lambda = user->param;
373: hx = 1.0/(PetscReal)(info->mx-1);
374: hy = 1.0/(PetscReal)(info->my-1);
375: sc = hx*hy*lambda;
376: hxdhy = hx/hy;
377: hydhx = hy/hx;
380: /*
381: Compute entries for the locally owned part of the Jacobian.
382: - Currently, all PETSc parallel matrix formats are partitioned by
383: contiguous chunks of rows across the processors.
384: - Each processor needs to insert only elements that it owns
385: locally (but any non-local elements will be sent to the
386: appropriate processor during matrix assembly).
387: - Here, we set all entries for a particular row at once.
388: - We can set matrix entries either using either
389: MatSetValuesLocal() or MatSetValues(), as discussed above.
390: */
391: for (j=info->ys; j<info->ys+info->ym; j++) {
392: for (i=info->xs; i<info->xs+info->xm; i++) {
393: row.j = j; row.i = i;
394: /* boundary points */
395: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
396: v[0] = 2.0*(hydhx + hxdhy);
397: MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
398: } else {
399: k = 0;
400: /* interior grid points */
401: if (j-1 != 0) {
402: v[k] = -hxdhy;
403: col[k].j = j - 1; col[k].i = i;
404: k++;
405: }
406: if (i-1 != 0) {
407: v[k] = -hydhx;
408: col[k].j = j; col[k].i = i-1;
409: k++;
410: }
412: v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;
414: if (i+1 != info->mx-1) {
415: v[k] = -hydhx;
416: col[k].j = j; col[k].i = i+1;
417: k++;
418: }
419: if (j+1 != info->mx-1) {
420: v[k] = -hxdhy;
421: col[k].j = j + 1; col[k].i = i;
422: k++;
423: }
424: MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
425: }
426: }
427: }
429: /*
430: Assemble matrix, using the 2-step process:
431: MatAssemblyBegin(), MatAssemblyEnd().
432: */
433: MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
434: MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
435: /*
436: Tell the matrix we will never add a new nonzero location to the
437: matrix. If we do, it will generate an error.
438: */
439: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
440: return(0);
441: }
443: #if defined(PETSC_HAVE_MATLAB_ENGINE)
446: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
447: {
448: AppCtx *user = (AppCtx*)ptr;
450: PetscInt Mx,My;
451: PetscReal lambda,hx,hy;
452: Vec localX,localF;
453: MPI_Comm comm;
454: DM da;
457: SNESGetDM(snes,&da);
458: DMGetLocalVector(da,&localX);
459: DMGetLocalVector(da,&localF);
460: PetscObjectSetName((PetscObject)localX,"localX");
461: PetscObjectSetName((PetscObject)localF,"localF");
462: 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);
464: lambda = user->param;
465: hx = 1.0/(PetscReal)(Mx-1);
466: hy = 1.0/(PetscReal)(My-1);
468: PetscObjectGetComm((PetscObject)snes,&comm);
469: /*
470: Scatter ghost points to local vector,using the 2-step process
471: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
472: By placing code between these two statements, computations can be
473: done while messages are in transition.
474: */
475: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
476: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
477: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
478: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
479: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
481: /*
482: Insert values into global vector
483: */
484: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
485: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
486: DMRestoreLocalVector(da,&localX);
487: DMRestoreLocalVector(da,&localF);
488: return(0);
489: }
490: #endif
492: /* ------------------------------------------------------------------- */
495: /*
496: Applies some sweeps on nonlinear Gauss-Seidel on each process
498: */
499: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
500: {
501: PetscInt i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
503: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
504: PetscScalar **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
505: PetscReal atol,rtol,stol;
506: DM da;
507: AppCtx *user;
508: Vec localX,localB;
511: tot_its = 0;
512: SNESNGSGetSweeps(snes,&sweeps);
513: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
514: SNESGetDM(snes,&da);
515: DMGetApplicationContext(da,(void**)&user);
517: 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);
519: lambda = user->param;
520: hx = 1.0/(PetscReal)(Mx-1);
521: hy = 1.0/(PetscReal)(My-1);
522: sc = hx*hy*lambda;
523: hxdhy = hx/hy;
524: hydhx = hy/hx;
527: DMGetLocalVector(da,&localX);
528: if (B) {
529: DMGetLocalVector(da,&localB);
530: }
531: for (l=0; l<sweeps; l++) {
533: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
534: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
535: if (B) {
536: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
537: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
538: }
539: /*
540: Get a pointer to vector data.
541: - For default PETSc vectors, VecGetArray() returns a pointer to
542: the data array. Otherwise, the routine is implementation dependent.
543: - You MUST call VecRestoreArray() when you no longer need access to
544: the array.
545: */
546: DMDAVecGetArray(da,localX,&x);
547: if (B) DMDAVecGetArray(da,localB,&b);
548: /*
549: Get local grid boundaries (for 2-dimensional DMDA):
550: xs, ys - starting grid indices (no ghost points)
551: xm, ym - widths of local grid (no ghost points)
552: */
553: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
555: for (j=ys; j<ys+ym; j++) {
556: for (i=xs; i<xs+xm; i++) {
557: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
558: /* boundary conditions are all zero Dirichlet */
559: x[j][i] = 0.0;
560: } else {
561: if (B) bij = b[j][i];
562: else bij = 0.;
564: u = x[j][i];
565: un = x[j-1][i];
566: us = x[j+1][i];
567: ue = x[j][i-1];
568: uw = x[j][i+1];
570: for (k=0; k<its; k++) {
571: eu = PetscExpScalar(u);
572: uxx = (2.0*u - ue - uw)*hydhx;
573: uyy = (2.0*u - un - us)*hxdhy;
574: F = uxx + uyy - sc*eu - bij;
575: if (k == 0) F0 = F;
576: J = 2.0*(hydhx + hxdhy) - sc*eu;
577: y = F/J;
578: u -= y;
579: tot_its++;
581: if (atol > PetscAbsReal(PetscRealPart(F)) ||
582: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
583: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
584: break;
585: }
586: }
587: x[j][i] = u;
588: }
589: }
590: }
591: /*
592: Restore vector
593: */
594: DMDAVecRestoreArray(da,localX,&x);
595: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
596: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
597: }
598: PetscLogFlops(tot_its*(21.0));
599: DMRestoreLocalVector(da,&localX);
600: if (B) {
601: DMDAVecRestoreArray(da,localB,&b);
602: DMRestoreLocalVector(da,&localB);
603: }
604: return(0);
605: }