Actual source code: ex1f.F90

petsc-3.10.2 2018-10-09
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: !  --------------------------------------------------------------------------

 36:       program main
 37:  #include <petsc/finclude/petscdraw.h>
 38:  #include <petsc/finclude/petscsnes.h>
 39:       use petscsnes
 40:       implicit none
 41: #if defined(PETSC_USING_F90) && !defined(PETSC_USE_FORTRANKIND)
 42:       external SNESCOMPUTEJACOBIANDEFAULTCOLOR
 43: #endif
 44: !
 45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 46: !                   Variable declarations
 47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48: !
 49: !  Variables:
 50: !     snes        - nonlinear solver
 51: !     x, r        - solution, residual vectors
 52: !     J           - Jacobian matrix
 53: !     its         - iterations for convergence
 54: !     matrix_free - flag - 1 indicates matrix-free version
 55: !     lambda      - nonlinearity parameter
 56: !     draw        - drawing context
 57: !
 58:       SNES               snes
 59:       MatColoring        mc
 60:       Vec                x,r
 61:       PetscDraw               draw
 62:       Mat                J
 63:       PetscBool  matrix_free,flg,fd_coloring
 64:       PetscErrorCode ierr
 65:       PetscInt   its,N, mx,my,i5
 66:       PetscMPIInt size,rank
 67:       PetscReal   lambda_max,lambda_min,lambda
 68:       MatFDColoring      fdcoloring
 69:       ISColoring         iscoloring

 71:       PetscScalar        lx_v(0:1)
 72:       PetscOffset        lx_i

 74: !  Store parameters in common block

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

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

 81:       external FormFunction,FormInitialGuess,FormJacobian

 83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84: !  Initialize program
 85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 87:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 88:       if (ierr .ne. 0) then
 89:         print*,'Unable to initialize PETSc'
 90:         stop
 91:       endif
 92:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 93:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

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

 97: !  Initialize problem parameters
 98:       i5 = 5
 99:       lambda_max = 6.81
100:       lambda_min = 0.0
101:       lambda     = 6.0
102:       mx         = 4
103:       my         = 4
104:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
105:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
106:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
107:       if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,1,'Lambda out of range '); endif
108:       N       = mx*my

110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: !  Create nonlinear solver context
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

114:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: !  Create vector data structures; set function evaluation routine
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

120:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
121:       call VecSetSizes(x,PETSC_DECIDE,N,ierr)
122:       call VecSetFromOptions(x,ierr)
123:       call VecDuplicate(x,r,ierr)

125: !  Set function evaluation routine and vector.  Whenever the nonlinear
126: !  solver needs to evaluate the nonlinear function, it will call this
127: !  routine.
128: !   - Note that the final routine argument is the user-defined
129: !     context that provides application-specific data for the
130: !     function evaluation routine.

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

134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: !  Create matrix data structure; set Jacobian evaluation routine
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

142:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)
143:       if (.not. matrix_free) then
144:         call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)
145:       endif

147: !
148: !     This option will cause the Jacobian to be computed via finite differences
149: !    efficiently using a coloring of the columns of the matrix.
150: !
151:       fd_coloring = .false.
152:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)
153:       if (fd_coloring) then

155: !
156: !      This initializes the nonzero structure of the Jacobian. This is artificial
157: !      because clearly if we had a routine to compute the Jacobian we won't need
158: !      to use finite differences.
159: !
160:         call FormJacobian(snes,x,J,J,0,ierr)
161: !
162: !       Color the matrix, i.e. determine groups of columns that share no common
163: !      rows. These columns in the Jacobian can all be computed simulataneously.
164: !
165:         call MatColoringCreate(J,mc,ierr)
166:         call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)
167:         call MatColoringSetFromOptions(mc,ierr)
168:         call MatColoringApply(mc,iscoloring,ierr)
169:         call MatColoringDestroy(mc,ierr)
170: !
171: !       Create the data structure that SNESComputeJacobianDefaultColor() uses
172: !       to compute the actual Jacobians via finite differences.
173: !
174:         call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)
175:         call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)
176:         call MatFDColoringSetFromOptions(fdcoloring,ierr)
177:         call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)
178: !
179: !        Tell SNES to use the routine SNESComputeJacobianDefaultColor()
180: !      to compute Jacobians.
181: !
182:         call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)
183:         call ISColoringDestroy(iscoloring,ierr)

185:       else if (.not. matrix_free) then

187: !  Set Jacobian matrix data structure and default Jacobian evaluation
188: !  routine.  Whenever the nonlinear solver needs to compute the
189: !  Jacobian matrix, it will call this routine.
190: !   - Note that the final routine argument is the user-defined
191: !     context that provides application-specific data for the
192: !     Jacobian evaluation routine.
193: !   - The user can override with:
194: !      -snes_fd : default finite differencing approximation of Jacobian
195: !      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
196: !                 (unless user explicitly sets preconditioner)
197: !      -snes_mf_operator : form preconditioning matrix as set by the user,
198: !                          but use matrix-free approx for Jacobian-vector
199: !                          products within Newton-Krylov method
200: !
201:         call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
202:       endif

