Actual source code: ex1f.F90

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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*/


 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/examples/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

 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

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

 94:       PetscScalar        lx_v(0:1)
 95:       PetscOffset        lx_i

 97: !  Store parameters in common block

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

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

104:       external FormFunction,FormInitialGuess,FormJacobian

106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: !  Initialize program
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

110:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
111:       if (ierr .ne. 0) then
112:         print*,'Unable to initialize PETSc'
113:         stop
114:       endif
115:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
116:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

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

120: !  Initialize problem parameters
121:       i5 = 5
122:       lambda_max = 6.81
123:       lambda_min = 0.0
124:       lambda     = 6.0
125:       mx         = 4
126:       my         = 4
127:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
128:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
129:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
130:       if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,1,'Lambda out of range '); endif
131:       N  = mx*my
132:       pc = PETSC_FALSE
133:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr);

135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: !  Create nonlinear solver context
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

139:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

141:       if (pc .eqv. PETSC_TRUE) then
142:         call SNESSetType(snes,SNESNEWTONTR,ierr)
143:         call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
144:       endif

146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: !  Create vector data structures; set function evaluation routine
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

150:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
151:       call VecSetSizes(x,PETSC_DECIDE,N,ierr)
152:       call VecSetFromOptions(x,ierr)
153:       call VecDuplicate(x,r,ierr)

155: !  Set function evaluation routine and vector.  Whenever the nonlinear
156: !  solver needs to evaluate the nonlinear function, it will call this
157: !  routine.
158: !   - Note that the final routine argument is the user-defined
159: !     context that provides application-specific data for the
160: !     function evaluation routine.

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

164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: !  Create matrix data structure; set Jacobian evaluation routine
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

172:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
173:       if (.not. matrix_free) then
174:         call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
175:       endif

177: !
178: !     This option will cause the Jacobian to be computed via finite differences
179: !    efficiently using a coloring of the columns of the matrix.
180: !
181:       fd_coloring = .false.
182:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
183:       if (fd_coloring) then

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

215:       else if (.not. matrix_free) then

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

234: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: !  Customize nonlinear solver; set runtime options
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

240:       call SNESSetFromOptions(snes,ierr)

242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: !  Evaluate initial guess; then solve nonlinear system.
244: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

251:       call FormInitialGuess(x,ierr)
252:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
253:       call SNESGetIterationNumber(snes,its,ierr);
254:       if (rank .eq. 0) then
255:          write(6,100) its
256:       endif
257:   100 format('Number of SNES iterations = ',i1)

259: !  PetscDraw contour plot of solution

261:       call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
262:       call PetscDrawSetFromOptions(draw,ierr)

264:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
265:       call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
266:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

268: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: !  Free work space.  All PETSc objects should be destroyed when they
270: !  are no longer needed.
271: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

273:       if (.not. matrix_free) call MatDestroy(J,ierr)
274:       if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)

276:       call VecDestroy(x,ierr)
277:       call VecDestroy(r,ierr)
278:       call SNESDestroy(snes,ierr)
279:       call PetscDrawDestroy(draw,ierr)
280:       call PetscFinalize(ierr)
281:       end

283: ! ---------------------------------------------------------------------
284: !
285: !  FormInitialGuess - Forms initial approximation.
286: !
287: !  Input Parameter:
288: !  X - vector
289: !
290: !  Output Parameters:
291: !  X - vector
292: !  ierr - error code
293: !
294: !  Notes:
295: !  This routine serves as a wrapper for the lower-level routine
296: !  "ApplicationInitialGuess", where the actual computations are
297: !  done using the standard Fortran style of treating the local
298: !  vector data as a multidimensional array over the local mesh.
299: !  This routine merely accesses the local vector data via
300: !  VecGetArray() and VecRestoreArray().
301: !
302:       subroutine FormInitialGuess(X,ierr)
303:       use petscsnes
304:       implicit none

306: !  Input/output variables:
307:       Vec           X
308:       PetscErrorCode    ierr

310: !  Declarations for use with local arrays:
311:       PetscScalar   lx_v(0:1)
312:       PetscOffset   lx_i

314:       0

316: !  Get a pointer to vector data.
317: !    - For default PETSc vectors, VecGetArray() returns a pointer to
318: !      the data array.  Otherwise, the routine is implementation dependent.
319: !    - You MUST call VecRestoreArray() when you no longer need access to
320: !      the array.
321: !    - Note that the Fortran interface to VecGetArray() differs from the
322: !      C version.  See the users manual for details.

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

326: !  Compute initial guess

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

330: !  Restore vector

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

334:       return
335:       end

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

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

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

366: !  Local variables:
367:       PetscInt     i,j
368:       PetscReal temp1,temp,hx,hy,one

370: !  Set parameters

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

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

389:       return
390:       end

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

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

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

