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

108:       character(len=PETSC_MAX_PATH_LEN) :: outputString

110:       PetscScalar,pointer :: lx_v(:)

112: !  Store parameters in common block

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

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

119:       external FormFunction,FormInitialGuess,FormJacobian

121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: !  Initialize program
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

125:       PetscCallA(PetscInitialize(ierr))
126:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
127:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))

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

131: !  Initialize problem parameters
132:       i5 = 5
133:       lambda_max = 6.81
134:       lambda_min = 0.0
135:       lambda     = 6.0
136:       mx         = 4
137:       my         = 4
138:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
139:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
140:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
141:       PetscCheckA(lambda .lt. lambda_max .and. lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range ')
142:       N  = mx*my
143:       pc = PETSC_FALSE
144:       PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))

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

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

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

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

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

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

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

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

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

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

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

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

226:       else if (.not. matrix_free) then

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

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

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

251:       PetscCallA(SNESSetFromOptions(snes,ierr))

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

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

262:       PetscCallA(FormInitialGuess(x,ierr))
263:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
264:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))
265:       write(outputString,*) its
266:       PetscCallA(PetscPrintf(PETSC_COMM_WORLD,'Number of SNES iterations = '//trim(outputString)//'\n',ierr))

268: !  PetscDraw contour plot of solution

270:       PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
271:       PetscCallA(PetscDrawSetFromOptions(draw,ierr))

273:       PetscCallA(VecGetArrayReadF90(x,lx_v,ierr))
274:       PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v,ierr))
275:       PetscCallA(VecRestoreArrayReadF90(x,lx_v,ierr))

277: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: !  Free work space.  All PETSc objects should be destroyed when they
279: !  are no longer needed.
280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

282:       if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
283:       if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))

285:       PetscCallA(VecDestroy(x,ierr))
286:       PetscCallA(VecDestroy(r,ierr))
287:       PetscCallA(SNESDestroy(snes,ierr))
288:       PetscCallA(PetscDrawDestroy(draw,ierr))
289:       PetscCallA(PetscFinalize(ierr))
290:       end

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

315: !  Input/output variables:
316:       Vec           X
317:       PetscErrorCode    ierr

319: !     Declarations for use with local arrays:
320:       PetscScalar,pointer :: lx_v(:)

322:       ierr   = 0

324: !  Get a pointer to vector data.
325: !    - VecGetArrayF90() returns a pointer to the data array.
326: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
327: !      the array.

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

331: !  Compute initial guess

333:       PetscCallA(ApplicationInitialGuess(lx_v,ierr))

335: !  Restore vector

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

339:       end

341: !  ApplicationInitialGuess - Computes initial approximation, called by
342: !  the higher level routine FormInitialGuess().
343: !
344: !  Input Parameter:
345: !  x - local vector data
346: !
347: !  Output Parameters:
348: !  f - local vector data, f(x)
349: !  ierr - error code
350: !
351: !  Notes:
352: !  This routine uses standard Fortran-style computations over a 2-dim array.
353: !
354:       subroutine ApplicationInitialGuess(x,ierr)
355:       use petscksp
356:       implicit none

358: !  Common blocks:
359:       PetscReal   lambda
360:       PetscInt     mx,my
361:       PetscBool         fd_coloring
362:       common      /params/ lambda,mx,my,fd_coloring

364: !  Input/output variables:
365:       PetscScalar x(mx,my)
366:       PetscErrorCode     ierr

368: !  Local variables:
369:       PetscInt     i,j
370:       PetscReal temp1,temp,hx,hy,one

372: !  Set parameters

374:       ierr   = 0
375:       one    = 1.0
376:       hx     = one/(mx-1)
377:       hy     = one/(my-1)
378:       temp1  = lambda/(lambda + one)

380:       do 20 j=1,my
381:          temp = min(j-1,my-j)*hy
382:          do 10 i=1,mx
383:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
384:               x(i,j) = 0.0
385:             else
386:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
387:             endif
388:  10      continue
389:  20   continue

391:       end