204: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: !  Customize nonlinear solver; set runtime options
206: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

210:       call SNESSetFromOptions(snes,ierr)

212: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: !  Evaluate initial guess; then solve nonlinear system.
214: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

221:       call FormInitialGuess(x,ierr)
222:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
223:       call SNESGetIterationNumber(snes,its,ierr);
224:       if (rank .eq. 0) then
225:          write(6,100) its
226:       endif
227:   100 format('Number of SNES iterations = ',i1)

229: !  PetscDraw contour plot of solution

231:       call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)
232:       call PetscDrawSetFromOptions(draw,ierr)

234:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
235:       call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)
236:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

238: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: !  Free work space.  All PETSc objects should be destroyed when they
240: !  are no longer needed.
241: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

243:       if (.not. matrix_free) call MatDestroy(J,ierr)
244:       if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr)

246:       call VecDestroy(x,ierr)
247:       call VecDestroy(r,ierr)
248:       call SNESDestroy(snes,ierr)
249:       call PetscDrawDestroy(draw,ierr)
250:       call PetscFinalize(ierr)
251:       end

253: ! ---------------------------------------------------------------------
254: !
255: !  FormInitialGuess - Forms initial approximation.
256: !
257: !  Input Parameter:
258: !  X - vector
259: !
260: !  Output Parameters:
261: !  X - vector
262: !  ierr - error code
263: !
264: !  Notes:
265: !  This routine serves as a wrapper for the lower-level routine
266: !  "ApplicationInitialGuess", where the actual computations are
267: !  done using the standard Fortran style of treating the local
268: !  vector data as a multidimensional array over the local mesh.
269: !  This routine merely accesses the local vector data via
270: !  VecGetArray() and VecRestoreArray().
271: !
272:       subroutine FormInitialGuess(X,ierr)
273:       use petscsnes
274:       implicit none

276: !  Input/output variables:
277:       Vec           X
278:       PetscErrorCode    ierr

280: !  Declarations for use with local arrays:
281:       PetscScalar   lx_v(0:1)
282:       PetscOffset   lx_i

284:       0

286: !  Get a pointer to vector data.
287: !    - For default PETSc vectors, VecGetArray() returns a pointer to
288: !      the data array.  Otherwise, the routine is implementation dependent.
289: !    - You MUST call VecRestoreArray() when you no longer need access to
290: !      the array.
291: !    - Note that the Fortran interface to VecGetArray() differs from the
292: !      C version.  See the users manual for details.

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

296: !  Compute initial guess

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

300: !  Restore vector

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

304:       return
305:       end

307: ! ---------------------------------------------------------------------
308: !
309: !  ApplicationInitialGuess - Computes initial approximation, called by
310: !  the higher level routine FormInitialGuess().
311: !
312: !  Input Parameter:
313: !  x - local vector data
314: !
315: !  Output Parameters:
316: !  f - local vector data, f(x)
317: !  ierr - error code
318: !
319: !  Notes:
320: !  This routine uses standard Fortran-style computations over a 2-dim array.
321: !
322:       subroutine ApplicationInitialGuess(x,ierr)
323:       use petscksp
324:       implicit none

326: !  Common blocks:
327:       PetscReal   lambda
328:       PetscInt     mx,my
329:       PetscBool         fd_coloring
330:       common      /params/ lambda,mx,my,fd_coloring

332: !  Input/output variables:
333:       PetscScalar x(mx,my)
334:       PetscErrorCode     ierr

336: !  Local variables:
337:       PetscInt     i,j
338:       PetscReal temp1,temp,hx,hy,one

340: !  Set parameters

342:       0
343:       one    = 1.0
344:       hx     = one/(mx-1)
345:       hy     = one/(my-1)
346:       temp1  = lambda/(lambda + one)

348:       do 20 j=1,my
349:          temp = min(j-1,my-j)*hy
350:          do 10 i=1,mx
351:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
352:               x(i,j) = 0.0
353:             else
354:               x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
355:             endif
356:  10      continue
357:  20   continue

359:       return
360:       end

362: ! ---------------------------------------------------------------------
363: !
364: !  FormFunction - Evaluates nonlinear function, F(x).
365: !
366: !  Input Parameters:
367: !  snes  - the SNES context
368: !  X     - input vector
369: !  dummy - optional user-defined context, as set by SNESSetFunction()
370: !          (not used here)
371: !
372: !  Output Parameter:
373: !  F     - vector with newly computed function
374: !
375: !  Notes:
376: !  This routine serves as a wrapper for the lower-level routine
377: !  "ApplicationFunction", where the actual computations are
378: !  done using the standard Fortran style of treating the local
379: !  vector data as a multidimensional array over the local mesh.
380: !  This routine merely accesses the local vector data via
381: !  VecGetArray() and VecRestoreArray().
382: !
383:       subroutine FormFunction(snes,X,F,fdcoloring,ierr)
384:       use petscsnes
385:       implicit none

