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 petscdraw
58: use petscsnes
59: implicit none
60: !
61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: ! Variable declarations
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: !
65: ! Variables:
66: ! snes - nonlinear solver
67: ! x, r - solution, residual vectors
68: ! J - Jacobian matrix
69: ! its - iterations for convergence
70: ! matrix_free - flag - 1 indicates matrix-free version
71: ! lambda - nonlinearity parameter
72: ! draw - drawing context
73: !
74: SNES snes
75: MatColoring mc
76: Vec x,r
77: PetscDraw draw
78: Mat J
79: PetscBool matrix_free,flg,fd_coloring
80: PetscErrorCode ierr
81: PetscInt its,N, mx,my,i5
82: PetscMPIInt size,rank
83: PetscReal lambda_max,lambda_min,lambda
84: MatFDColoring fdcoloring
85: ISColoring iscoloring
86: PetscBool pc
87: integer4 imx,imy
88: external postcheck
89: character(len=PETSC_MAX_PATH_LEN) :: outputString
90: PetscScalar,pointer :: lx_v(:)
91: integer4 xl,yl,width,height
93: ! Store parameters in common block
95: common /params/ lambda,mx,my,fd_coloring
97: ! Note: Any user-defined Fortran routines (such as FormJacobian)
98: ! MUST be declared as external.
100: external FormFunction,FormInitialGuess,FormJacobian
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: ! Initialize program
104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: PetscCallA(PetscInitialize(ierr))
107: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
108: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
110: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')
112: ! Initialize problem parameters
113: i5 = 5
114: lambda_max = 6.81
115: lambda_min = 0.0
116: lambda = 6.0
117: mx = 4
118: my = 4
119: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
120: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
121: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
122: PetscCheckA(lambda .lt. lambda_max .and. lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range ')
123: N = mx*my
124: pc = PETSC_FALSE
125: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! Create nonlinear solver context
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
133: if (pc .eqv. PETSC_TRUE) then
134: PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
135: PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
136: endif
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! Create vector data structures; set function evaluation routine
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
143: PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
144: PetscCallA(VecSetFromOptions(x,ierr))
145: PetscCallA(VecDuplicate(x,r,ierr))
147: ! Set function evaluation routine and vector. Whenever the nonlinear
148: ! solver needs to evaluate the nonlinear function, it will call this
149: ! routine.
150: ! - Note that the final routine argument is the user-defined
151: ! context that provides application-specific data for the
152: ! function evaluation routine.
154: PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Create matrix data structure; set Jacobian evaluation routine
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Create matrix. Here we only approximately preallocate storage space
161: ! for the Jacobian. See the users manual for a discussion of better
162: ! techniques for preallocating matrix memory.
164: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
165: if (.not. matrix_free) then
166: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER_ARRAY,J,ierr))
167: endif
169: !
170: ! This option will cause the Jacobian to be computed via finite differences
171: ! efficiently using a coloring of the columns of the matrix.
172: !
173: fd_coloring = .false.
174: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
175: if (fd_coloring) then
177: !
178: ! This initializes the nonzero structure of the Jacobian. This is artificial
179: ! because clearly if we had a routine to compute the Jacobian we won't need
180: ! to use finite differences.
181: !
182: PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
183: !
184: ! Color the matrix, i.e. determine groups of columns that share no common
185: ! rows. These columns in the Jacobian can all be computed simultaneously.
186: !
187: PetscCallA(MatColoringCreate(J,mc,ierr))
188: PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
189: PetscCallA(MatColoringSetFromOptions(mc,ierr))
190: PetscCallA(MatColoringApply(mc,iscoloring,ierr))
191: PetscCallA(MatColoringDestroy(mc,ierr))
192: !
193: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
194: ! to compute the actual Jacobians via finite differences.
195: !
196: PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
197: PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
198: PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
199: PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
200: !
201: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
202: ! to compute Jacobians.
203: !
204: PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
205: PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr))
206: PetscCallA(SNESGetJacobian(snes,J,PETSC_NULL_MAT,PETSC_NULL_FUNCTION,fdcoloring,ierr))
207: PetscCallA(ISColoringDestroy(iscoloring,ierr))
209: else if (.not. matrix_free) then
211: ! Set Jacobian matrix data structure and default Jacobian evaluation
212: ! routine. Whenever the nonlinear solver needs to compute the
213: ! Jacobian matrix, it will call this routine.
214: ! - Note that the final routine argument is the user-defined
215: ! context that provides application-specific data for the
216: ! Jacobian evaluation routine.
217: ! - The user can override with:
218: ! -snes_fd : default finite differencing approximation of Jacobian
219: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
220: ! (unless user explicitly sets preconditioner)
221: ! -snes_mf_operator : form preconditioning matrix as set by the user,
222: ! but use matrix-free approx for Jacobian-vector
223: ! products within Newton-Krylov method
224: !
225: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
226: endif
228: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: ! Customize nonlinear solver; set runtime options
230: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
234: PetscCallA(SNESSetFromOptions(snes,ierr))
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: ! Evaluate initial guess; then solve nonlinear system.
238: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240: ! Note: The user should initialize the vector, x, with the initial guess
241: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
242: ! to employ an initial guess of zero, the user should explicitly set
243: ! this vector to zero by calling VecSet().
245: PetscCallA(FormInitialGuess(x,ierr))
246: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
247: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
248: write(outputString,*) its
249: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,'Number of SNES iterations = '//trim(outputString)//'\n',ierr))
251: ! PetscDraw contour plot of solution
253: xl = 300
254: yl = 0
255: width = 300
256: height = 300
257: PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',xl,yl,width,height,draw,ierr))
258: PetscCallA(PetscDrawSetFromOptions(draw,ierr))
260: PetscCallA(VecGetArrayRead(x,lx_v,ierr))
261: imx = int(mx, kind=kind(imx))
262: imy = int(my, kind=kind(imy))
263: PetscCallA(PetscDrawTensorContour(draw,imx,imy,PETSC_NULL_REAL_ARRAY,PETSC_NULL_REAL_ARRAY,lx_v,ierr))
264: PetscCallA(VecRestoreArrayRead(x,lx_v,ierr))
266: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267: ! Free work space. All PETSc objects should be destroyed when they
268: ! are no longer needed.
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271: if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
272: if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))
274: PetscCallA(VecDestroy(x,ierr))
275: PetscCallA(VecDestroy(r,ierr))
276: PetscCallA(SNESDestroy(snes,ierr))
277: PetscCallA(PetscDrawDestroy(draw,ierr))
278: PetscCallA(PetscFinalize(ierr))
279: end
281: ! ---------------------------------------------------------------------
282: !
283: ! FormInitialGuess - Forms initial approximation.
284: !
285: ! Input Parameter:
286: ! X - vector
287: !
288: ! Output Parameters:
289: ! X - vector
290: ! ierr - error code
291: !
292: ! Notes:
293: ! This routine serves as a wrapper for the lower-level routine
294: ! "ApplicationInitialGuess", where the actual computations are
295: ! done using the standard Fortran style of treating the local
296: ! vector data as a multidimensional array over the local mesh.
297: ! This routine merely accesses the local vector data via
298: ! VecGetArray() and VecRestoreArray().
299: !
300: subroutine FormInitialGuess(X,ierr)
301: use petscsnes
302: implicit none
304: ! Input/output variables:
305: Vec X
306: PetscErrorCode ierr
308: ! Declarations for use with local arrays:
309: PetscScalar,pointer :: lx_v(:)
311: ierr = 0
313: ! Get a pointer to vector data.
314: ! - VecGetArray() returns a pointer to the data array.
315: ! - You MUST call VecRestoreArray() when you no longer need access to
316: ! the array.
318: PetscCallA(VecGetArray(X,lx_v,ierr))
320: ! Compute initial guess
322: PetscCallA(ApplicationInitialGuess(lx_v,ierr))
324: ! Restore vector
326: PetscCallA(VecRestoreArray(X,lx_v,ierr))
328: end
330: ! ApplicationInitialGuess - Computes initial approximation, called by
331: ! the higher level routine FormInitialGuess().
332: !
333: ! Input Parameter:
334: ! x - local vector data
335: !
336: ! Output Parameters:
337: ! f - local vector data, f(x)
338: ! ierr - error code
339: !
340: ! Notes:
341: ! This routine uses standard Fortran-style computations over a 2-dim array.
342: !
343: subroutine ApplicationInitialGuess(x,ierr)
344: use petscksp
345: implicit none
347: ! Common blocks:
348: PetscReal lambda
349: PetscInt mx,my
350: PetscBool fd_coloring
351: common /params/ lambda,mx,my,fd_coloring
353: ! Input/output variables:
354: PetscScalar x(mx,my)
355: PetscErrorCode ierr
357: ! Local variables:
358: PetscInt i,j
359: PetscReal temp1,temp,hx,hy,one
361: ! Set parameters
363: ierr = 0
364: one = 1.0
365: hx = one/(mx-1)
366: hy = one/(my-1)
367: temp1 = lambda/(lambda + one)
369: do 20 j=1,my
370: temp = min(j-1,my-j)*hy
371: do 10 i=1,mx
372: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
373: x(i,j) = 0.0
374: else
375: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
376: endif
377: 10 continue
378: 20 continue
380: end
382: ! ---------------------------------------------------------------------
383: !
384: ! FormFunction - Evaluates nonlinear function, F(x).
385: !
386: ! Input Parameters:
387: ! snes - the SNES context
388: ! X - input vector
389: ! dummy - optional user-defined context, as set by SNESSetFunction()
390: ! (not used here)
391: !
392: ! Output Parameter:
393: ! F - vector with newly computed function
394: !
395: ! Notes:
396: ! This routine serves as a wrapper for the lower-level routine
397: ! "ApplicationFunction", where the actual computations are
398: ! done using the standard Fortran style of treating the local
399: ! vector data as a multidimensional array over the local mesh.
400: ! This routine merely accesses the local vector data via
401: ! VecGetArray() and VecRestoreArray().
402: !
403: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
404: use petscsnes
405: implicit none
407: ! Input/output variables:
408: SNES snes
409: Vec X,F
410: PetscErrorCode ierr
411: MatFDColoring fdcoloring
413: ! Common blocks:
414: PetscReal lambda
415: PetscInt mx,my
416: PetscBool fd_coloring
417: common /params/ lambda,mx,my,fd_coloring
419: ! Declarations for use with local arrays:
420: PetscScalar,pointer :: lx_v(:), lf_v(:)
421: PetscInt, pointer :: indices(:)
423: ! Get pointers to vector data.
424: ! - VecGetArray() returns a pointer to the data array.
425: ! - You MUST call VecRestoreArray() when you no longer need access to
426: ! the array.
428: PetscCallA(VecGetArrayRead(X,lx_v,ierr))
429: PetscCallA(VecGetArray(F,lf_v,ierr))
431: ! Compute function
433: PetscCallA(ApplicationFunction(lx_v,lf_v,ierr))
435: ! Restore vectors
437: PetscCallA(VecRestoreArrayRead(X,lx_v,ierr))
438: PetscCallA(VecRestoreArray(F,lf_v,ierr))
440: PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
441: !
442: ! fdcoloring is in the common block and used here ONLY to test the
443: ! calls to MatFDColoringGetPerturbedColumns() and MatFDColoringRestorePerturbedColumns()
444: !
445: if (fd_coloring) then
446: PetscCallA(MatFDColoringGetPerturbedColumns(fdcoloring,PETSC_NULL_INTEGER,indices,ierr))
447: print*,'Indices from GetPerturbedColumns'
448: write(*,1000) indices
449: 1000 format(50i4)
450: PetscCallA(MatFDColoringRestorePerturbedColumns(fdcoloring,PETSC_NULL_INTEGER,indices,ierr))
451: endif
452: end
454: ! ---------------------------------------------------------------------
455: !
456: ! ApplicationFunction - Computes nonlinear function, called by
457: ! the higher level routine FormFunction().
458: !
459: ! Input Parameter:
460: ! x - local vector data
461: !
462: ! Output Parameters:
463: ! f - local vector data, f(x)
464: ! ierr - error code
465: !
466: ! Notes:
467: ! This routine uses standard Fortran-style computations over a 2-dim array.
468: !
469: subroutine ApplicationFunction(x,f,ierr)
470: use petscsnes
471: implicit none
473: ! Common blocks:
474: PetscReal lambda
475: PetscInt mx,my
476: PetscBool fd_coloring
477: common /params/ lambda,mx,my,fd_coloring
479: ! Input/output variables:
480: PetscScalar x(mx,my),f(mx,my)
481: PetscErrorCode ierr
483: ! Local variables:
484: PetscScalar two,one,hx,hy
485: PetscScalar hxdhy,hydhx,sc
486: PetscScalar u,uxx,uyy
487: PetscInt i,j
489: ierr = 0
490: one = 1.0
491: two = 2.0
492: hx = one/(mx-1)
493: hy = one/(my-1)
494: sc = hx*hy*lambda
495: hxdhy = hx/hy
496: hydhx = hy/hx
498: ! Compute function
500: do 20 j=1,my
501: do 10 i=1,mx
502: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
503: f(i,j) = x(i,j)
504: else
505: u = x(i,j)
506: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
507: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
508: f(i,j) = uxx + uyy - sc*exp(u)
509: endif
510: 10 continue
511: 20 continue
513: end
515: ! ---------------------------------------------------------------------
516: !
517: ! FormJacobian - Evaluates Jacobian matrix.
518: !
519: ! Input Parameters:
520: ! snes - the SNES context
521: ! x - input vector
522: ! dummy - optional user-defined context, as set by SNESSetJacobian()
523: ! (not used here)
524: !
525: ! Output Parameters:
526: ! jac - Jacobian matrix
527: ! jac_prec - optionally different preconditioning matrix (not used here)
528: !
529: ! Notes:
530: ! This routine serves as a wrapper for the lower-level routine
531: ! "ApplicationJacobian", where the actual computations are
532: ! done using the standard Fortran style of treating the local
533: ! vector data as a multidimensional array over the local mesh.
534: ! This routine merely accesses the local vector data via
535: ! VecGetArray() and VecRestoreArray().
536: !
537: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
538: use petscsnes
539: implicit none
541: ! Input/output variables:
542: SNES snes
543: Vec X
544: Mat jac,jac_prec
545: PetscErrorCode ierr
546: integer dummy
548: ! Common blocks:
549: PetscReal lambda
550: PetscInt mx,my
551: PetscBool fd_coloring
552: common /params/ lambda,mx,my,fd_coloring
554: ! Declarations for use with local array:
555: PetscScalar,pointer :: lx_v(:)
557: ! Get a pointer to vector data
559: PetscCallA(VecGetArrayRead(X,lx_v,ierr))
561: ! Compute Jacobian entries
563: PetscCallA(ApplicationJacobian(lx_v,jac,jac_prec,ierr))
565: ! Restore vector
567: PetscCallA(VecRestoreArrayRead(X,lx_v,ierr))
569: ! Assemble matrix
571: PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
572: PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
574: end
576: ! ---------------------------------------------------------------------
577: !
578: ! ApplicationJacobian - Computes Jacobian matrix, called by
579: ! the higher level routine FormJacobian().
580: !
581: ! Input Parameters:
582: ! x - local vector data
583: !
584: ! Output Parameters:
585: ! jac - Jacobian matrix
586: ! jac_prec - optionally different preconditioning matrix (not used here)
587: ! ierr - error code
588: !
589: ! Notes:
590: ! This routine uses standard Fortran-style computations over a 2-dim array.
591: !
592: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
593: use petscsnes
594: implicit none
596: ! Common blocks:
597: PetscReal lambda
598: PetscInt mx,my
599: PetscBool fd_coloring
600: common /params/ lambda,mx,my,fd_coloring
602: ! Input/output variables:
603: PetscScalar x(mx,my)
604: Mat jac,jac_prec
605: PetscErrorCode ierr
607: ! Local variables:
608: PetscInt i,j,row(1),col(5),i1,i5
609: PetscScalar two,one, hx,hy
610: PetscScalar hxdhy,hydhx,sc,v(5)
612: ! Set parameters
614: i1 = 1
615: i5 = 5
616: one = 1.0
617: two = 2.0
618: hx = one/(mx-1)
619: hy = one/(my-1)
620: sc = hx*hy
621: hxdhy = hx/hy
622: hydhx = hy/hx
624: ! Compute entries of the Jacobian matrix
625: ! - Here, we set all entries for a particular row at once.
626: ! - Note that MatSetValues() uses 0-based row and column numbers
627: ! in Fortran as well as in C.
629: do 20 j=1,my
630: row(1) = (j-1)*mx - 1
631: do 10 i=1,mx
632: row(1) = row(1) + 1
633: ! boundary points
634: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
635: PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,[one],INSERT_VALUES,ierr))
636: ! interior grid points
637: else
638: v(1) = -hxdhy
639: v(2) = -hydhx
640: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
641: v(4) = -hydhx
642: v(5) = -hxdhy
643: col(1) = row(1) - mx
644: col(2) = row(1) - 1
645: col(3) = row(1)
646: col(4) = row(1) + 1
647: col(5) = row(1) + mx
648: PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
649: endif
650: 10 continue
651: 20 continue
653: end
655: !
656: !/*TEST
657: !
658: ! build:
659: ! requires: !single !complex
660: !
661: ! test:
662: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
663: !
664: ! test:
665: ! suffix: 2
666: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
667: !
668: ! test:
669: ! suffix: 3
670: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
671: ! filter: sort -b
672: ! filter_output: sort -b
673: !
674: ! test:
675: ! suffix: 4
676: ! args: -pc -par 6.807 -nox
677: !
678: !TEST*/