429: !  Declarations for use with local arrays:
430:       PetscScalar       lx_v(0:1),lf_v(0:1)
431:       PetscOffset       lx_i,lf_i

433:       PetscInt, pointer :: indices(:)

435: !  Get pointers to vector data.
436: !    - For default PETSc vectors, VecGetArray() returns a pointer to
437: !      the data array.  Otherwise, the routine is implementation dependent.
438: !    - You MUST call VecRestoreArray() when you no longer need access to
439: !      the array.
440: !    - Note that the Fortran interface to VecGetArray() differs from the
441: !      C version.  See the Fortran chapter of the users manual for details.

443:       call VecGetArrayRead(X,lx_v,lx_i,ierr)
444:       call VecGetArray(F,lf_v,lf_i,ierr)

446: !  Compute function

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

450: !  Restore vectors

452:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
453:       call VecRestoreArray(F,lf_v,lf_i,ierr)

455:       call PetscLogFlops(11.0d0*mx*my,ierr)
456: !
457: !     fdcoloring is in the common block and used here ONLY to test the
458: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
459: !
460:       if (fd_coloring) then
461:          call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
462:          print*,'Indices from GetPerturbedColumnsF90'
463:          print*,indices
464:          call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
465:       endif
466:       return
467:       end

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

488: !  Common blocks:
489:       PetscReal      lambda
490:       PetscInt        mx,my
491:       PetscBool         fd_coloring
492:       common         /params/ lambda,mx,my,fd_coloring

494: !  Input/output variables:
495:       PetscScalar    x(mx,my),f(mx,my)
496:       PetscErrorCode       ierr

498: !  Local variables:
499:       PetscScalar    two,one,hx,hy
500:       PetscScalar    hxdhy,hydhx,sc
501:       PetscScalar    u,uxx,uyy
502:       PetscInt        i,j

504:       0
505:       one    = 1.0
506:       two    = 2.0
507:       hx     = one/(mx-1)
508:       hy     = one/(my-1)
509:       sc     = hx*hy*lambda
510:       hxdhy  = hx/hy
511:       hydhx  = hy/hx

513: !  Compute function

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

528:       return
529:       end

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

558: !  Input/output variables:
559:       SNES          snes
560:       Vec           X
561:       Mat           jac,jac_prec
562:       PetscErrorCode      ierr
563:       integer dummy

565: !  Common blocks:
566:       PetscReal     lambda
567:       PetscInt       mx,my
568:       PetscBool         fd_coloring
569:       common        /params/ lambda,mx,my,fd_coloring

571: !  Declarations for use with local array:
572:       PetscScalar   lx_v(0:1)
573:       PetscOffset   lx_i

575: !  Get a pointer to vector data

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

579: !  Compute Jacobian entries

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

583: !  Restore vector

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

587: !  Assemble matrix

589:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
590:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)


593:       return
594:       end

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

616: !  Common blocks:
617:       PetscReal    lambda
618:       PetscInt      mx,my
619:       PetscBool         fd_coloring
620:       common       /params/ lambda,mx,my,fd_coloring

622: !  Input/output variables:
623:       PetscScalar  x(mx,my)
624:       Mat          jac,jac_prec
625:       PetscErrorCode      ierr

627: !  Local variables:
628:       PetscInt      i,j,row(1),col(5),i1,i5
629:       PetscScalar  two,one, hx,hy
630:       PetscScalar  hxdhy,hydhx,sc,v(5)

632: !  Set parameters

634:       i1     = 1
635:       i5     = 5
636:       one    = 1.0
637:       two    = 2.0
638:       hx     = one/(mx-1)
639:       hy     = one/(my-1)
640:       sc     = hx*hy
641:       hxdhy  = hx/hy
642:       hydhx  = hy/hx

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

649:       do 20 j=1,my
650:          row(1) = (j-1)*mx - 1
651:          do 10 i=1,mx
652:             row(1) = row(1) + 1
653: !           boundary points
654:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
655:                call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)
656: !           interior grid points
657:             else
658:                v(1) = -hxdhy
659:                v(2) = -hydhx
660:                v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
661:                v(4) = -hydhx
662:                v(5) = -hxdhy
663:                col(1) = row(1) - mx
664:                col(2) = row(1) - 1
665:                col(3) = row(1)
666:                col(4) = row(1) + 1
667:                col(5) = row(1) + mx
668:                call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)
669:             endif
670:  10      continue
671:  20   continue

673:       return
674:       end

676: !
677: !/*TEST
678: !
679: !   build:
680: !      requires: !single
681: !
682: !   test:
683: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
684: !
685: !   test:
686: !      suffix: 2
687: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
688: !
689: !   test:
690: !      suffix: 3
691: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
692: !      filter: sort -b
693: !      filter_output: sort -b
694: !
695: !   test:
696: !     suffix: 4
697: !     args: -pc -par 6.807 -nox
698: !
699: !TEST*/