Actual source code: ex1f.F90

  1: !
  2: !  Description: This example solves a nonlinear system on 1 processor with SNES.
  3: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
  4: !  domain.  The command line options include:
  5: !    -par <parameter>, where <parameter> indicates the nonlinearity of the problem
  6: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
  7: !    -mx <xg>, where <xg> = number of grid points in the x-direction
  8: !    -my <yg>, where <yg> = number of grid points in the y-direction
  9: !

 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:       MPI_Comm       comm

 43:       character(len=PETSC_MAX_PATH_LEN) :: outputString

 45:       PetscCallA(PetscObjectGetComm(snes,comm,ierr))
 46:       PetscCallA(VecDuplicate(x,tmp,ierr))
 47:       mone = -1.0
 48:       PetscCallA(VecWAXPY(tmp,mone,x,w,ierr))
 49:       PetscCallA(VecNorm(tmp,NORM_2,norm,ierr))
 50:       PetscCallA(VecDestroy(tmp,ierr))
 51:       write(outputString,*) norm
 52:       PetscCallA(PetscPrintf(comm,'Norm of search step '//trim(outputString)//'\n',ierr))
 53:       end

 55:       program main
 56: #include <petsc/finclude/petscdraw.h>
 57:       use petscsnes
 58:       implicit none
 59: !
 60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61: !                   Variable declarations
 62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63: !
 64: !  Variables:
 65: !     snes        - nonlinear solver
 66: !     x, r        - solution, residual vectors
 67: !     J           - Jacobian matrix
 68: !     its         - iterations for convergence
 69: !     matrix_free - flag - 1 indicates matrix-free version
 70: !     lambda      - nonlinearity parameter
 71: !     draw        - drawing context
 72: !
 73:       SNES               snes
 74:       MatColoring        mc
 75:       Vec                x,r
 76:       PetscDraw               draw
 77:       Mat                J
 78:       PetscBool  matrix_free,flg,fd_coloring
 79:       PetscErrorCode ierr
 80:       PetscInt   its,N, mx,my,i5
 81:       PetscMPIInt size,rank
 82:       PetscReal   lambda_max,lambda_min,lambda
 83:       MatFDColoring      fdcoloring
 84:       ISColoring         iscoloring
 85:       PetscBool          pc
 86:       external           postcheck

 88:       character(len=PETSC_MAX_PATH_LEN) :: outputString

 90:       PetscScalar,pointer :: lx_v(:)

 92: !  Store parameters in common block

 94:       common /params/ lambda,mx,my,fd_coloring

 96: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 97: !  MUST be declared as external.

 99:       external FormFunction,FormInitialGuess,FormJacobian

101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: !  Initialize program
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

105:       PetscCallA(PetscInitialize(ierr))
106:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
107:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))

109:       PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')

111: !  Initialize problem parameters
112:       i5 = 5
113:       lambda_max = 6.81
114:       lambda_min = 0.0
115:       lambda     = 6.0
116:       mx         = 4
117:       my         = 4
118:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
119:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
120:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
121:       PetscCheckA(lambda .lt. lambda_max .and. lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range ')
122:       N  = mx*my
123:       pc = PETSC_FALSE
124:       PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))

126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: !  Create nonlinear solver context
128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

130:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

132:       if (pc .eqv. PETSC_TRUE) then
133:         PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
134:         PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
135:       endif

137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: !  Create vector data structures; set function evaluation routine
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

141:       PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
142:       PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
143:       PetscCallA(VecSetFromOptions(x,ierr))
144:       PetscCallA(VecDuplicate(x,r,ierr))

146: !  Set function evaluation routine and vector.  Whenever the nonlinear
147: !  solver needs to evaluate the nonlinear function, it will call this
148: !  routine.
149: !   - Note that the final routine argument is the user-defined
150: !     context that provides application-specific data for the
151: !     function evaluation routine.

153:       PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))

