Actual source code: ex1f.F90
1: !
2: ! Description: This example solves a nonlinear system on 1 processor with SNES.
3: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
4: ! domain. The command line options include:
5: ! -par <parameter>, where <parameter> indicates the nonlinearity of the problem
6: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
7: ! -mx <xg>, where <xg> = number of grid points in the x-direction
8: ! -my <yg>, where <yg> = number of grid points in the y-direction
9: !
11: !
12: ! --------------------------------------------------------------------------
13: !
14: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
15: ! the partial differential equation
16: !
17: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
18: !
19: ! with boundary conditions
20: !
21: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
22: !
23: ! A finite difference approximation with the usual 5-point stencil
24: ! is used to discretize the boundary value problem to obtain a nonlinear
25: ! system of equations.
26: !
27: ! The parallel version of this code is snes/tutorials/ex5f.F
28: !
29: ! --------------------------------------------------------------------------
30: subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr)
31: #include <petsc/finclude/petscsnes.h>
32: use petscsnes
33: implicit none
34: SNES snes
35: PetscReal norm
36: Vec tmp,x,y,w
37: PetscBool changed_w,changed_y
38: PetscErrorCode ierr
39: PetscInt ctx
40: PetscScalar mone
42: call VecDuplicate(x,tmp,ierr)
43: mone = -1.0
44: call VecWAXPY(tmp,mone,x,w,ierr)
45: call VecNorm(tmp,NORM_2,norm,ierr)
46: call VecDestroy(tmp,ierr)
47: print*, 'Norm of search step ',norm
48: changed_y = PETSC_FALSE
49: changed_w = PETSC_FALSE
50: return
51: end
53: program main
54: #include <petsc/finclude/petscdraw.h>
55: use petscsnes
56: implicit none
57: interface SNESSetJacobian
58: subroutine SNESSetJacobian1(a,b,c,d,e,z)
59: use petscsnes
60: SNES a
61: Mat b
62: Mat c
63: external d
64: MatFDColoring e
65: PetscErrorCode z
66: end subroutine
67: subroutine SNESSetJacobian2(a,b,c,d,e,z)
68: use petscsnes
69: SNES a
70: Mat b
71: Mat c
72: external d
73: integer e
74: PetscErrorCode z
75: end subroutine
76: end interface
77: !
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: ! Variable declarations
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: !
82: ! Variables:
83: ! snes - nonlinear solver
84: ! x, r - solution, residual vectors
85: ! J - Jacobian matrix
86: ! its - iterations for convergence
87: ! matrix_free - flag - 1 indicates matrix-free version
88: ! lambda - nonlinearity parameter
89: ! draw - drawing context
90: !
91: SNES snes
92: MatColoring mc
93: Vec x,r
94: PetscDraw draw
95: Mat J
96: PetscBool matrix_free,flg,fd_coloring
97: PetscErrorCode ierr
98: PetscInt its,N, mx,my,i5
99: PetscMPIInt size,rank
100: PetscReal lambda_max,lambda_min,lambda
101: MatFDColoring fdcoloring
102: ISColoring iscoloring
103: PetscBool pc
104: external postcheck
106: PetscScalar lx_v(0:1)
107: PetscOffset lx_i
109: ! Store parameters in common block
111: common /params/ lambda,mx,my,fd_coloring
113: ! Note: Any user-defined Fortran routines (such as FormJacobian)
114: ! MUST be declared as external.
116: external FormFunction,FormInitialGuess,FormJacobian
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: ! Initialize program
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
123: if (ierr .ne. 0) then
124: print*,'Unable to initialize PETSc'
125: stop
126: endif
127: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
128: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
130: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
132: ! Initialize problem parameters
133: i5 = 5
134: lambda_max = 6.81
135: lambda_min = 0.0
136: lambda = 6.0
137: mx = 4
138: my = 4
139: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
140: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
141: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
142: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
143: N = mx*my
144: pc = PETSC_FALSE
145: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr);
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: ! Create nonlinear solver context
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
153: if (pc .eqv. PETSC_TRUE) then
154: call SNESSetType(snes,SNESNEWTONTR,ierr)
155: call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
156: endif
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Create vector data structures; set function evaluation routine
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: call VecCreate(PETSC_COMM_WORLD,x,ierr)
163: call VecSetSizes(x,PETSC_DECIDE,N,ierr)
164: call VecSetFromOptions(x,ierr)
165: call VecDuplicate(x,r,ierr)
167: ! Set function evaluation routine and vector. Whenever the nonlinear
168: ! solver needs to evaluate the nonlinear function, it will call this
169: ! routine.
170: ! - Note that the final routine argument is the user-defined
171: ! context that provides application-specific data for the
172: ! function evaluation routine.
174: call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: ! Create matrix data structure; set Jacobian evaluation routine
178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: ! Create matrix. Here we only approximately preallocate storage space
181: ! for the Jacobian. See the users manual for a discussion of better
182: ! techniques for preallocating matrix memory.
184: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
185: if (.not. matrix_free) then
186: call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
187: endif
189: !
190: ! This option will cause the Jacobian to be computed via finite differences
191: ! efficiently using a coloring of the columns of the matrix.
192: !
193: fd_coloring = .false.
194: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
195: if (fd_coloring) then
197: !
198: ! This initializes the nonzero structure of the Jacobian. This is artificial
199: ! because clearly if we had a routine to compute the Jacobian we won't need
200: ! to use finite differences.
201: !
202: call FormJacobian(snes,x,J,J,0,ierr)
203: !
204: ! Color the matrix, i.e. determine groups of columns that share no common
205: ! rows. These columns in the Jacobian can all be computed simultaneously.
206: !
207: call MatColoringCreate(J,mc,ierr)
208: call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
209: call MatColoringSetFromOptions(mc,ierr)
210: call MatColoringApply(mc,iscoloring,ierr)
211: call MatColoringDestroy(mc,ierr)
212: !
213: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
214: ! to compute the actual Jacobians via finite differences.
215: !
216: call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
217: call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
218: call MatFDColoringSetFromOptions(fdcoloring,ierr)
219: call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
220: !
221: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
222: ! to compute Jacobians.
223: !
224: call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
225: call ISColoringDestroy(iscoloring,ierr)
227: else if (.not. matrix_free) then
229: ! Set Jacobian matrix data structure and default Jacobian evaluation
230: ! routine. Whenever the nonlinear solver needs to compute the
231: ! Jacobian matrix, it will call this routine.
232: ! - Note that the final routine argument is the user-defined
233: ! context that provides application-specific data for the
234: ! Jacobian evaluation routine.
235: ! - The user can override with:
236: ! -snes_fd : default finite differencing approximation of Jacobian
237: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
238: ! (unless user explicitly sets preconditioner)
239: ! -snes_mf_operator : form preconditioning matrix as set by the user,
240: ! but use matrix-free approx for Jacobian-vector
241: ! products within Newton-Krylov method
242: !
243: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
244: endif
246: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247: ! Customize nonlinear solver; set runtime options
248: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
252: call SNESSetFromOptions(snes,ierr)
254: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: ! Evaluate initial guess; then solve nonlinear system.
256: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: ! Note: The user should initialize the vector, x, with the initial guess
259: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
260: ! to employ an initial guess of zero, the user should explicitly set
261: ! this vector to zero by calling VecSet().
263: call FormInitialGuess(x,ierr)
264: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
265: call SNESGetIterationNumber(snes,its,ierr);
266: if (rank .eq. 0) then
267: write(6,100) its
268: endif
269: 100 format('Number of SNES iterations = ',i1)
271: ! PetscDraw contour plot of solution
273: call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
274: call PetscDrawSetFromOptions(draw,ierr)
276: call VecGetArrayRead(x,lx_v,lx_i,ierr)
277: call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
278: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: ! Free work space. All PETSc objects should be destroyed when they
282: ! are no longer needed.
283: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285: if (.not. matrix_free) call MatDestroy(J,ierr)
286: if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)
288: call VecDestroy(x,ierr)
289: call VecDestroy(r,ierr)
290: call SNESDestroy(snes,ierr)
291: call PetscDrawDestroy(draw,ierr)
292: call PetscFinalize(ierr)
293: end
295: ! ---------------------------------------------------------------------
296: !
297: ! FormInitialGuess - Forms initial approximation.
298: !
299: ! Input Parameter:
300: ! X - vector
301: !
302: ! Output Parameters:
303: ! X - vector
304: ! ierr - error code
305: !
306: ! Notes:
307: ! This routine serves as a wrapper for the lower-level routine
308: ! "ApplicationInitialGuess", where the actual computations are
309: ! done using the standard Fortran style of treating the local
310: ! vector data as a multidimensional array over the local mesh.
311: ! This routine merely accesses the local vector data via
312: ! VecGetArray() and VecRestoreArray().
313: !
314: subroutine FormInitialGuess(X,ierr)
315: use petscsnes
316: implicit none
318: ! Input/output variables:
319: Vec X
320: PetscErrorCode ierr
322: ! Declarations for use with local arrays:
323: PetscScalar lx_v(0:1)
324: PetscOffset lx_i
326: 0
328: ! Get a pointer to vector data.
329: ! - For default PETSc vectors, VecGetArray() returns a pointer to
330: ! the data array. Otherwise, the routine is implementation dependent.
331: ! - You MUST call VecRestoreArray() when you no longer need access to
332: ! the array.
333: ! - Note that the Fortran interface to VecGetArray() differs from the
334: ! C version. See the users manual for details.
336: call VecGetArray(X,lx_v,lx_i,ierr)
338: ! Compute initial guess
340: call ApplicationInitialGuess(lx_v(lx_i),ierr)
342: ! Restore vector
344: call VecRestoreArray(X,lx_v,lx_i,ierr)
346: return
347: end
349: ! ---------------------------------------------------------------------
350: !
351: ! ApplicationInitialGuess - Computes initial approximation, called by
352: ! the higher level routine FormInitialGuess().
353: !
354: ! Input Parameter:
355: ! x - local vector data
356: !
357: ! Output Parameters:
358: ! f - local vector data, f(x)
359: ! ierr - error code
360: !
361: ! Notes:
362: ! This routine uses standard Fortran-style computations over a 2-dim array.
363: !
364: subroutine ApplicationInitialGuess(x,ierr)
365: use petscksp
366: implicit none
368: ! Common blocks:
369: PetscReal lambda
370: PetscInt mx,my
371: PetscBool fd_coloring
372: common /params/ lambda,mx,my,fd_coloring
374: ! Input/output variables:
375: PetscScalar x(mx,my)
376: PetscErrorCode ierr
378: ! Local variables:
379: PetscInt i,j
380: PetscReal temp1,temp,hx,hy,one
382: ! Set parameters
384: 0
385: one = 1.0
386: hx = one/(mx-1)
387: hy = one/(my-1)
388: temp1 = lambda/(lambda + one)
390: do 20 j=1,my
391: temp = min(j-1,my-j)*hy
392: do 10 i=1,mx
393: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
394: x(i,j) = 0.0
395: else
396: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
397: endif
398: 10 continue
399: 20 continue
401: return
402: end
404: ! ---------------------------------------------------------------------
405: !
406: ! FormFunction - Evaluates nonlinear function, F(x).
407: !
408: ! Input Parameters:
409: ! snes - the SNES context
410: ! X - input vector
411: ! dummy - optional user-defined context, as set by SNESSetFunction()
412: ! (not used here)
413: !
414: ! Output Parameter:
415: ! F - vector with newly computed function
416: !
417: ! Notes:
418: ! This routine serves as a wrapper for the lower-level routine
419: ! "ApplicationFunction", where the actual computations are
420: ! done using the standard Fortran style of treating the local
421: ! vector data as a multidimensional array over the local mesh.
422: ! This routine merely accesses the local vector data via
423: ! VecGetArray() and VecRestoreArray().
424: !
425: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
426: use petscsnes
427: implicit none
429: ! Input/output variables:
430: SNES snes
431: Vec X,F
432: PetscErrorCode ierr
433: MatFDColoring fdcoloring
435: ! Common blocks:
436: PetscReal lambda
437: PetscInt mx,my
438: PetscBool fd_coloring
439: common /params/ lambda,mx,my,fd_coloring
441: ! Declarations for use with local arrays:
442: PetscScalar lx_v(0:1),lf_v(0:1)
443: PetscOffset lx_i,lf_i
445: PetscInt, pointer :: indices(:)
447: ! Get pointers to vector data.
448: ! - For default PETSc vectors, VecGetArray() returns a pointer to
449: ! the data array. Otherwise, the routine is implementation dependent.
450: ! - You MUST call VecRestoreArray() when you no longer need access to
451: ! the array.
452: ! - Note that the Fortran interface to VecGetArray() differs from the
453: ! C version. See the Fortran chapter of the users manual for details.
455: call VecGetArrayRead(X,lx_v,lx_i,ierr)
456: call VecGetArray(F,lf_v,lf_i,ierr)
458: ! Compute function
460: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)
462: ! Restore vectors
464: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
465: call VecRestoreArray(F,lf_v,lf_i,ierr)
467: call PetscLogFlops(11.0d0*mx*my,ierr)
468: !
469: ! fdcoloring is in the common block and used here ONLY to test the
470: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
471: !
472: if (fd_coloring) then
473: call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
474: print*,'Indices from GetPerturbedColumnsF90'
475: write(*,1000) indices
476: 1000 format(50i4)
477: call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
478: endif
479: return
480: end
482: ! ---------------------------------------------------------------------
483: !
484: ! ApplicationFunction - Computes nonlinear function, called by
485: ! the higher level routine FormFunction().
486: !
487: ! Input Parameter:
488: ! x - local vector data
489: !
490: ! Output Parameters:
491: ! f - local vector data, f(x)
492: ! ierr - error code
493: !
494: ! Notes:
495: ! This routine uses standard Fortran-style computations over a 2-dim array.
496: !
497: subroutine ApplicationFunction(x,f,ierr)
498: use petscsnes
499: implicit none
501: ! Common blocks:
502: PetscReal lambda
503: PetscInt mx,my
504: PetscBool fd_coloring
505: common /params/ lambda,mx,my,fd_coloring
507: ! Input/output variables:
508: PetscScalar x(mx,my),f(mx,my)
509: PetscErrorCode ierr
511: ! Local variables:
512: PetscScalar two,one,hx,hy
513: PetscScalar hxdhy,hydhx,sc
514: PetscScalar u,uxx,uyy
515: PetscInt i,j
517: 0
518: one = 1.0
519: two = 2.0
520: hx = one/(mx-1)
521: hy = one/(my-1)
522: sc = hx*hy*lambda
523: hxdhy = hx/hy
524: hydhx = hy/hx
526: ! Compute function
528: do 20 j=1,my
529: do 10 i=1,mx
530: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
531: f(i,j) = x(i,j)
532: else
533: u = x(i,j)
534: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
535: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
536: f(i,j) = uxx + uyy - sc*exp(u)
537: endif
538: 10 continue
539: 20 continue
541: return
542: end
544: ! ---------------------------------------------------------------------
545: !
546: ! FormJacobian - Evaluates Jacobian matrix.
547: !
548: ! Input Parameters:
549: ! snes - the SNES context
550: ! x - input vector
551: ! dummy - optional user-defined context, as set by SNESSetJacobian()
552: ! (not used here)
553: !
554: ! Output Parameters:
555: ! jac - Jacobian matrix
556: ! jac_prec - optionally different preconditioning matrix (not used here)
557: ! flag - flag indicating matrix structure
558: !
559: ! Notes:
560: ! This routine serves as a wrapper for the lower-level routine
561: ! "ApplicationJacobian", where the actual computations are
562: ! done using the standard Fortran style of treating the local
563: ! vector data as a multidimensional array over the local mesh.
564: ! This routine merely accesses the local vector data via
565: ! VecGetArray() and VecRestoreArray().
566: !
567: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
568: use petscsnes
569: implicit none
571: ! Input/output variables:
572: SNES snes
573: Vec X
574: Mat jac,jac_prec
575: PetscErrorCode ierr
576: integer dummy
578: ! Common blocks:
579: PetscReal lambda
580: PetscInt mx,my
581: PetscBool fd_coloring
582: common /params/ lambda,mx,my,fd_coloring
584: ! Declarations for use with local array:
585: PetscScalar lx_v(0:1)
586: PetscOffset lx_i
588: ! Get a pointer to vector data
590: call VecGetArrayRead(X,lx_v,lx_i,ierr)
592: ! Compute Jacobian entries
594: call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)
596: ! Restore vector
598: call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
600: ! Assemble matrix
602: call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
603: call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
605: return
606: end
608: ! ---------------------------------------------------------------------
609: !
610: ! ApplicationJacobian - Computes Jacobian matrix, called by
611: ! the higher level routine FormJacobian().
612: !
613: ! Input Parameters:
614: ! x - local vector data
615: !
616: ! Output Parameters:
617: ! jac - Jacobian matrix
618: ! jac_prec - optionally different preconditioning matrix (not used here)
619: ! ierr - error code
620: !
621: ! Notes:
622: ! This routine uses standard Fortran-style computations over a 2-dim array.
623: !
624: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
625: use petscsnes
626: implicit none
628: ! Common blocks:
629: PetscReal lambda
630: PetscInt mx,my
631: PetscBool fd_coloring
632: common /params/ lambda,mx,my,fd_coloring
634: ! Input/output variables:
635: PetscScalar x(mx,my)
636: Mat jac,jac_prec
637: PetscErrorCode ierr
639: ! Local variables:
640: PetscInt i,j,row(1),col(5),i1,i5
641: PetscScalar two,one, hx,hy
642: PetscScalar hxdhy,hydhx,sc,v(5)
644: ! Set parameters
646: i1 = 1
647: i5 = 5
648: one = 1.0
649: two = 2.0
650: hx = one/(mx-1)
651: hy = one/(my-1)
652: sc = hx*hy
653: hxdhy = hx/hy
654: hydhx = hy/hx
656: ! Compute entries of the Jacobian matrix
657: ! - Here, we set all entries for a particular row at once.
658: ! - Note that MatSetValues() uses 0-based row and column numbers
659: ! in Fortran as well as in C.
661: do 20 j=1,my
662: row(1) = (j-1)*mx - 1
663: do 10 i=1,mx
664: row(1) = row(1) + 1
665: ! boundary points
666: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
667: call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
668: ! interior grid points
669: else
670: v(1) = -hxdhy
671: v(2) = -hydhx
672: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
673: v(4) = -hydhx
674: v(5) = -hxdhy
675: col(1) = row(1) - mx
676: col(2) = row(1) - 1
677: col(3) = row(1)
678: col(4) = row(1) + 1
679: col(5) = row(1) + mx
680: call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
681: endif
682: 10 continue
683: 20 continue
685: return
686: end
688: !
689: !/*TEST
690: !
691: ! build:
692: ! requires: !single
693: !
694: ! test:
695: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
696: !
697: ! test:
698: ! suffix: 2
699: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
700: !
701: ! test:
702: ! suffix: 3
703: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
704: ! filter: sort -b
705: ! filter_output: sort -b
706: !
707: ! test:
708: ! suffix: 4
709: ! args: -pc -par 6.807 -nox
710: !
711: !TEST*/