393: ! ---------------------------------------------------------------------
394: !
395: !  FormFunction - Evaluates nonlinear function, F(x).
396: !
397: !  Input Parameters:
398: !  snes  - the SNES context
399: !  X     - input vector
400: !  dummy - optional user-defined context, as set by SNESSetFunction()
401: !          (not used here)
402: !
403: !  Output Parameter:
404: !  F     - vector with newly computed function
405: !
406: !  Notes:
407: !  This routine serves as a wrapper for the lower-level routine
408: !  "ApplicationFunction", where the actual computations are
409: !  done using the standard Fortran style of treating the local
410: !  vector data as a multidimensional array over the local mesh.
411: !  This routine merely accesses the local vector data via
412: !  VecGetArrayF90() and VecRestoreArrayF90().
413: !
414:       subroutine FormFunction(snes,X,F,fdcoloring,ierr)
415:       use petscsnes
416:       implicit none

418: !  Input/output variables:
419:       SNES              snes
420:       Vec               X,F
421:       PetscErrorCode          ierr
422:       MatFDColoring fdcoloring

424: !  Common blocks:
425:       PetscReal         lambda
426:       PetscInt          mx,my
427:       PetscBool         fd_coloring
428:       common            /params/ lambda,mx,my,fd_coloring

430: !  Declarations for use with local arrays:
431:       PetscScalar,pointer :: lx_v(:), lf_v(:)
432:       PetscInt, pointer :: indices(:)

434: !  Get pointers to vector data.
435: !    - VecGetArrayF90() returns a pointer to the data array.
436: !    - You MUST call VecRestoreArrayF90() when you no longer need access to
437: !      the array.

439:       PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
440:       PetscCallA(VecGetArrayF90(F,lf_v,ierr))

442: !  Compute function

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

446: !  Restore vectors

448:       PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
449:       PetscCallA(VecRestoreArrayF90(F,lf_v,ierr))

451:       PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
452: !
453: !     fdcoloring is in the common block and used here ONLY to test the
454: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
455: !
456:       if (fd_coloring) then
457:          PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
458:          print*,'Indices from GetPerturbedColumnsF90'
459:          write(*,1000) indices
460:  1000    format(50i4)
461:          PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
462:       endif
463:       end

465: ! ---------------------------------------------------------------------
466: !
467: !  ApplicationFunction - Computes nonlinear function, called by
468: !  the higher level routine FormFunction().
469: !
470: !  Input Parameter:
471: !  x    - local vector data
472: !
473: !  Output Parameters:
474: !  f    - local vector data, f(x)
475: !  ierr - error code
476: !
477: !  Notes:
478: !  This routine uses standard Fortran-style computations over a 2-dim array.
479: !
480:       subroutine ApplicationFunction(x,f,ierr)
481:       use petscsnes
482:       implicit none

484: !  Common blocks:
485:       PetscReal      lambda
486:       PetscInt        mx,my
487:       PetscBool         fd_coloring
488:       common         /params/ lambda,mx,my,fd_coloring

490: !  Input/output variables:
491:       PetscScalar    x(mx,my),f(mx,my)
492:       PetscErrorCode       ierr

494: !  Local variables:
495:       PetscScalar    two,one,hx,hy
496:       PetscScalar    hxdhy,hydhx,sc
497:       PetscScalar    u,uxx,uyy
498:       PetscInt        i,j

500:       ierr   = 0
501:       one    = 1.0
502:       two    = 2.0
503:       hx     = one/(mx-1)
504:       hy     = one/(my-1)
505:       sc     = hx*hy*lambda
506:       hxdhy  = hx/hy
507:       hydhx  = hy/hx

509: !  Compute function

511:       do 20 j=1,my
512:          do 10 i=1,mx
513:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
514:                f(i,j) = x(i,j)
515:             else
516:                u = x(i,j)
517:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
518:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
519:                f(i,j) = uxx + uyy - sc*exp(u)
520:             endif
521:  10      continue
522:  20   continue

524:       end

