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