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