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

 42:       call VecDuplicate(x,tmp,ierr)
 43:       mone = -1.0
 44:       call VecWAXPY(tmp,mone,x,w,ierr)
 45:       call VecNorm(tmp,NORM_2,norm,ierr)
 46:       call VecDestroy(tmp,ierr)
 47:       print*, 'Norm of search step ',norm
 48:       changed_y = PETSC_FALSE
 49:       changed_w = PETSC_FALSE
 50:       return
 51:       end

 53:       program main
 54: #include <petsc/finclude/petscdraw.h>
 55:       use petscsnes
 56:       implicit none
 57:       interface SNESSetJacobian
 58:       subroutine SNESSetJacobian1(a,b,c,d,e,z)
 59:        use petscsnes
 60:        SNES a
 61:        Mat b
 62:        Mat c
 63:        external d
 64:        MatFDColoring e
 65:        PetscErrorCode z
 66:       end subroutine
 67:       subroutine SNESSetJacobian2(a,b,c,d,e,z)
 68:        use petscsnes
 69:        SNES a
 70:        Mat b
 71:        Mat c
 72:        external d
 73:        integer e
 74:        PetscErrorCode z
 75:       end subroutine
 76:       end interface
 77: !
 78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79: !                   Variable declarations
 80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81: !
 82: !  Variables:
 83: !     snes        - nonlinear solver
 84: !     x, r        - solution, residual vectors
 85: !     J           - Jacobian matrix
 86: !     its         - iterations for convergence
 87: !     matrix_free - flag - 1 indicates matrix-free version
 88: !     lambda      - nonlinearity parameter
 89: !     draw        - drawing context
 90: !
 91:       SNES               snes
 92:       MatColoring        mc
 93:       Vec                x,r
 94:       PetscDraw               draw
 95:       Mat                J
 96:       PetscBool  matrix_free,flg,fd_coloring
 97:       PetscErrorCode ierr
 98:       PetscInt   its,N, mx,my,i5
 99:       PetscMPIInt size,rank
100:       PetscReal   lambda_max,lambda_min,lambda
101:       MatFDColoring      fdcoloring
102:       ISColoring         iscoloring
103:       PetscBool          pc
104:       external           postcheck

106:       PetscScalar        lx_v(0:1)
107:       PetscOffset        lx_i

109: !  Store parameters in common block

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

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

116:       external FormFunction,FormInitialGuess,FormJacobian

118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: !  Initialize program
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

122:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
123:       if (ierr .ne. 0) then
124:         print*,'Unable to initialize PETSc'
125:         stop
126:       endif
127:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
128:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

130:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif

132: !  Initialize problem parameters
133:       i5 = 5
134:       lambda_max = 6.81
135:       lambda_min = 0.0
136:       lambda     = 6.0
137:       mx         = 4
138:       my         = 4
139:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
140:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
141:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
142:       if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif
143:       N  = mx*my
144:       pc = PETSC_FALSE
145:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr);

147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: !  Create nonlinear solver context
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

151:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

153:       if (pc .eqv. PETSC_TRUE) then
154:         call SNESSetType(snes,SNESNEWTONTR,ierr)
155:         call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
156:       endif

158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: !  Create vector data structures; set function evaluation routine
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

162:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
163:       call VecSetSizes(x,PETSC_DECIDE,N,ierr)
164:       call VecSetFromOptions(x,ierr)
165:       call VecDuplicate(x,r,ierr)

167: !  Set function evaluation routine and vector.  Whenever the nonlinear
168: !  solver needs to evaluate the nonlinear function, it will call this
169: !  routine.
170: !   - Note that the final routine argument is the user-defined
171: !     context that provides application-specific data for the
172: !     function evaluation routine.

174:       call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)

176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: !  Create matrix data structure; set Jacobian evaluation routine
178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

184:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
185:       if (.not. matrix_free) then
186:         call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
187:       endif

189: !
190: !     This option will cause the Jacobian to be computed via finite differences
191: !    efficiently using a coloring of the columns of the matrix.
192: !
193:       fd_coloring = .false.
194:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
195:       if (fd_coloring) then

197: !
198: !      This initializes the nonzero structure of the Jacobian. This is artificial
199: !      because clearly if we had a routine to compute the Jacobian we won't need
200: !      to use finite differences.
201: !
202:         call FormJacobian(snes,x,J,J,0,ierr)
203: !
204: !       Color the matrix, i.e. determine groups of columns that share no common
205: !      rows. These columns in the Jacobian can all be computed simultaneously.
206: !
207:         call MatColoringCreate(J,mc,ierr)
208:         call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
209:         call MatColoringSetFromOptions(mc,ierr)
210:         call MatColoringApply(mc,iscoloring,ierr)
211:         call MatColoringDestroy(mc,ierr)
212: !
213: !       Create the data structure that SNESComputeJacobianDefaultColor() uses
214: !       to compute the actual Jacobians via finite differences.
215: !
216:         call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
217:         call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
218:         call MatFDColoringSetFromOptions(fdcoloring,ierr)
219:         call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
220: !
221: !        Tell SNES to use the routine SNESComputeJacobianDefaultColor()
222: !      to compute Jacobians.
223: !
224:         call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
225:         call ISColoringDestroy(iscoloring,ierr)