155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: !  Create matrix data structure; set Jacobian evaluation routine
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

159: !  Create matrix. Here we only approximately preallocate storage space
160: !  for the Jacobian.  See the users manual for a discussion of better
161: !  techniques for preallocating matrix memory.

163:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
164:       if (.not. matrix_free) then
165:         PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER_ARRAY,J,ierr))
166:       endif

168: !
169: !     This option will cause the Jacobian to be computed via finite differences
170: !    efficiently using a coloring of the columns of the matrix.
171: !
172:       fd_coloring = .false.
173:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
174:       if (fd_coloring) then

176: !
177: !      This initializes the nonzero structure of the Jacobian. This is artificial
178: !      because clearly if we had a routine to compute the Jacobian we won't need
179: !      to use finite differences.
180: !
181:         PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
182: !
183: !       Color the matrix, i.e. determine groups of columns that share no common
184: !      rows. These columns in the Jacobian can all be computed simultaneously.
185: !
186:         PetscCallA(MatColoringCreate(J,mc,ierr))
187:         PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
188:         PetscCallA(MatColoringSetFromOptions(mc,ierr))
189:         PetscCallA(MatColoringApply(mc,iscoloring,ierr))
190:         PetscCallA(MatColoringDestroy(mc,ierr))
191: !
192: !       Create the data structure that SNESComputeJacobianDefaultColor() uses
193: !       to compute the actual Jacobians via finite differences.
194: !
195:         PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
196:         PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
197:         PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
198:         PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
199: !
200: !        Tell SNES to use the routine SNESComputeJacobianDefaultColor()
201: !      to compute Jacobians.
202: !
203:         PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
204:         PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr))
205:         PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,fdcoloring,ierr))
206:         PetscCallA(ISColoringDestroy(iscoloring,ierr))

208:       else if (.not. matrix_free) then

210: !  Set Jacobian matrix data structure and default Jacobian evaluation
211: !  routine.  Whenever the nonlinear solver needs to compute the
212: !  Jacobian matrix, it will call this routine.
213: !   - Note that the final routine argument is the user-defined
214: !     context that provides application-specific data for the
215: !     Jacobian evaluation routine.
216: !   - The user can override with:
217: !      -snes_fd : default finite differencing approximation of Jacobian
218: !      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
219: !                 (unless user explicitly sets preconditioner)
220: !      -snes_mf_operator : form preconditioning matrix as set by the user,
221: !                          but use matrix-free approx for Jacobian-vector
222: !                          products within Newton-Krylov method
223: !
224:         PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
225:       endif

227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: !  Customize nonlinear solver; set runtime options
229: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

231: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)

233:       PetscCallA(SNESSetFromOptions(snes,ierr))

235: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: !  Evaluate initial guess; then solve nonlinear system.
237: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

239: !  Note: The user should initialize the vector, x, with the initial guess
240: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
241: !  to employ an initial guess of zero, the user should explicitly set
242: !  this vector to zero by calling VecSet().

244:       PetscCallA(FormInitialGuess(x,ierr))
245:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
246:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))
247:       write(outputString,*) its
248:       PetscCallA(PetscPrintf(PETSC_COMM_WORLD,'Number of SNES iterations = '//trim(outputString)//'\n',ierr))

250: !  PetscDraw contour plot of solution

252:       PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
253:       PetscCallA(PetscDrawSetFromOptions(draw,ierr))

255:       PetscCallA(VecGetArrayReadF90(x,lx_v,ierr))
256:       PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v,ierr))
257:       PetscCallA(VecRestoreArrayReadF90(x,lx_v,ierr))

259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: !  Free work space.  All PETSc objects should be destroyed when they
261: !  are no longer needed.
262: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

264:       if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
265:       if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))