526: ! ---------------------------------------------------------------------
527: !
528: !  FormJacobian - Evaluates Jacobian matrix.
529: !
530: !  Input Parameters:
531: !  snes    - the SNES context
532: !  x       - input vector
533: !  dummy   - optional user-defined context, as set by SNESSetJacobian()
534: !            (not used here)
535: !
536: !  Output Parameters:
537: !  jac      - Jacobian matrix
538: !  jac_prec - optionally different preconditioning matrix (not used here)
539: !  flag     - flag indicating matrix structure
540: !
541: !  Notes:
542: !  This routine serves as a wrapper for the lower-level routine
543: !  "ApplicationJacobian", where the actual computations are
544: !  done using the standard Fortran style of treating the local
545: !  vector data as a multidimensional array over the local mesh.
546: !  This routine merely accesses the local vector data via
547: !  VecGetArrayF90() and VecRestoreArrayF90().
548: !
549:       subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
550:       use petscsnes
551:       implicit none

553: !  Input/output variables:
554:       SNES          snes
555:       Vec           X
556:       Mat           jac,jac_prec
557:       PetscErrorCode      ierr
558:       integer dummy

560: !  Common blocks:
561:       PetscReal     lambda
562:       PetscInt       mx,my
563:       PetscBool         fd_coloring
564:       common        /params/ lambda,mx,my,fd_coloring

566: !  Declarations for use with local array:
567:       PetscScalar,pointer :: lx_v(:)

569: !  Get a pointer to vector data

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

573: !  Compute Jacobian entries

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

577: !  Restore vector

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

581: !  Assemble matrix

583:       PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
584:       PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))

586:       end

588: ! ---------------------------------------------------------------------
589: !
590: !  ApplicationJacobian - Computes Jacobian matrix, called by
591: !  the higher level routine FormJacobian().
592: !
593: !  Input Parameters:
594: !  x        - local vector data
595: !
596: !  Output Parameters:
597: !  jac      - Jacobian matrix
598: !  jac_prec - optionally different preconditioning matrix (not used here)
599: !  ierr     - error code
600: !
601: !  Notes:
602: !  This routine uses standard Fortran-style computations over a 2-dim array.
603: !
604:       subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
605:       use petscsnes
606:       implicit none

608: !  Common blocks:
609:       PetscReal    lambda
610:       PetscInt      mx,my
611:       PetscBool         fd_coloring
612:       common       /params/ lambda,mx,my,fd_coloring

614: !  Input/output variables:
615:       PetscScalar  x(mx,my)
616:       Mat          jac,jac_prec
617:       PetscErrorCode      ierr

619: !  Local variables:
620:       PetscInt      i,j,row(1),col(5),i1,i5
621:       PetscScalar  two,one, hx,hy
622:       PetscScalar  hxdhy,hydhx,sc,v(5)

624: !  Set parameters

626:       i1     = 1
627:       i5     = 5
628:       one    = 1.0
629:       two    = 2.0
630:       hx     = one/(mx-1)
631:       hy     = one/(my-1)
632:       sc     = hx*hy
633:       hxdhy  = hx/hy
634:       hydhx  = hy/hx

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

641:       do 20 j=1,my
642:          row(1) = (j-1)*mx - 1
643:          do 10 i=1,mx
644:             row(1) = row(1) + 1
645: !           boundary points
646:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
647:                PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr))
648: !           interior grid points
649:             else
650:                v(1) = -hxdhy
651:                v(2) = -hydhx
652:                v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
653:                v(4) = -hydhx
654:                v(5) = -hxdhy
655:                col(1) = row(1) - mx
656:                col(2) = row(1) - 1
657:                col(3) = row(1)
658:                col(4) = row(1) + 1
659:                col(5) = row(1) + mx
660:                PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
661:             endif
662:  10      continue
663:  20   continue

665:       end

667: !
668: !/*TEST
669: !
670: !   build:
671: !      requires: !single
672: !
673: !   test:
674: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
675: !
676: !   test:
677: !      suffix: 2
678: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
679: !
680: !   test:
681: !      suffix: 3
682: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
683: !      filter: sort -b
684: !      filter_output: sort -b
685: !
686: !   test:
687: !     suffix: 4
688: !     args: -pc -par 6.807 -nox
689: !
690: !TEST*/