227:       else if (.not. matrix_free) then

229: !  Set Jacobian matrix data structure and default Jacobian evaluation
230: !  routine.  Whenever the nonlinear solver needs to compute the
231: !  Jacobian matrix, it will call this routine.
232: !   - Note that the final routine argument is the user-defined
233: !     context that provides application-specific data for the
234: !     Jacobian evaluation routine.
235: !   - The user can override with:
236: !      -snes_fd : default finite differencing approximation of Jacobian
237: !      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
238: !                 (unless user explicitly sets preconditioner)
239: !      -snes_mf_operator : form preconditioning matrix as set by the user,
240: !                          but use matrix-free approx for Jacobian-vector
241: !                          products within Newton-Krylov method
242: !
243:         call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
244:       endif

246: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247: !  Customize nonlinear solver; set runtime options
248: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

252:       call SNESSetFromOptions(snes,ierr)

254: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: !  Evaluate initial guess; then solve nonlinear system.
256: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

263:       call FormInitialGuess(x,ierr)
264:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
265:       call SNESGetIterationNumber(snes,its,ierr);
266:       if (rank .eq. 0) then
267:          write(6,100) its
268:       endif
269:   100 format('Number of SNES iterations = ',i1)

271: !  PetscDraw contour plot of solution

273:       call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
274:       call PetscDrawSetFromOptions(draw,ierr)

276:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
277:       call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
278:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: !  Free work space.  All PETSc objects should be destroyed when they
282: !  are no longer needed.
283: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

285:       if (.not. matrix_free) call MatDestroy(J,ierr)
286:       if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)

288:       call VecDestroy(x,ierr)
289:       call VecDestroy(r,ierr)
290:       call SNESDestroy(snes,ierr)
291:       call PetscDrawDestroy(draw,ierr)
292:       call PetscFinalize(ierr)
293:       end

295: ! ---------------------------------------------------------------------
296: !
297: !  FormInitialGuess - Forms initial approximation.
298: !
299: !  Input Parameter:
300: !  X - vector
301: !
302: !  Output Parameters:
303: !  X - vector
304: !  ierr - error code
305: !
306: !  Notes:
307: !  This routine serves as a wrapper for the lower-level routine
308: !  "ApplicationInitialGuess", where the actual computations are
309: !  done using the standard Fortran style of treating the local
310: !  vector data as a multidimensional array over the local mesh.
311: !  This routine merely accesses the local vector data via
312: !  VecGetArray() and VecRestoreArray().
313: !
314:       subroutine FormInitialGuess(X,ierr)
315:       use petscsnes
316:       implicit none

318: !  Input/output variables:
319:       Vec           X
320:       PetscErrorCode    ierr

322: !  Declarations for use with local arrays:
323:       PetscScalar   lx_v(0:1)
324:       PetscOffset   lx_i

326:       0

328: !  Get a pointer to vector data.
329: !    - For default PETSc vectors, VecGetArray() returns a pointer to
330: !      the data array.  Otherwise, the routine is implementation dependent.
331: !    - You MUST call VecRestoreArray() when you no longer need access to
332: !      the array.
333: !    - Note that the Fortran interface to VecGetArray() differs from the
334: !      C version.  See the users manual for details.

336:       call VecGetArray(X,lx_v,lx_i,ierr)

338: !  Compute initial guess

340:       call ApplicationInitialGuess(lx_v(lx_i),ierr)

342: !  Restore vector

344:       call VecRestoreArray(X,lx_v,lx_i,ierr)

346:       return
347:       end

349: ! ---------------------------------------------------------------------
350: !
351: !  ApplicationInitialGuess - Computes initial approximation, called by
352: !  the higher level routine FormInitialGuess().
353: !
354: !  Input Parameter:
355: !  x - local vector data
356: !
357: !  Output Parameters:
358: !  f - local vector data, f(x)
359: !  ierr - error code
360: !
361: !  Notes:
362: !  This routine uses standard Fortran-style computations over a 2-dim array.
363: !
364:       subroutine ApplicationInitialGuess(x,ierr)
365:       use petscksp
366:       implicit none

368: !  Common blocks:
369:       PetscReal   lambda
370:       PetscInt     mx,my
371:       PetscBool         fd_coloring
372:       common      /params/ lambda,mx,my,fd_coloring

