xref: /petsc/src/snes/tests/ex1f.F90 (revision daa037dfd3c3bec8dc8659548d2b20b07c1dc6de)
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
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
41
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
52
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
105
106      PetscScalar        lx_v(0:1)
107      PetscOffset        lx_i
108
109!  Store parameters in common block
110
111      common /params/ lambda,mx,my,fd_coloring
112
113!  Note: Any user-defined Fortran routines (such as FormJacobian)
114!  MUST be declared as external.
115
116      external FormFunction,FormInitialGuess,FormJacobian
117
118! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119!  Initialize program
120! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121
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)
129
130      if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
131
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);
146
147! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148!  Create nonlinear solver context
149! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150
151      call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
152
153      if (pc .eqv. PETSC_TRUE) then
154        call SNESSetType(snes,SNESNEWTONTR,ierr)
155        call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
156      endif
157
158! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159!  Create vector data structures; set function evaluation routine
160! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161
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)
166
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.
173
174      call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)
175
176! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177!  Create matrix data structure; set Jacobian evaluation routine
178! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179
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.
183
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
188
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
196
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)
226
227      else if (.not. matrix_free) then
228
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
245
246! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247!  Customize nonlinear solver; set runtime options
248! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249
250!  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
251
252      call SNESSetFromOptions(snes,ierr)
253
254! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255!  Evaluate initial guess; then solve nonlinear system.
256! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257
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().
262
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)
270
271!  PetscDraw contour plot of solution
272
273      call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
274      call PetscDrawSetFromOptions(draw,ierr)
275
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)
279
280! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281!  Free work space.  All PETSc objects should be destroyed when they
282!  are no longer needed.
283! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284
285      if (.not. matrix_free) call MatDestroy(J,ierr)
286      if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)
287
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
294
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
317
318!  Input/output variables:
319      Vec           X
320      PetscErrorCode    ierr
321
322!  Declarations for use with local arrays:
323      PetscScalar   lx_v(0:1)
324      PetscOffset   lx_i
325
326      ierr   = 0
327
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.
335
336      call VecGetArray(X,lx_v,lx_i,ierr)
337
338!  Compute initial guess
339
340      call ApplicationInitialGuess(lx_v(lx_i),ierr)
341
342!  Restore vector
343
344      call VecRestoreArray(X,lx_v,lx_i,ierr)
345
346      return
347      end
348
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
367
368!  Common blocks:
369      PetscReal   lambda
370      PetscInt     mx,my
371      PetscBool         fd_coloring
372      common      /params/ lambda,mx,my,fd_coloring
373
374!  Input/output variables:
375      PetscScalar x(mx,my)
376      PetscErrorCode     ierr
377
378!  Local variables:
379      PetscInt     i,j
380      PetscReal temp1,temp,hx,hy,one
381
382!  Set parameters
383
384      ierr   = 0
385      one    = 1.0
386      hx     = one/(mx-1)
387      hy     = one/(my-1)
388      temp1  = lambda/(lambda + one)
389
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
400
401      return
402      end
403
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
428
429!  Input/output variables:
430      SNES              snes
431      Vec               X,F
432      PetscErrorCode          ierr
433      MatFDColoring fdcoloring
434
435!  Common blocks:
436      PetscReal         lambda
437      PetscInt          mx,my
438      PetscBool         fd_coloring
439      common            /params/ lambda,mx,my,fd_coloring
440
441!  Declarations for use with local arrays:
442      PetscScalar       lx_v(0:1),lf_v(0:1)
443      PetscOffset       lx_i,lf_i
444
445      PetscInt, pointer :: indices(:)
446
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.
454
455      call VecGetArrayRead(X,lx_v,lx_i,ierr)
456      call VecGetArray(F,lf_v,lf_i,ierr)
457
458!  Compute function
459
460      call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)
461
462!  Restore vectors
463
464      call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
465      call VecRestoreArray(F,lf_v,lf_i,ierr)
466
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
481
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
500
501!  Common blocks:
502      PetscReal      lambda
503      PetscInt        mx,my
504      PetscBool         fd_coloring
505      common         /params/ lambda,mx,my,fd_coloring
506
507!  Input/output variables:
508      PetscScalar    x(mx,my),f(mx,my)
509      PetscErrorCode       ierr
510
511!  Local variables:
512      PetscScalar    two,one,hx,hy
513      PetscScalar    hxdhy,hydhx,sc
514      PetscScalar    u,uxx,uyy
515      PetscInt        i,j
516
517      ierr   = 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
525
526!  Compute function
527
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
540
541      return
542      end
543
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
570
571!  Input/output variables:
572      SNES          snes
573      Vec           X
574      Mat           jac,jac_prec
575      PetscErrorCode      ierr
576      integer dummy
577
578!  Common blocks:
579      PetscReal     lambda
580      PetscInt       mx,my
581      PetscBool         fd_coloring
582      common        /params/ lambda,mx,my,fd_coloring
583
584!  Declarations for use with local array:
585      PetscScalar   lx_v(0:1)
586      PetscOffset   lx_i
587
588!  Get a pointer to vector data
589
590      call VecGetArrayRead(X,lx_v,lx_i,ierr)
591
592!  Compute Jacobian entries
593
594      call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)
595
596!  Restore vector
597
598      call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
599
600!  Assemble matrix
601
602      call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
603      call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
604
605      return
606      end
607
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
627
628!  Common blocks:
629      PetscReal    lambda
630      PetscInt      mx,my
631      PetscBool         fd_coloring
632      common       /params/ lambda,mx,my,fd_coloring
633
634!  Input/output variables:
635      PetscScalar  x(mx,my)
636      Mat          jac,jac_prec
637      PetscErrorCode      ierr
638
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)
643
644!  Set parameters
645
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
655
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.
660
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
684
685      return
686      end
687
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*/
712