Actual source code: ex1f.F90
1: !
2: ! Description: This example solves a nonlinear system on 1 processor with SNES.
3: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
4: ! domain. The command line options include:
5: ! -par <parameter>, where <parameter> indicates the nonlinearity of the problem
6: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
7: ! -mx <xg>, where <xg> = number of grid points in the x-direction
8: ! -my <yg>, where <yg> = number of grid points in the y-direction
9: !
11: !
12: ! --------------------------------------------------------------------------
13: !
14: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
15: ! the partial differential equation
16: !
17: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
18: !
19: ! with boundary conditions
20: !
21: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
22: !
23: ! A finite difference approximation with the usual 5-point stencil
24: ! is used to discretize the boundary value problem to obtain a nonlinear
25: ! system of equations.
26: !
27: ! The parallel version of this code is snes/tutorials/ex5f.F
28: !
29: ! --------------------------------------------------------------------------
30: subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr)
31: #include <petsc/finclude/petscsnes.h>
32: use petscsnes
33: implicit none
34: SNES snes
35: PetscReal norm
36: Vec tmp,x,y,w
37: PetscBool changed_w,changed_y
38: PetscErrorCode ierr
39: PetscInt ctx
40: PetscScalar mone
41: MPI_Comm comm
43: character(len=PETSC_MAX_PATH_LEN) :: outputString
45: PetscCallA(PetscObjectGetComm(snes,comm,ierr))
46: PetscCallA(VecDuplicate(x,tmp,ierr))
47: mone = -1.0
48: PetscCallA(VecWAXPY(tmp,mone,x,w,ierr))
49: PetscCallA(VecNorm(tmp,NORM_2,norm,ierr))
50: PetscCallA(VecDestroy(tmp,ierr))
51: write(outputString,*) norm
52: PetscCallA(PetscPrintf(comm,'Norm of search step '//trim(outputString)//'\n',ierr))
53: end
55: program main
56: #include <petsc/finclude/petscdraw.h>
57: use petscsnes
58: implicit none
59: interface SNESSetJacobian
60: subroutine SNESSetJacobian1(a,b,c,d,e,z)
61: use petscsnes
62: SNES a
63: Mat b
64: Mat c
65: external d
66: MatFDColoring e
67: PetscErrorCode z
68: end subroutine
69: subroutine SNESSetJacobian2(a,b,c,d,e,z)
70: use petscsnes
71: SNES a
72: Mat b
73: Mat c
74: external d
75: integer e
76: PetscErrorCode z
77: end subroutine
78: end interface
79: !
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: ! Variable declarations
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: !
84: ! Variables:
85: ! snes - nonlinear solver
86: ! x, r - solution, residual vectors
87: ! J - Jacobian matrix
88: ! its - iterations for convergence
89: ! matrix_free - flag - 1 indicates matrix-free version
90: ! lambda - nonlinearity parameter
91: ! draw - drawing context
92: !
93: SNES snes
94: MatColoring mc
95: Vec x,r
96: PetscDraw draw
97: Mat J
98: PetscBool matrix_free,flg,fd_coloring
99: PetscErrorCode ierr
100: PetscInt its,N, mx,my,i5
101: PetscMPIInt size,rank
102: PetscReal lambda_max,lambda_min,lambda
103: MatFDColoring fdcoloring
104: ISColoring iscoloring
105: PetscBool pc
106: external postcheck
108: character(len=PETSC_MAX_PATH_LEN) :: outputString
110: PetscScalar,pointer :: lx_v(:)
112: ! Store parameters in common block
114: common /params/ lambda,mx,my,fd_coloring
116: ! Note: Any user-defined Fortran routines (such as FormJacobian)
117: ! MUST be declared as external.
119: external FormFunction,FormInitialGuess,FormJacobian
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! Initialize program
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: PetscCallA(PetscInitialize(ierr))
126: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
127: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
129: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')
131: ! Initialize problem parameters
132: i5 = 5
133: lambda_max = 6.81
134: lambda_min = 0.0
135: lambda = 6.0
136: mx = 4
137: my = 4
138: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
139: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
140: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr))
141: PetscCheckA(lambda .lt. lambda_max .and. lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range ')
142: N = mx*my
143: pc = PETSC_FALSE
144: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr))
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Create nonlinear solver context
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
152: if (pc .eqv. PETSC_TRUE) then
153: PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr))
154: PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr))
155: endif
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: ! Create vector data structures; set function evaluation routine
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
162: PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr))
163: PetscCallA(VecSetFromOptions(x,ierr))
164: PetscCallA(VecDuplicate(x,r,ierr))
166: ! Set function evaluation routine and vector. Whenever the nonlinear
167: ! solver needs to evaluate the nonlinear function, it will call this
168: ! routine.
169: ! - Note that the final routine argument is the user-defined
170: ! context that provides application-specific data for the
171: ! function evaluation routine.
173: PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr))
175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: ! Create matrix data structure; set Jacobian evaluation routine
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: ! Create matrix. Here we only approximately preallocate storage space
180: ! for the Jacobian. See the users manual for a discussion of better
181: ! techniques for preallocating matrix memory.
183: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr))
184: if (.not. matrix_free) then
185: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr))
186: endif
188: !
189: ! This option will cause the Jacobian to be computed via finite differences
190: ! efficiently using a coloring of the columns of the matrix.
191: !
192: fd_coloring = .false.
193: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr))
194: if (fd_coloring) then
196: !
197: ! This initializes the nonzero structure of the Jacobian. This is artificial
198: ! because clearly if we had a routine to compute the Jacobian we won't need
199: ! to use finite differences.
200: !
201: PetscCallA(FormJacobian(snes,x,J,J,0,ierr))
202: !
203: ! Color the matrix, i.e. determine groups of columns that share no common
204: ! rows. These columns in the Jacobian can all be computed simultaneously.
205: !
206: PetscCallA(MatColoringCreate(J,mc,ierr))
207: PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr))
208: PetscCallA(MatColoringSetFromOptions(mc,ierr))
209: PetscCallA(MatColoringApply(mc,iscoloring,ierr))
210: PetscCallA(MatColoringDestroy(mc,ierr))
211: !
212: ! Create the data structure that SNESComputeJacobianDefaultColor() uses
213: ! to compute the actual Jacobians via finite differences.
214: !
215: PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr))
216: PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr))
217: PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr))
218: PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr))
219: !
220: ! Tell SNES to use the routine SNESComputeJacobianDefaultColor()
221: ! to compute Jacobians.
222: !
223: PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr))
224: PetscCallA(ISColoringDestroy(iscoloring,ierr))
226: else if (.not. matrix_free) then
228: ! Set Jacobian matrix data structure and default Jacobian evaluation
229: ! routine. Whenever the nonlinear solver needs to compute the
230: ! Jacobian matrix, it will call this routine.
231: ! - Note that the final routine argument is the user-defined
232: ! context that provides application-specific data for the
233: ! Jacobian evaluation routine.
234: ! - The user can override with:
235: ! -snes_fd : default finite differencing approximation of Jacobian
236: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
237: ! (unless user explicitly sets preconditioner)
238: ! -snes_mf_operator : form preconditioning matrix as set by the user,
239: ! but use matrix-free approx for Jacobian-vector
240: ! products within Newton-Krylov method
241: !
242: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
243: endif
245: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
246: ! Customize nonlinear solver; set runtime options
247: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
251: PetscCallA(SNESSetFromOptions(snes,ierr))
253: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254: ! Evaluate initial guess; then solve nonlinear system.
255: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257: ! Note: The user should initialize the vector, x, with the initial guess
258: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
259: ! to employ an initial guess of zero, the user should explicitly set
260: ! this vector to zero by calling VecSet().
262: PetscCallA(FormInitialGuess(x,ierr))
263: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
264: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
265: write(outputString,*) its
266: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,'Number of SNES iterations = '//trim(outputString)//'\n',ierr))
268: ! PetscDraw contour plot of solution
270: PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr))
271: PetscCallA(PetscDrawSetFromOptions(draw,ierr))
273: PetscCallA(VecGetArrayReadF90(x,lx_v,ierr))
274: PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v,ierr))
275: PetscCallA(VecRestoreArrayReadF90(x,lx_v,ierr))
277: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: ! Free work space. All PETSc objects should be destroyed when they
279: ! are no longer needed.
280: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
282: if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr))
283: if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr))
285: PetscCallA(VecDestroy(x,ierr))
286: PetscCallA(VecDestroy(r,ierr))
287: PetscCallA(SNESDestroy(snes,ierr))
288: PetscCallA(PetscDrawDestroy(draw,ierr))
289: PetscCallA(PetscFinalize(ierr))
290: end
292: ! ---------------------------------------------------------------------
293: !
294: ! FormInitialGuess - Forms initial approximation.
295: !
296: ! Input Parameter:
297: ! X - vector
298: !
299: ! Output Parameters:
300: ! X - vector
301: ! ierr - error code
302: !
303: ! Notes:
304: ! This routine serves as a wrapper for the lower-level routine
305: ! "ApplicationInitialGuess", where the actual computations are
306: ! done using the standard Fortran style of treating the local
307: ! vector data as a multidimensional array over the local mesh.
308: ! This routine merely accesses the local vector data via
309: ! VecGetArrayF90() and VecRestoreArrayF90().
310: !
311: subroutine FormInitialGuess(X,ierr)
312: use petscsnes
313: implicit none
315: ! Input/output variables:
316: Vec X
317: PetscErrorCode ierr
319: ! Declarations for use with local arrays:
320: PetscScalar,pointer :: lx_v(:)
322: ierr = 0
324: ! Get a pointer to vector data.
325: ! - VecGetArrayF90() returns a pointer to the data array.
326: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
327: ! the array.
329: PetscCallA(VecGetArrayF90(X,lx_v,ierr))
331: ! Compute initial guess
333: PetscCallA(ApplicationInitialGuess(lx_v,ierr))
335: ! Restore vector
337: PetscCallA(VecRestoreArrayF90(X,lx_v,ierr))
339: end
341: ! ApplicationInitialGuess - Computes initial approximation, called by
342: ! the higher level routine FormInitialGuess().
343: !
344: ! Input Parameter:
345: ! x - local vector data
346: !
347: ! Output Parameters:
348: ! f - local vector data, f(x)
349: ! ierr - error code
350: !
351: ! Notes:
352: ! This routine uses standard Fortran-style computations over a 2-dim array.
353: !
354: subroutine ApplicationInitialGuess(x,ierr)
355: use petscksp
356: implicit none
358: ! Common blocks:
359: PetscReal lambda
360: PetscInt mx,my
361: PetscBool fd_coloring
362: common /params/ lambda,mx,my,fd_coloring
364: ! Input/output variables:
365: PetscScalar x(mx,my)
366: PetscErrorCode ierr
368: ! Local variables:
369: PetscInt i,j
370: PetscReal temp1,temp,hx,hy,one
372: ! Set parameters
374: ierr = 0
375: one = 1.0
376: hx = one/(mx-1)
377: hy = one/(my-1)
378: temp1 = lambda/(lambda + one)
380: do 20 j=1,my
381: temp = min(j-1,my-j)*hy
382: do 10 i=1,mx
383: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
384: x(i,j) = 0.0
385: else
386: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp))
387: endif
388: 10 continue
389: 20 continue
391: end
393: ! ---------------------------------------------------------------------
394: !
395: ! FormFunction - Evaluates nonlinear function, F(x).
396: !
397: ! Input Parameters:
398: ! snes - the SNES context
399: ! X - input vector
400: ! dummy - optional user-defined context, as set by SNESSetFunction()
401: ! (not used here)
402: !
403: ! Output Parameter:
404: ! F - vector with newly computed function
405: !
406: ! Notes:
407: ! This routine serves as a wrapper for the lower-level routine
408: ! "ApplicationFunction", where the actual computations are
409: ! done using the standard Fortran style of treating the local
410: ! vector data as a multidimensional array over the local mesh.
411: ! This routine merely accesses the local vector data via
412: ! VecGetArrayF90() and VecRestoreArrayF90().
413: !
414: subroutine FormFunction(snes,X,F,fdcoloring,ierr)
415: use petscsnes
416: implicit none
418: ! Input/output variables:
419: SNES snes
420: Vec X,F
421: PetscErrorCode ierr
422: MatFDColoring fdcoloring
424: ! Common blocks:
425: PetscReal lambda
426: PetscInt mx,my
427: PetscBool fd_coloring
428: common /params/ lambda,mx,my,fd_coloring
430: ! Declarations for use with local arrays:
431: PetscScalar,pointer :: lx_v(:), lf_v(:)
432: PetscInt, pointer :: indices(:)
434: ! Get pointers to vector data.
435: ! - VecGetArrayF90() returns a pointer to the data array.
436: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
437: ! the array.
439: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
440: PetscCallA(VecGetArrayF90(F,lf_v,ierr))
442: ! Compute function
444: PetscCallA(ApplicationFunction(lx_v,lf_v,ierr))
446: ! Restore vectors
448: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
449: PetscCallA(VecRestoreArrayF90(F,lf_v,ierr))
451: PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr))
452: !
453: ! fdcoloring is in the common block and used here ONLY to test the
454: ! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90()
455: !
456: if (fd_coloring) then
457: PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr))
458: print*,'Indices from GetPerturbedColumnsF90'
459: write(*,1000) indices
460: 1000 format(50i4)
461: PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr))
462: endif
463: end
465: ! ---------------------------------------------------------------------
466: !
467: ! ApplicationFunction - Computes nonlinear function, called by
468: ! the higher level routine FormFunction().
469: !
470: ! Input Parameter:
471: ! x - local vector data
472: !
473: ! Output Parameters:
474: ! f - local vector data, f(x)
475: ! ierr - error code
476: !
477: ! Notes:
478: ! This routine uses standard Fortran-style computations over a 2-dim array.
479: !
480: subroutine ApplicationFunction(x,f,ierr)
481: use petscsnes
482: implicit none
484: ! Common blocks:
485: PetscReal lambda
486: PetscInt mx,my
487: PetscBool fd_coloring
488: common /params/ lambda,mx,my,fd_coloring
490: ! Input/output variables:
491: PetscScalar x(mx,my),f(mx,my)
492: PetscErrorCode ierr
494: ! Local variables:
495: PetscScalar two,one,hx,hy
496: PetscScalar hxdhy,hydhx,sc
497: PetscScalar u,uxx,uyy
498: PetscInt i,j
500: ierr = 0
501: one = 1.0
502: two = 2.0
503: hx = one/(mx-1)
504: hy = one/(my-1)
505: sc = hx*hy*lambda
506: hxdhy = hx/hy
507: hydhx = hy/hx
509: ! Compute function
511: do 20 j=1,my
512: do 10 i=1,mx
513: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
514: f(i,j) = x(i,j)
515: else
516: u = x(i,j)
517: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
518: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
519: f(i,j) = uxx + uyy - sc*exp(u)
520: endif
521: 10 continue
522: 20 continue
524: end
526: ! ---------------------------------------------------------------------
527: !
528: ! FormJacobian - Evaluates Jacobian matrix.
529: !
530: ! Input Parameters:
531: ! snes - the SNES context
532: ! x - input vector
533: ! dummy - optional user-defined context, as set by SNESSetJacobian()
534: ! (not used here)
535: !
536: ! Output Parameters:
537: ! jac - Jacobian matrix
538: ! jac_prec - optionally different preconditioning matrix (not used here)
539: ! flag - flag indicating matrix structure
540: !
541: ! Notes:
542: ! This routine serves as a wrapper for the lower-level routine
543: ! "ApplicationJacobian", where the actual computations are
544: ! done using the standard Fortran style of treating the local
545: ! vector data as a multidimensional array over the local mesh.
546: ! This routine merely accesses the local vector data via
547: ! VecGetArrayF90() and VecRestoreArrayF90().
548: !
549: subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr)
550: use petscsnes
551: implicit none
553: ! Input/output variables:
554: SNES snes
555: Vec X
556: Mat jac,jac_prec
557: PetscErrorCode ierr
558: integer dummy
560: ! Common blocks:
561: PetscReal lambda
562: PetscInt mx,my
563: PetscBool fd_coloring
564: common /params/ lambda,mx,my,fd_coloring
566: ! Declarations for use with local array:
567: PetscScalar,pointer :: lx_v(:)
569: ! Get a pointer to vector data
571: PetscCallA(VecGetArrayReadF90(X,lx_v,ierr))
573: ! Compute Jacobian entries
575: PetscCallA(ApplicationJacobian(lx_v,jac,jac_prec,ierr))
577: ! Restore vector
579: PetscCallA(VecRestoreArrayReadF90(X,lx_v,ierr))
581: ! Assemble matrix
583: PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
584: PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr))
586: end
588: ! ---------------------------------------------------------------------
589: !
590: ! ApplicationJacobian - Computes Jacobian matrix, called by
591: ! the higher level routine FormJacobian().
592: !
593: ! Input Parameters:
594: ! x - local vector data
595: !
596: ! Output Parameters:
597: ! jac - Jacobian matrix
598: ! jac_prec - optionally different preconditioning matrix (not used here)
599: ! ierr - error code
600: !
601: ! Notes:
602: ! This routine uses standard Fortran-style computations over a 2-dim array.
603: !
604: subroutine ApplicationJacobian(x,jac,jac_prec,ierr)
605: use petscsnes
606: implicit none
608: ! Common blocks:
609: PetscReal lambda
610: PetscInt mx,my
611: PetscBool fd_coloring
612: common /params/ lambda,mx,my,fd_coloring
614: ! Input/output variables:
615: PetscScalar x(mx,my)
616: Mat jac,jac_prec
617: PetscErrorCode ierr
619: ! Local variables:
620: PetscInt i,j,row(1),col(5),i1,i5
621: PetscScalar two,one, hx,hy
622: PetscScalar hxdhy,hydhx,sc,v(5)
624: ! Set parameters
626: i1 = 1
627: i5 = 5
628: one = 1.0
629: two = 2.0
630: hx = one/(mx-1)
631: hy = one/(my-1)
632: sc = hx*hy
633: hxdhy = hx/hy
634: hydhx = hy/hx
636: ! Compute entries of the Jacobian matrix
637: ! - Here, we set all entries for a particular row at once.
638: ! - Note that MatSetValues() uses 0-based row and column numbers
639: ! in Fortran as well as in C.
641: do 20 j=1,my
642: row(1) = (j-1)*mx - 1
643: do 10 i=1,mx
644: row(1) = row(1) + 1
645: ! boundary points
646: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
647: PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr))
648: ! interior grid points
649: else
650: v(1) = -hxdhy
651: v(2) = -hydhx
652: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
653: v(4) = -hydhx
654: v(5) = -hxdhy
655: col(1) = row(1) - mx
656: col(2) = row(1) - 1
657: col(3) = row(1)
658: col(4) = row(1) + 1
659: col(5) = row(1) + mx
660: PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr))
661: endif
662: 10 continue
663: 20 continue
665: end
667: !
668: !/*TEST
669: !
670: ! build:
671: ! requires: !single
672: !
673: ! test:
674: ! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
675: !
676: ! test:
677: ! suffix: 2
678: ! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always
679: !
680: ! test:
681: ! suffix: 3
682: ! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always
683: ! filter: sort -b
684: ! filter_output: sort -b
685: !
686: ! test:
687: ! suffix: 4
688: ! args: -pc -par 6.807 -nox
689: !
690: !TEST*/