374: !  Input/output variables:
375:       PetscScalar x(mx,my)
376:       PetscErrorCode     ierr

378: !  Local variables:
379:       PetscInt     i,j
380:       PetscReal temp1,temp,hx,hy,one

382: !  Set parameters

384:       0
385:       one    = 1.0
386:       hx     = one/(mx-1)
387:       hy     = one/(my-1)
388:       temp1  = lambda/(lambda + one)

390:       do 20 j=1,my
391:          temp = min(j-1,my-j)*hy
392:          do 10 i=1,mx
393:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
394:               x(i,j) = 0.0
395:             else
396:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
397:             endif
398:  10      continue
399:  20   continue

401:       return
402:       end

404: ! ---------------------------------------------------------------------
405: !
406: !  FormFunction - Evaluates nonlinear function, F(x).
407: !
408: !  Input Parameters:
409: !  snes  - the SNES context
410: !  X     - input vector
411: !  dummy - optional user-defined context, as set by SNESSetFunction()
412: !          (not used here)
413: !
414: !  Output Parameter:
415: !  F     - vector with newly computed function
416: !
417: !  Notes:
418: !  This routine serves as a wrapper for the lower-level routine
419: !  "ApplicationFunction", where the actual computations are
420: !  done using the standard Fortran style of treating the local
421: !  vector data as a multidimensional array over the local mesh.
422: !  This routine merely accesses the local vector data via
423: !  VecGetArray() and VecRestoreArray().
424: !
425:       subroutine FormFunction(snes,X,F,fdcoloring,ierr)
426:       use petscsnes
427:       implicit none

429: !  Input/output variables:
430:       SNES              snes
431:       Vec               X,F
432:       PetscErrorCode          ierr
433:       MatFDColoring fdcoloring

435: !  Common blocks:
436:       PetscReal         lambda
437:       PetscInt          mx,my
438:       PetscBool         fd_coloring
439:       common            /params/ lambda,mx,my,fd_coloring

441: !  Declarations for use with local arrays:
442:       PetscScalar       lx_v(0:1),lf_v(0:1)
443:       PetscOffset       lx_i,lf_i

445:       PetscInt, pointer :: indices(:)

447: !  Get pointers to vector data.
448: !    - For default PETSc vectors, VecGetArray() returns a pointer to
449: !      the data array.  Otherwise, the routine is implementation dependent.
450: !    - You MUST call VecRestoreArray() when you no longer need access to
451: !      the array.
452: !    - Note that the Fortran interface to VecGetArray() differs from the
453: !      C version.  See the Fortran chapter of the users manual for details.

455:       call VecGetArrayRead(X,lx_v,lx_i,ierr)
456:       call VecGetArray(F,lf_v,lf_i,ierr)

458: !  Compute function

460:       call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)

462: !  Restore vectors

464:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
465:       call VecRestoreArray(F,lf_v,lf_i,ierr)

467:       call PetscLogFlops(11.0d0*mx*my,ierr)
468: !
469: !     fdcoloring is in the common block and used here ONLY to test the
470: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
471: !
472:       if (fd_coloring) then
473:          call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
474:          print*,'Indices from GetPerturbedColumnsF90'
475:          write(*,1000) indices
476:  1000    format(50i4)
477:          call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
478:       endif
479:       return
480:       end

482: ! ---------------------------------------------------------------------
483: !
484: !  ApplicationFunction - Computes nonlinear function, called by
485: !  the higher level routine FormFunction().
486: !
487: !  Input Parameter:
488: !  x    - local vector data
489: !
490: !  Output Parameters:
491: !  f    - local vector data, f(x)
492: !  ierr - error code
493: !
494: !  Notes:
495: !  This routine uses standard Fortran-style computations over a 2-dim array.
496: !
497:       subroutine ApplicationFunction(x,f,ierr)
498:       use petscsnes
499:       implicit none

501: !  Common blocks:
502:       PetscReal      lambda
503:       PetscInt        mx,my
504:       PetscBool         fd_coloring
505:       common         /params/ lambda,mx,my,fd_coloring

507: !  Input/output variables:
508:       PetscScalar    x(mx,my),f(mx,my)
509:       PetscErrorCode       ierr

511: !  Local variables:
512:       PetscScalar    two,one,hx,hy
513:       PetscScalar    hxdhy,hydhx,sc
514:       PetscScalar    u,uxx,uyy
515:       PetscInt        i,j

517:       0
518:       one    = 1.0
519:       two    = 2.0
520:       hx     = one/(mx-1)
521:       hy     = one/(my-1)
522:       sc     = hx*hy*lambda
523:       hxdhy  = hx/hy
524:       hydhx  = hy/hx

526: !  Compute function