267:       PetscCallA(VecDestroy(x,ierr))
268:       PetscCallA(VecDestroy(r,ierr))
269:       PetscCallA(SNESDestroy(snes,ierr))
270:       PetscCallA(PetscDrawDestroy(draw,ierr))
271:       PetscCallA(PetscFinalize(ierr))
272:       end

274: ! ---------------------------------------------------------------------
275: !
276: !  FormInitialGuess - Forms initial approximation.
277: !
278: !  Input Parameter:
279: !  X - vector
280: !
281: !  Output Parameters:
282: !  X - vector
283: !  ierr - error code
284: !
285: !  Notes:
286: !  This routine serves as a wrapper for the lower-level routine
287: !  "ApplicationInitialGuess", where the actual computations are
288: !  done using the standard Fortran style of treating the local
289: !  vector data as a multidimensional array over the local mesh.
290: !  This routine merely accesses the local vector data via
291: !  VecGetArrayF90() and VecRestoreArrayF90().
292: !
293:       subroutine FormInitialGuess(X,ierr)
294:       use petscsnes
295:       implicit none

297: !  Input/output variables:
298:       Vec           X
299:       PetscErrorCode    ierr

301: !     Declarations for use with local arrays:
302:       PetscScalar,pointer :: lx_v(:)

304:       ierr   = 0

306: !  Get a pointer to vector data.
307: !    - VecGetArrayF90() returns a pointer to the data array.
308: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
309: !      the array.

311:       PetscCallA(VecGetArrayF90(X,lx_v,ierr))

313: !  Compute initial guess

315:       PetscCallA(ApplicationInitialGuess(lx_v,ierr))

317: !  Restore vector

319:       PetscCallA(VecRestoreArrayF90(X,lx_v,ierr))

321:       end

323: !  ApplicationInitialGuess - Computes initial approximation, called by
324: !  the higher level routine FormInitialGuess().
325: !
326: !  Input Parameter:
327: !  x - local vector data
328: !
329: !  Output Parameters:
330: !  f - local vector data, f(x)
331: !  ierr - error code
332: !
333: !  Notes:
334: !  This routine uses standard Fortran-style computations over a 2-dim array.
335: !
336:       subroutine ApplicationInitialGuess(x,ierr)
337:       use petscksp
338:       implicit none

340: !  Common blocks:
341:       PetscReal   lambda
342:       PetscInt     mx,my
343:       PetscBool         fd_coloring
344:       common      /params/ lambda,mx,my,fd_coloring

346: !  Input/output variables:
347:       PetscScalar x(mx,my)
348:       PetscErrorCode     ierr

350: !  Local variables:
351:       PetscInt     i,j
352:       PetscReal temp1,temp,hx,hy,one

354: !  Set parameters

356:       ierr   = 0
357:       one    = 1.0
358:       hx     = one/(mx-1)
359:       hy     = one/(my-1)
360:       temp1  = lambda/(lambda + one)

362:       do 20 j=1,my
363:          temp = min(j-1,my-j)*hy
364:          do 10 i=1,mx
365:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
366:               x(i,j) = 0.0
367:             else
368:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
369:             endif
370:  10      continue
371:  20   continue

373:       end

375: ! ---------------------------------------------------------------------
376: !
377: !  FormFunction - Evaluates nonlinear function, F(x).
378: !
379: !  Input Parameters:
380: !  snes  - the SNES context
381: !  X     - input vector
382: !  dummy - optional user-defined context, as set by SNESSetFunction()
383: !          (not used here)
384: !
385: !  Output Parameter:
386: !  F     - vector with newly computed function
387: !
388: !  Notes:
389: !  This routine serves as a wrapper for the lower-level routine
390: !  "ApplicationFunction", where the actual computations are
391: !  done using the standard Fortran style of treating the local
392: !  vector data as a multidimensional array over the local mesh.
393: !  This routine merely accesses the local vector data via
394: !  VecGetArrayF90() and VecRestoreArrayF90().
395: !
396:       subroutine FormFunction(snes,X,F,fdcoloring,ierr)
397:       use petscsnes
398:       implicit none