387: !  Input/output variables:
388:       SNES              snes
389:       Vec               X,F
390:       PetscErrorCode          ierr
391:       MatFDColoring fdcoloring

393: !  Common blocks:
394:       PetscReal         lambda
395:       PetscInt          mx,my
396:       PetscBool         fd_coloring
397:       common            /params/ lambda,mx,my,fd_coloring

399: !  Declarations for use with local arrays:
400:       PetscScalar       lx_v(0:1),lf_v(0:1)
401:       PetscOffset       lx_i,lf_i

403:       PetscInt, pointer :: indices(:)

405: !  Get pointers to vector data.
406: !    - For default PETSc vectors, VecGetArray() returns a pointer to
407: !      the data array.  Otherwise, the routine is implementation dependent.
408: !    - You MUST call VecRestoreArray() when you no longer need access to
409: !      the array.
410: !    - Note that the Fortran interface to VecGetArray() differs from the
411: !      C version.  See the Fortran chapter of the users manual for details.

413:       call VecGetArrayRead(X,lx_v,lx_i,ierr)
414:       call VecGetArray(F,lf_v,lf_i,ierr)

416: !  Compute function

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

420: !  Restore vectors

422:       call VecRestoreArrayRead(X,lx_v,lx_i,ierr)
423:       call VecRestoreArray(F,lf_v,lf_i,ierr)

425:       call PetscLogFlops(11.0d0*mx*my,ierr)
426: !
427: !     fdcoloring is in the common block and used here ONLY to test the
428: !     calls to MatFDColoringGetPerturbedColumnsF90() and  MatFDColoringRestorePerturbedColumnsF90()
429: !
430:       if (fd_coloring) then
431:          call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)
432:          print*,'Indices from GetPerturbedColumnsF90'
433:          print*,indices
434:          call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)
435:       endif
436:       return
437:       end

439: ! ---------------------------------------------------------------------
440: !
441: !  ApplicationFunction - Computes nonlinear function, called by
442: !  the higher level routine FormFunction().
443: !
444: !  Input Parameter:
445: !  x    - local vector data
446: !
447: !  Output Parameters:
448: !  f    - local vector data, f(x)
449: !  ierr - error code
450: !
451: !  Notes:
452: !  This routine uses standard Fortran-style computations over a 2-dim array.
453: !
454:       subroutine ApplicationFunction(x,f,ierr)
455:       use petscsnes
456:       implicit none

458: !  Common blocks:
459:       PetscReal      lambda
460:       PetscInt        mx,my
461:       PetscBool         fd_coloring
462:       common         /params/ lambda,mx,my,fd_coloring

464: !  Input/output variables:
465:       PetscScalar    x(mx,my),f(mx,my)
466:       PetscErrorCode       ierr

468: !  Local variables:
469:       PetscScalar    two,one,hx,hy
470:       PetscScalar    hxdhy,hydhx,sc
471:       PetscScalar    u,uxx,uyy
472:       PetscInt        i,j

474:       0
475:       one    = 1.0
476:       two    = 2.0
477:       hx     = one/(mx-1)
478:       hy     = one/(my-1)
479:       sc     = hx*hy*lambda
480:       hxdhy  = hx/hy
481:       hydhx  = hy/hx

483: !  Compute function

485:       do 20 j=1,my
486:          do 10 i=1,mx
487:             if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
488:                f(i,j) = x(i,j)
489:             else
490:                u = x(i,j)
491:                uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
492:                uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
493:                f(i,j) = uxx + uyy - sc*exp(u)
494:             endif
495:  10      continue
496:  20   continue

498:       return
499:       end

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

528: !  Input/output variables:
529:       SNES          snes
530:       Vec           X
531:       Mat           jac,jac_prec
532:       PetscErrorCode      ierr
533:       integer dummy

535: !  Common blocks:
536:       PetscReal     lambda
537:       PetscInt       mx,my
538:       PetscBool         fd_coloring
539:       common        /params/ lambda,mx,my,fd_coloring

541: !  Declarations for use with local array:
542:       PetscScalar   lx_v(0:1)
543:       PetscOffset   lx_i

545: !  Get a pointer to vector data

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

549: !  Compute Jacobian entries

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

553: !  Restore vector

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

557: !  Assemble matrix

559:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)
560:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)


563:       return
564:       end

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

586: !  Common blocks:
587:       PetscReal    lambda
588:       PetscInt      mx,my
589:       PetscBool         fd_coloring
590:       common       /params/ lambda,mx,my,fd_coloring

592: !  Input/output variables:
593:       PetscScalar  x(mx,my)
594:       Mat          jac,jac_prec
595:       PetscErrorCode      ierr

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

602: !  Set parameters

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

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

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

643:       return
644:       end

646: !
647: !/*TEST
648: !
649: !   build:
650: !      requires: !single
651: !
652: !   test:
653: !      args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
654: !
655: !   test:
656: !      suffix: 2
657: !      args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
658: !
659: !   test:
660: !      suffix: 3
661: !      args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
662: !      filter: sort -b
663: !      filter_output: sort -b
664: !
665: !TEST*/