xref: /petsc/src/snes/tests/ex1f.F90 (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce)
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!  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
45
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
56
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
109
110      PetscScalar        lx_v(0:1)
111      PetscOffset        lx_i
112
113!  Store parameters in common block
114
115      common /params/ lambda,mx,my,fd_coloring
116
117!  Note: Any user-defined Fortran routines (such as FormJacobian)
118!  MUST be declared as external.
119
120      external FormFunction,FormInitialGuess,FormJacobian
121
122! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123!  Initialize program
124! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125
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)
133
134      if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
135
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);
150
151! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152!  Create nonlinear solver context
153! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154
155      call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
156
157      if (pc .eqv. PETSC_TRUE) then
158        call SNESSetType(snes,SNESNEWTONTR,ierr)
159        call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
160      endif
161
162! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163!  Create vector data structures; set function evaluation routine
164! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165
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)
170
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.
177
178      call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)
179
180! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181!  Create matrix data structure; set Jacobian evaluation routine
182! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183
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.
187
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
192
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
200
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 simulataneously.
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)
230
231      else if (.not. matrix_free) then
232
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
249
250! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251!  Customize nonlinear solver; set runtime options
252! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253
254!  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
255
256      call SNESSetFromOptions(snes,ierr)
257
258! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259!  Evaluate initial guess; then solve nonlinear system.
260! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
261
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().
266
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)
274
275!  PetscDraw contour plot of solution
276
277      call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
278      call PetscDrawSetFromOptions(draw,ierr)
279
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)
283
284! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285!  Free work space.  All PETSc objects should be destroyed when they
286!  are no longer needed.
287! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288
289      if (.not. matrix_free) call MatDestroy(J,ierr)
290      if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)
291
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
298
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
321
322!  Input/output variables:
323      Vec           X
324      PetscErrorCode    ierr
325
326!  Declarations for use with local arrays:
327      PetscScalar   lx_v(0:1)
328      PetscOffset   lx_i
329
330      ierr   = 0
331
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.
339
340      call VecGetArray(X,lx_v,lx_i,ierr)
341
342!  Compute initial guess
343
344      call ApplicationInitialGuess(lx_v(lx_i),ierr)
345
346!  Restore vector
347
348      call VecRestoreArray(X,lx_v,lx_i,ierr)
349
350      return
351      end
352
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
371
372!  Common blocks:
373      PetscReal   lambda
374      PetscInt     mx,my
375      PetscBool         fd_coloring
376      common      /params/ lambda,mx,my,fd_coloring
377
378!  Input/output variables:
379      PetscScalar x(mx,my)
380      PetscErrorCode     ierr
381
382!  Local variables:
383      PetscInt     i,j
384      PetscReal temp1,temp,hx,hy,one
385
386!  Set parameters
387
388      ierr   = 0
389      one    = 1.0
390      hx     = one/(mx-1)
391      hy     = one/(my-1)
392      temp1  = lambda/(lambda + one)
393
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
404
405      return
406      end
407
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
432
433!  Input/output variables:
434      SNES              snes
435      Vec               X,F
436      PetscErrorCode          ierr
437      MatFDColoring fdcoloring
438
439!  Common blocks:
440      PetscReal         lambda
441      PetscInt          mx,my
442      PetscBool         fd_coloring
443      common            /params/ lambda,mx,my,fd_coloring
444
445!  Declarations for use with local arrays:
446      PetscScalar       lx_v(0:1),lf_v(0:1)
447      PetscOffset       lx_i,lf_i
448
449      PetscInt, pointer :: indices(:)
450
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.
458
459      call VecGetArrayRead(X,lx_v,lx_i,ierr)
460      call VecGetArray(F,lf_v,lf_i,ierr)
461
462!  Compute function
463
464      call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)
465
466!  Restore vectors
467
468      call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
469      call VecRestoreArray(F,lf_v,lf_i,ierr)
470
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
485
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
504
505!  Common blocks:
506      PetscReal      lambda
507      PetscInt        mx,my
508      PetscBool         fd_coloring
509      common         /params/ lambda,mx,my,fd_coloring
510
511!  Input/output variables:
512      PetscScalar    x(mx,my),f(mx,my)
513      PetscErrorCode       ierr
514
515!  Local variables:
516      PetscScalar    two,one,hx,hy
517      PetscScalar    hxdhy,hydhx,sc
518      PetscScalar    u,uxx,uyy
519      PetscInt        i,j
520
521      ierr   = 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
529
530!  Compute function
531
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
544
545      return
546      end
547
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
574
575!  Input/output variables:
576      SNES          snes
577      Vec           X
578      Mat           jac,jac_prec
579      PetscErrorCode      ierr
580      integer dummy
581
582!  Common blocks:
583      PetscReal     lambda
584      PetscInt       mx,my
585      PetscBool         fd_coloring
586      common        /params/ lambda,mx,my,fd_coloring
587
588!  Declarations for use with local array:
589      PetscScalar   lx_v(0:1)
590      PetscOffset   lx_i
591
592!  Get a pointer to vector data
593
594      call VecGetArrayRead(X,lx_v,lx_i,ierr)
595
596!  Compute Jacobian entries
597
598      call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)
599
600!  Restore vector
601
602      call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
603
604!  Assemble matrix
605
606      call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
607      call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
608
609      return
610      end
611
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
631
632!  Common blocks:
633      PetscReal    lambda
634      PetscInt      mx,my
635      PetscBool         fd_coloring
636      common       /params/ lambda,mx,my,fd_coloring
637
638!  Input/output variables:
639      PetscScalar  x(mx,my)
640      Mat          jac,jac_prec
641      PetscErrorCode      ierr
642
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)
647
648!  Set parameters
649
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
659
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.
664
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
688
689      return
690      end
691
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*/
716