400: !  Input/output variables:
401:       SNES              snes
402:       Vec               X,F
403:       PetscErrorCode          ierr
404:       MatFDColoring fdcoloring

406: !  Common blocks:
407:       PetscReal         lambda
408:       PetscInt          mx,my
409:       PetscBool         fd_coloring
410:       common            /params/ lambda,mx,my,fd_coloring

412: !  Declarations for use with local arrays:
413:       PetscScalar,pointer :: lx_v(:), lf_v(:)
414:       PetscInt, pointer :: indices(:)

416: !  Get pointers to vector data.
417: !    - VecGetArrayF90() returns a pointer to the data array.
418: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
419: !      the array.

421:       PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
422:       PetscCallA(VecGetArrayF90(F,lf_v,ierr))

424: !  Compute function

426:       PetscCallA(ApplicationFunction(lx_v,lf_v,ierr))

428: !  Restore vectors

430:       PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
431:       PetscCallA(VecRestoreArrayF90(F,lf_v,ierr))

433:       PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
434: !
435: !     fdcoloring is in the common block and used here ONLY to test the
436: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
437: !
438:       if (fd_coloring) then
439:          PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
440:          print*,'Indices from GetPerturbedColumnsF90'
441:          write(*,1000) indices
442:  1000    format(50i4)
443:          PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
444:       endif
445:       end

447: ! ---------------------------------------------------------------------
448: !
449: !  ApplicationFunction - Computes nonlinear function, called by
450: !  the higher level routine FormFunction().
451: !
452: !  Input Parameter:
453: !  x    - local vector data
454: !
455: !  Output Parameters:
456: !  f    - local vector data, f(x)
457: !  ierr - error code
458: !
459: !  Notes:
460: !  This routine uses standard Fortran-style computations over a 2-dim array.
461: !
462:       subroutine ApplicationFunction(x,f,ierr)
463:       use petscsnes
464:       implicit none

466: !  Common blocks:
467:       PetscReal      lambda
468:       PetscInt        mx,my
469:       PetscBool         fd_coloring
470:       common         /params/ lambda,mx,my,fd_coloring

472: !  Input/output variables:
473:       PetscScalar    x(mx,my),f(mx,my)
474:       PetscErrorCode       ierr

476: !  Local variables:
477:       PetscScalar    two,one,hx,hy
478:       PetscScalar    hxdhy,hydhx,sc
479:       PetscScalar    u,uxx,uyy
480:       PetscInt        i,j

482:       ierr   = 0
483:       one    = 1.0
484:       two    = 2.0
485:       hx     = one/(mx-1)
486:       hy     = one/(my-1)
487:       sc     = hx*hy*lambda
488:       hxdhy  = hx/hy
489:       hydhx  = hy/hx

491: !  Compute function

493:       do 20 j=1,my
494:          do 10 i=1,mx
495:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
496:                f(i,j) = x(i,j)
497:             else
498:                u = x(i,j)
499:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
500:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
501:                f(i,j) = uxx + uyy - sc*exp(u)
502:             endif
503:  10      continue
504:  20   continue

506:       end

508: ! ---------------------------------------------------------------------
509: !
510: !  FormJacobian - Evaluates Jacobian matrix.
511: !
512: !  Input Parameters:
513: !  snes    - the SNES context
514: !  x       - input vector
515: !  dummy   - optional user-defined context, as set by SNESSetJacobian()
516: !            (not used here)
517: !
518: !  Output Parameters:
519: !  jac      - Jacobian matrix
520: !  jac_prec - optionally different preconditioning matrix (not used here)
521: !
522: !  Notes:
523: !  This routine serves as a wrapper for the lower-level routine
524: !  "ApplicationJacobian", where the actual computations are
525: !  done using the standard Fortran style of treating the local
526: !  vector data as a multidimensional array over the local mesh.
527: !  This routine merely accesses the local vector data via
528: !  VecGetArrayF90() and VecRestoreArrayF90().
529: !
530:       subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
531:       use petscsnes
532:       implicit none

