Actual source code: ex1f.F90

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

111:       PetscScalar        lx_v(0:1)
112:       PetscOffset        lx_i

114: !  Store parameters in common block

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

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

121:       external FormFunction,FormInitialGuess,FormJacobian

123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: !  Initialize program
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

127:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
128:       if (ierr .ne. 0) then
129:         print*,'Unable to initialize PETSc'
130:         stop
131:       endif
132:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
133:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

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

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

152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: !  Create nonlinear solver context
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

156:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

158:       if (pc .eqv. PETSC_TRUE) then
159:         call SNESSetType(snes,SNESNEWTONTR,ierr)
160:         call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)
161:       endif

163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: !  Create vector data structures; set function evaluation routine
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

167:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
168:       call VecSetSizes(x,PETSC_DECIDE,N,ierr)
169:       call VecSetFromOptions(x,ierr)
170:       call VecDuplicate(x,r,ierr)

172: !  Set function evaluation routine and vector.  Whenever the nonlinear
173: !  solver needs to evaluate the nonlinear function, it will call this
174: !  routine.
175: !   - Note that the final routine argument is the user-defined
176: !     context that provides Section 1.5 Writing Application Codes with PETSc-specific data for the
177: !     function evaluation routine.

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

181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: !  Create matrix data structure; set Jacobian evaluation routine
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

189:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
190:       if (.not. matrix_free) then
191:         call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
192:       endif

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

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

232:       else if (.not. matrix_free) then

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

251: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: !  Customize nonlinear solver; set runtime options
253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

257:       call SNESSetFromOptions(snes,ierr)

259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: !  Evaluate initial guess; then solve nonlinear system.
261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

268:       call FormInitialGuess(x,ierr)
269:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
270:       call SNESGetIterationNumber(snes,its,ierr);
271:       if (rank .eq. 0) then
272:          write(6,100) its
273:       endif
274:   100 format('Number of SNES iterations = ',i1)

276: !  PetscDraw contour plot of solution

278:       call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
279:       call PetscDrawSetFromOptions(draw,ierr)

281:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
282:       call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
283:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

285: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: !  Free work space.  All PETSc objects should be destroyed when they
287: !  are no longer needed.
288: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

290:       if (.not. matrix_free) call MatDestroy(J,ierr)
291:       if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)

293:       call VecDestroy(x,ierr)
294:       call VecDestroy(r,ierr)
295:       call SNESDestroy(snes,ierr)
296:       call PetscDrawDestroy(draw,ierr)
297:       call PetscFinalize(ierr)
298:       end

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

323: !  Input/output variables:
324:       Vec           X
325:       PetscErrorCode    ierr

327: !  Declarations for use with local arrays:
328:       PetscScalar   lx_v(0:1)
329:       PetscOffset   lx_i

331:       0

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

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

343: !  Compute initial guess

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

347: !  Restore vector

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

351:       return
352:       end

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

373: !  Common blocks:
374:       PetscReal   lambda
375:       PetscInt     mx,my
376:       PetscBool         fd_coloring
377:       common      /params/ lambda,mx,my,fd_coloring

379: !  Input/output variables:
380:       PetscScalar x(mx,my)
381:       PetscErrorCode     ierr

383: !  Local variables:
384:       PetscInt     i,j
385:       PetscReal temp1,temp,hx,hy,one

387: !  Set parameters

389:       0
390:       one    = 1.0
391:       hx     = one/(mx-1)
392:       hy     = one/(my-1)
393:       temp1  = lambda/(lambda + one)

395:       do 20 j=1,my
396:          temp = min(j-1,my-j)*hy
397:          do 10 i=1,mx
398:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
399:               x(i,j) = 0.0
400:             else
401:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
402:             endif
403:  10      continue
404:  20   continue

406:       return
407:       end

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

434: !  Input/output variables:
435:       SNES              snes
436:       Vec               X,F
437:       PetscErrorCode          ierr
438:       MatFDColoring fdcoloring

440: !  Common blocks:
441:       PetscReal         lambda
442:       PetscInt          mx,my
443:       PetscBool         fd_coloring
444:       common            /params/ lambda,mx,my,fd_coloring

446: !  Declarations for use with local arrays:
447:       PetscScalar       lx_v(0:1),lf_v(0:1)
448:       PetscOffset       lx_i,lf_i

450:       PetscInt, pointer :: indices(:)

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

460:       call VecGetArrayRead(X,lx_v,lx_i,ierr)
461:       call VecGetArray(F,lf_v,lf_i,ierr)

463: !  Compute function

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

467: !  Restore vectors

469:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
470:       call VecRestoreArray(F,lf_v,lf_i,ierr)

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

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

506: !  Common blocks:
507:       PetscReal      lambda
508:       PetscInt        mx,my
509:       PetscBool         fd_coloring
510:       common         /params/ lambda,mx,my,fd_coloring

512: !  Input/output variables:
513:       PetscScalar    x(mx,my),f(mx,my)
514:       PetscErrorCode       ierr

516: !  Local variables:
517:       PetscScalar    two,one,hx,hy
518:       PetscScalar    hxdhy,hydhx,sc
519:       PetscScalar    u,uxx,uyy
520:       PetscInt        i,j

522:       0
523:       one    = 1.0
524:       two    = 2.0
525:       hx     = one/(mx-1)
526:       hy     = one/(my-1)
527:       sc     = hx*hy*lambda
528:       hxdhy  = hx/hy
529:       hydhx  = hy/hx

531: !  Compute function

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

546:       return
547:       end

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

576: !  Input/output variables:
577:       SNES          snes
578:       Vec           X
579:       Mat           jac,jac_prec
580:       PetscErrorCode      ierr
581:       integer dummy

583: !  Common blocks:
584:       PetscReal     lambda
585:       PetscInt       mx,my
586:       PetscBool         fd_coloring
587:       common        /params/ lambda,mx,my,fd_coloring

589: !  Declarations for use with local array:
590:       PetscScalar   lx_v(0:1)
591:       PetscOffset   lx_i

593: !  Get a pointer to vector data

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

597: !  Compute Jacobian entries

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

601: !  Restore vector

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

605: !  Assemble matrix

607:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
608:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)


611:       return
612:       end

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

634: !  Common blocks:
635:       PetscReal    lambda
636:       PetscInt      mx,my
637:       PetscBool         fd_coloring
638:       common       /params/ lambda,mx,my,fd_coloring

640: !  Input/output variables:
641:       PetscScalar  x(mx,my)
642:       Mat          jac,jac_prec
643:       PetscErrorCode      ierr

645: !  Local variables:
646:       PetscInt      i,j,row(1),col(5),i1,i5
647:       PetscScalar  two,one, hx,hy
648:       PetscScalar  hxdhy,hydhx,sc,v(5)

650: !  Set parameters

652:       i1     = 1
653:       i5     = 5
654:       one    = 1.0
655:       two    = 2.0
656:       hx     = one/(mx-1)
657:       hy     = one/(my-1)
658:       sc     = hx*hy
659:       hxdhy  = hx/hy
660:       hydhx  = hy/hx

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

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

691:       return
692:       end

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