Actual source code: ex1f.F90
1: !
2: !
3: ! Description: Uses the Newton method to solve a two-variable system.
4: !
6: program main
7: #include <petsc/finclude/petsc.h>
8: use petsc
9: implicit none
11: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
12: ! Variable declarations
13: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
14: !
15: ! Variables:
16: ! snes - nonlinear solver
17: ! ksp - linear solver
18: ! pc - preconditioner context
19: ! ksp - Krylov subspace method context
20: ! x, r - solution, residual vectors
21: ! J - Jacobian matrix
22: ! its - iterations for convergence
23: !
24: SNES snes
25: PC pc
26: KSP ksp
27: Vec x,r
28: Mat J
29: SNESLineSearch linesearch
30: PetscErrorCode ierr
31: PetscInt its,i2,i20
32: PetscMPIInt size,rank
33: PetscScalar pfive
34: PetscReal tol
35: PetscBool setls
36: #if defined(PETSC_USE_LOG)
37: PetscViewer viewer
38: #endif
39: double precision threshold,oldthreshold
41: ! Note: Any user-defined Fortran routines (such as FormJacobian)
42: ! MUST be declared as external.
44: external FormFunction, FormJacobian, MyLineSearch
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: ! Macro definitions
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: !
50: ! Macros to make clearer the process of setting values in vectors and
51: ! getting values from vectors. These vectors are used in the routines
52: ! FormFunction() and FormJacobian().
53: ! - The element lx_a(ib) is element ib in the vector x
54: !
55: #define lx_a(ib) lx_v(lx_i + (ib))
56: #define lf_a(ib) lf_v(lf_i + (ib))
57: !
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: ! Beginning of program
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
63: if (ierr .ne. 0) then
64: print*,'Unable to initialize PETSc'
65: stop
66: endif
67: call PetscLogNestedBegin(ierr);CHKERRA(ierr)
68: threshold = 1.0
69: call PetscLogSetThreshold(threshold,oldthreshold,ierr)
70: ! dummy test of logging a reduction
71: #if defined(PETSC_USE_LOG)
72: PetscAReduce()
73: #endif
74: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
75: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
76: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif
78: i2 = 2
79: i20 = 20
80: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
81: ! Create nonlinear solver context
82: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
84: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
86: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: ! Create matrix and vector data structures; set corresponding routines
88: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: ! Create vectors for solution and nonlinear function
92: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
93: call VecDuplicate(x,r,ierr)
95: ! Create Jacobian matrix data structure
97: call MatCreate(PETSC_COMM_SELF,J,ierr)
98: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
99: call MatSetFromOptions(J,ierr)
100: call MatSetUp(J,ierr)
102: ! Set function evaluation routine and vector
104: call SNESSetFunction(snes,r,FormFunction,0,ierr)
106: ! Set Jacobian matrix data structure and Jacobian evaluation routine
108: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Customize nonlinear solver; set runtime options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! Set linear solver defaults for this problem. By extracting the
115: ! KSP, KSP, and PC contexts from the SNES context, we can then
116: ! directly call any KSP, KSP, and PC routines to set various options.
118: call SNESGetKSP(snes,ksp,ierr)
119: call KSPGetPC(ksp,pc,ierr)
120: call PCSetType(pc,PCNONE,ierr)
121: tol = 1.e-4
122: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
123: & PETSC_DEFAULT_REAL,i20,ierr)
125: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
126: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
127: ! These options will override those specified above as long as
128: ! SNESSetFromOptions() is called _after_ any other customization
129: ! routines.
131: call SNESSetFromOptions(snes,ierr)
133: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
134: & '-setls',setls,ierr)
136: if (setls) then
137: call SNESGetLineSearch(snes, linesearch, ierr)
138: call SNESLineSearchSetType(linesearch, 'shell', ierr)
139: call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch, &
140: & 0, ierr)
141: endif
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! Evaluate initial guess; then solve nonlinear system
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Note: The user should initialize the vector, x, with the initial guess
148: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
149: ! to employ an initial guess of zero, the user should explicitly set
150: ! this vector to zero by calling VecSet().
152: pfive = 0.5
153: call VecSet(x,pfive,ierr)
154: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
156: ! View solver converged reason; we could instead use the option -snes_converged_reason
157: call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr)
159: call SNESGetIterationNumber(snes,its,ierr);
160: if (rank .eq. 0) then
161: write(6,100) its
162: endif
163: 100 format('Number of SNES iterations = ',i5)
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: ! Free work space. All PETSc objects should be destroyed when they
167: ! are no longer needed.
168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: call VecDestroy(x,ierr)
171: call VecDestroy(r,ierr)
172: call MatDestroy(J,ierr)
173: call SNESDestroy(snes,ierr)
174: #if defined(PETSC_USE_LOG)
175: call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)
176: call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)
177: call PetscLogView(viewer,ierr)
178: call PetscViewerDestroy(viewer,ierr)
179: #endif
180: call PetscFinalize(ierr)
181: end
182: !
183: ! ------------------------------------------------------------------------
184: !
185: ! FormFunction - Evaluates nonlinear function, F(x).
186: !
187: ! Input Parameters:
188: ! snes - the SNES context
189: ! x - input vector
190: ! dummy - optional user-defined context (not used here)
191: !
192: ! Output Parameter:
193: ! f - function vector
194: !
195: subroutine FormFunction(snes,x,f,dummy,ierr)
196: use petscsnes
197: implicit none
199: SNES snes
200: Vec x,f
201: PetscErrorCode ierr
202: integer dummy(*)
204: ! Declarations for use with local arrays
206: PetscScalar lx_v(2),lf_v(2)
207: PetscOffset lx_i,lf_i
209: ! Get pointers to vector data.
210: ! - For default PETSc vectors, VecGetArray() returns a pointer to
211: ! the data array. Otherwise, the routine is implementation dependent.
212: ! - You MUST call VecRestoreArray() when you no longer need access to
213: ! the array.
214: ! - Note that the Fortran interface to VecGetArray() differs from the
215: ! C version. See the Fortran chapter of the users manual for details.
217: call VecGetArrayRead(x,lx_v,lx_i,ierr)
218: call VecGetArray(f,lf_v,lf_i,ierr)
220: ! Compute function
222: lf_a(1) = lx_a(1)*lx_a(1) &
223: & + lx_a(1)*lx_a(2) - 3.0
224: lf_a(2) = lx_a(1)*lx_a(2) &
225: & + lx_a(2)*lx_a(2) - 6.0
227: ! Restore vectors
229: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
230: call VecRestoreArray(f,lf_v,lf_i,ierr)
232: return
233: end
235: ! ---------------------------------------------------------------------
236: !
237: ! FormJacobian - Evaluates Jacobian matrix.
238: !
239: ! Input Parameters:
240: ! snes - the SNES context
241: ! x - input vector
242: ! dummy - optional user-defined context (not used here)
243: !
244: ! Output Parameters:
245: ! A - Jacobian matrix
246: ! B - optionally different preconditioning matrix
247: !
248: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
249: use petscsnes
250: implicit none
252: SNES snes
253: Vec X
254: Mat jac,B
255: PetscScalar A(4)
256: PetscErrorCode ierr
257: PetscInt idx(2),i2
258: integer dummy(*)
260: ! Declarations for use with local arrays
262: PetscScalar lx_v(2)
263: PetscOffset lx_i
265: ! Get pointer to vector data
267: i2 = 2
268: call VecGetArrayRead(x,lx_v,lx_i,ierr)
270: ! Compute Jacobian entries and insert into matrix.
271: ! - Since this is such a small problem, we set all entries for
272: ! the matrix at once.
273: ! - Note that MatSetValues() uses 0-based row and column numbers
274: ! in Fortran as well as in C (as set here in the array idx).
276: idx(1) = 0
277: idx(2) = 1
278: A(1) = 2.0*lx_a(1) + lx_a(2)
279: A(2) = lx_a(1)
280: A(3) = lx_a(2)
281: A(4) = lx_a(1) + 2.0*lx_a(2)
282: call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
284: ! Restore vector
286: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
288: ! Assemble matrix
290: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
291: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
292: if (B .ne. jac) then
293: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
294: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
295: endif
297: return
298: end
300: subroutine MyLineSearch(linesearch, lctx, ierr)
301: use petscsnes
302: implicit none
304: SNESLineSearch linesearch
305: SNES snes
306: integer lctx
307: Vec x, f,g, y, w
308: PetscReal ynorm,gnorm,xnorm
309: PetscErrorCode ierr
311: PetscScalar mone
313: mone = -1.0
314: call SNESLineSearchGetSNES(linesearch, snes, ierr)
315: call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
316: call VecNorm(y,NORM_2,ynorm,ierr)
317: call VecAXPY(x,mone,y,ierr)
318: call SNESComputeFunction(snes,x,f,ierr)
319: call VecNorm(f,NORM_2,gnorm,ierr)
320: call VecNorm(x,NORM_2,xnorm,ierr)
321: call VecNorm(y,NORM_2,ynorm,ierr)
322: call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, &
323: & ierr)
324: return
325: end
327: !/*TEST
328: !
329: ! test:
330: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
331: ! requires: !single
332: !
333: !TEST*/