528:       do 20 j=1,my
529:          do 10 i=1,mx
530:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
531:                f(i,j) = x(i,j)
532:             else
533:                u = x(i,j)
534:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
535:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
536:                f(i,j) = uxx + uyy - sc*exp(u)
537:             endif
538:  10      continue
539:  20   continue

541:       return
542:       end

544: ! ---------------------------------------------------------------------
545: !
546: !  FormJacobian - Evaluates Jacobian matrix.
547: !
548: !  Input Parameters:
549: !  snes    - the SNES context
550: !  x       - input vector
551: !  dummy   - optional user-defined context, as set by SNESSetJacobian()
552: !            (not used here)
553: !
554: !  Output Parameters:
555: !  jac      - Jacobian matrix
556: !  jac_prec - optionally different preconditioning matrix (not used here)
557: !  flag     - flag indicating matrix structure
558: !
559: !  Notes:
560: !  This routine serves as a wrapper for the lower-level routine
561: !  "ApplicationJacobian", where the actual computations are
562: !  done using the standard Fortran style of treating the local
563: !  vector data as a multidimensional array over the local mesh.
564: !  This routine merely accesses the local vector data via
565: !  VecGetArray() and VecRestoreArray().
566: !
567:       subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
568:       use petscsnes
569:       implicit none

571: !  Input/output variables:
572:       SNES          snes
573:       Vec           X
574:       Mat           jac,jac_prec
575:       PetscErrorCode      ierr
576:       integer dummy

578: !  Common blocks:
579:       PetscReal     lambda
580:       PetscInt       mx,my
581:       PetscBool         fd_coloring
582:       common        /params/ lambda,mx,my,fd_coloring

584: !  Declarations for use with local array:
585:       PetscScalar   lx_v(0:1)
586:       PetscOffset   lx_i

588: !  Get a pointer to vector data

590:       call VecGetArrayRead(X,lx_v,lx_i,ierr)

592: !  Compute Jacobian entries

594:       call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)

596: !  Restore vector

598:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)

600: !  Assemble matrix

602:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
603:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)

605:       return
606:       end

608: ! ---------------------------------------------------------------------
609: !
610: !  ApplicationJacobian - Computes Jacobian matrix, called by
611: !  the higher level routine FormJacobian().
612: !
613: !  Input Parameters:
614: !  x        - local vector data
615: !
616: !  Output Parameters:
617: !  jac      - Jacobian matrix
618: !  jac_prec - optionally different preconditioning matrix (not used here)
619: !  ierr     - error code
620: !
621: !  Notes:
622: !  This routine uses standard Fortran-style computations over a 2-dim array.
623: !
624:       subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
625:       use petscsnes
626:       implicit none

628: !  Common blocks:
629:       PetscReal    lambda
630:       PetscInt      mx,my
631:       PetscBool         fd_coloring
632:       common       /params/ lambda,mx,my,fd_coloring

634: !  Input/output variables:
635:       PetscScalar  x(mx,my)
636:       Mat          jac,jac_prec
637:       PetscErrorCode      ierr

639: !  Local variables:
640:       PetscInt      i,j,row(1),col(5),i1,i5
641:       PetscScalar  two,one, hx,hy
642:       PetscScalar  hxdhy,hydhx,sc,v(5)

644: !  Set parameters

646:       i1     = 1
647:       i5     = 5
648:       one    = 1.0
649:       two    = 2.0
650:       hx     = one/(mx-1)
651:       hy     = one/(my-1)
652:       sc     = hx*hy
653:       hxdhy  = hx/hy
654:       hydhx  = hy/hx

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

661:       do 20 j=1,my
662:          row(1) = (j-1)*mx - 1
663:          do 10 i=1,mx
664:             row(1) = row(1) + 1
665: !           boundary points
666:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
667:                call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
668: !           interior grid points
669:             else
670:                v(1) = -hxdhy
671:                v(2) = -hydhx
672:                v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
673:                v(4) = -hydhx
674:                v(5) = -hxdhy
675:                col(1) = row(1) - mx
676:                col(2) = row(1) - 1
677:                col(3) = row(1)
678:                col(4) = row(1) + 1
679:                col(5) = row(1) + mx
680:                call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
681:             endif
682:  10      continue
683:  20   continue

685:       return
686:       end

688: !
689: !/*TEST
690: !
691: !   build:
692: !      requires: !single
693: !
694: !   test:
695: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
696: !
697: !   test:
698: !      suffix: 2
699: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
700: !
701: !   test:
702: !      suffix: 3
703: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
704: !      filter: sort -b
705: !      filter_output: sort -b
706: !
707: !   test:
708: !     suffix: 4
709: !     args: -pc -par 6.807 -nox
710: !
711: !TEST*/