534: !  Input/output variables:
535:       SNES          snes
536:       Vec           X
537:       Mat           jac,jac_prec
538:       PetscErrorCode      ierr
539:       integer dummy

541: !  Common blocks:
542:       PetscReal     lambda
543:       PetscInt       mx,my
544:       PetscBool         fd_coloring
545:       common        /params/ lambda,mx,my,fd_coloring

547: !  Declarations for use with local array:
548:       PetscScalar,pointer :: lx_v(:)

550: !  Get a pointer to vector data

552:       PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))

554: !  Compute Jacobian entries

556:       PetscCallA(ApplicationJacobian(lx_v,jac,jac_prec,ierr))

558: !  Restore vector

560:       PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))

562: !  Assemble matrix

564:       PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
565:       PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))

567:       end

569: ! ---------------------------------------------------------------------
570: !
571: !  ApplicationJacobian - Computes Jacobian matrix, called by
572: !  the higher level routine FormJacobian().
573: !
574: !  Input Parameters:
575: !  x        - local vector data
576: !
577: !  Output Parameters:
578: !  jac      - Jacobian matrix
579: !  jac_prec - optionally different preconditioning matrix (not used here)
580: !  ierr     - error code
581: !
582: !  Notes:
583: !  This routine uses standard Fortran-style computations over a 2-dim array.
584: !
585:       subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
586:       use petscsnes
587:       implicit none

589: !  Common blocks:
590:       PetscReal    lambda
591:       PetscInt      mx,my
592:       PetscBool         fd_coloring
593:       common       /params/ lambda,mx,my,fd_coloring

595: !  Input/output variables:
596:       PetscScalar  x(mx,my)
597:       Mat          jac,jac_prec
598:       PetscErrorCode      ierr

600: !  Local variables:
601:       PetscInt      i,j,row(1),col(5),i1,i5
602:       PetscScalar  two,one, hx,hy
603:       PetscScalar  hxdhy,hydhx,sc,v(5)

605: !  Set parameters

607:       i1     = 1
608:       i5     = 5
609:       one    = 1.0
610:       two    = 2.0
611:       hx     = one/(mx-1)
612:       hy     = one/(my-1)
613:       sc     = hx*hy
614:       hxdhy  = hx/hy
615:       hydhx  = hy/hx

617: !  Compute entries of the Jacobian matrix
618: !   - Here, we set all entries for a particular row at once.
619: !   - Note that MatSetValues() uses 0-based row and column numbers
620: !     in Fortran as well as in C.

622:       do 20 j=1,my
623:          row(1) = (j-1)*mx - 1
624:          do 10 i=1,mx
625:             row(1) = row(1) + 1
626: !           boundary points
627:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
628:                PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,[one],INSERT_VALUES,ierr))
629: !           interior grid points
630:             else
631:                v(1) = -hxdhy
632:                v(2) = -hydhx
633:                v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
634:                v(4) = -hydhx
635:                v(5) = -hxdhy
636:                col(1) = row(1) - mx
637:                col(2) = row(1) - 1
638:                col(3) = row(1)
639:                col(4) = row(1) + 1
640:                col(5) = row(1) + mx
641:                PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
642:             endif
643:  10      continue
644:  20   continue

646:       end

648: !
649: !/*TEST
650: !
651: !   build:
652: !      requires: !single
653: !
654: !   test:
655: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
656: !
657: !   test:
658: !      suffix: 2
659: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
660: !
661: !   test:
662: !      suffix: 3
663: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
664: !      filter: sort -b
665: !      filter_output: sort -b
666: !
667: !   test:
668: !     suffix: 4
669: !     args: -pc -par 6.807 -nox
670: !
671: !TEST*/