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