Actual source code: ex1f.F
petsc-3.7.5 2017-01-01
1: !
2: !
3: ! Description: Uses the Newton method to solve a two-variable system.
4: !
5: !/*T
6: ! Concepts: SNES^basic uniprocessor example
7: ! Processors: 1
8: !T*/
9: !
10: ! -----------------------------------------------------------------------
12: program main
13: implicit none
15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16: ! Include files
17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18: !
19: ! The following include statements are generally used in SNES Fortran
20: ! programs:
21: ! petscsys.h - base PETSc routines
22: ! petscvec.h - vectors
23: ! petscmat.h - matrices
24: ! petscksp.h - Krylov subspace methods
25: ! petscpc.h - preconditioners
26: ! petscsnes.h - SNES interface
27: ! Other include statements may be needed if using additional PETSc
28: ! routines in a Fortran program, e.g.,
29: ! petscviewer.h - viewers
30: ! petscis.h - index sets
31: !
32: #include <petsc/finclude/petscsys.h>
33: #include <petsc/finclude/petscvec.h>
34: #include <petsc/finclude/petscmat.h>
35: #include <petsc/finclude/petscksp.h>
36: #include <petsc/finclude/petscpc.h>
37: #include <petsc/finclude/petscsnes.h>
38: !
39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: ! Variable declarations
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: !
43: ! Variables:
44: ! snes - nonlinear solver
45: ! ksp - linear solver
46: ! pc - preconditioner context
47: ! ksp - Krylov subspace method context
48: ! x, r - solution, residual vectors
49: ! J - Jacobian matrix
50: ! its - iterations for convergence
51: !
52: SNES snes
53: PC pc
54: KSP ksp
55: Vec x,r
56: Mat J
57: SNESLineSearch linesearch
58: PetscErrorCode ierr
59: PetscInt its,i2,i20
60: PetscMPIInt size,rank
61: PetscScalar pfive
62: PetscReal tol
63: PetscBool setls
65: ! Note: Any user-defined Fortran routines (such as FormJacobian)
66: ! MUST be declared as external.
68: external FormFunction, FormJacobian, MyLineSearch
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: ! Macro definitions
72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: !
74: ! Macros to make clearer the process of setting values in vectors and
75: ! getting values from vectors. These vectors are used in the routines
76: ! FormFunction() and FormJacobian().
77: ! - The element lx_a(ib) is element ib in the vector x
78: !
79: #define lx_a(ib) lx_v(lx_i + (ib))
80: #define lf_a(ib) lf_v(lf_i + (ib))
81: !
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! Beginning of program
84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
87: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
88: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
89: if (size .ne. 1) then
90: if (rank .eq. 0) then
91: write(6,*) 'This is a uniprocessor example only!'
92: endif
93: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
94: endif
96: i2 = 2
97: i20 = 20
98: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
99: ! Create nonlinear solver context
100: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
102: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: ! Create matrix and vector data structures; set corresponding routines
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: ! Create vectors for solution and nonlinear function
110: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
111: call VecDuplicate(x,r,ierr)
113: ! Create Jacobian matrix data structure
115: call MatCreate(PETSC_COMM_SELF,J,ierr)
116: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
117: call MatSetFromOptions(J,ierr)
118: call MatSetUp(J,ierr)
120: ! Set function evaluation routine and vector
122: call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr)
124: ! Set Jacobian matrix data structure and Jacobian evaluation routine
126: call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT, &
127: & ierr)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: ! Customize nonlinear solver; set runtime options
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! Set linear solver defaults for this problem. By extracting the
134: ! KSP, KSP, and PC contexts from the SNES context, we can then
135: ! directly call any KSP, KSP, and PC routines to set various options.
137: call SNESGetKSP(snes,ksp,ierr)
138: call KSPGetPC(ksp,pc,ierr)
139: call PCSetType(pc,PCNONE,ierr)
140: tol = 1.e-4
141: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, &
142: & PETSC_DEFAULT_REAL,i20,ierr)
144: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
145: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
146: ! These options will override those specified above as long as
147: ! SNESSetFromOptions() is called _after_ any other customization
148: ! routines.
151: call SNESSetFromOptions(snes,ierr)
153: call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
154: & '-setls',setls,ierr)
156: if (setls) then
157: call SNESGetLineSearch(snes, linesearch, ierr)
158: call SNESLineSearchSetType(linesearch, 'shell', ierr)
159: call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch, &
160: &PETSC_NULL_OBJECT, ierr)
161: endif
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Evaluate initial guess; then solve nonlinear system
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: ! Note: The user should initialize the vector, x, with the initial guess
168: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
169: ! to employ an initial guess of zero, the user should explicitly set
170: ! this vector to zero by calling VecSet().
172: pfive = 0.5
173: call VecSet(x,pfive,ierr)
174: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
175: call SNESGetIterationNumber(snes,its,ierr);
176: if (rank .eq. 0) then
177: write(6,100) its
178: endif
179: 100 format('Number of SNES iterations = ',i5)
181: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: ! Free work space. All PETSc objects should be destroyed when they
183: ! are no longer needed.
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: call VecDestroy(x,ierr)
187: call VecDestroy(r,ierr)
188: call MatDestroy(J,ierr)
189: call SNESDestroy(snes,ierr)
190: call PetscFinalize(ierr)
191: end
192: !
193: ! ------------------------------------------------------------------------
194: !
195: ! FormFunction - Evaluates nonlinear function, F(x).
196: !
197: ! Input Parameters:
198: ! snes - the SNES context
199: ! x - input vector
200: ! dummy - optional user-defined context (not used here)
201: !
202: ! Output Parameter:
203: ! f - function vector
204: !
205: subroutine FormFunction(snes,x,f,dummy,ierr)
206: implicit none
208: #include <petsc/finclude/petscsys.h>
209: #include <petsc/finclude/petscvec.h>
210: #include <petsc/finclude/petscsnes.h>
212: SNES snes
213: Vec x,f
214: PetscErrorCode ierr
215: integer dummy(*)
217: ! Declarations for use with local arrays
219: PetscScalar lx_v(2),lf_v(2)
220: PetscOffset lx_i,lf_i
222: ! Get pointers to vector data.
223: ! - For default PETSc vectors, VecGetArray() returns a pointer to
224: ! the data array. Otherwise, the routine is implementation dependent.
225: ! - You MUST call VecRestoreArray() when you no longer need access to
226: ! the array.
227: ! - Note that the Fortran interface to VecGetArray() differs from the
228: ! C version. See the Fortran chapter of the users manual for details.
230: call VecGetArrayRead(x,lx_v,lx_i,ierr)
231: call VecGetArray(f,lf_v,lf_i,ierr)
233: ! Compute function
235: lf_a(1) = lx_a(1)*lx_a(1) &
236: & + lx_a(1)*lx_a(2) - 3.0
237: lf_a(2) = lx_a(1)*lx_a(2) &
238: & + lx_a(2)*lx_a(2) - 6.0
240: ! Restore vectors
242: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
243: call VecRestoreArray(f,lf_v,lf_i,ierr)
245: return
246: end
248: ! ---------------------------------------------------------------------
249: !
250: ! FormJacobian - Evaluates Jacobian matrix.
251: !
252: ! Input Parameters:
253: ! snes - the SNES context
254: ! x - input vector
255: ! dummy - optional user-defined context (not used here)
256: !
257: ! Output Parameters:
258: ! A - Jacobian matrix
259: ! B - optionally different preconditioning matrix
260: ! flag - flag indicating matrix structure
261: !
262: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
263: implicit none
265: #include <petsc/finclude/petscsys.h>
266: #include <petsc/finclude/petscvec.h>
267: #include <petsc/finclude/petscmat.h>
268: #include <petsc/finclude/petscpc.h>
269: #include <petsc/finclude/petscsnes.h>
271: SNES snes
272: Vec X
273: Mat jac,B
274: PetscScalar A(4)
275: PetscErrorCode ierr
276: PetscInt idx(2),i2
277: integer dummy(*)
279: ! Declarations for use with local arrays
281: PetscScalar lx_v(2)
282: PetscOffset lx_i
284: ! Get pointer to vector data
286: i2 = 2
287: call VecGetArrayRead(x,lx_v,lx_i,ierr)
289: ! Compute Jacobian entries and insert into matrix.
290: ! - Since this is such a small problem, we set all entries for
291: ! the matrix at once.
292: ! - Note that MatSetValues() uses 0-based row and column numbers
293: ! in Fortran as well as in C (as set here in the array idx).
295: idx(1) = 0
296: idx(2) = 1
297: A(1) = 2.0*lx_a(1) + lx_a(2)
298: A(2) = lx_a(1)
299: A(3) = lx_a(2)
300: A(4) = lx_a(1) + 2.0*lx_a(2)
301: call MatSetValues(jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
303: ! Restore vector
305: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
307: ! Assemble matrix
309: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
310: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
312: return
313: end
316: subroutine MyLineSearch(linesearch, lctx, ierr)
317: #include <petsc/finclude/petscsys.h>
318: #include <petsc/finclude/petscvec.h>
319: #include <petsc/finclude/petscmat.h>
320: #include <petsc/finclude/petscksp.h>
321: #include <petsc/finclude/petscpc.h>
322: #include <petsc/finclude/petscsnes.h>
324: SNES snes
325: integer lctx
326: Vec x, f,g, y, w
327: PetscReal ynorm,gnorm,xnorm
328: PetscBool flag
329: PetscErrorCode ierr
331: PetscScalar mone
333: mone = -1.0
334: call SNESLineSearchGetSNES(linesearch, snes, ierr)
335: call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
336: call VecNorm(y,NORM_2,ynorm,ierr)
337: call VecAXPY(x,mone,y,ierr)
338: call SNESComputeFunction(snes,x,f,ierr)
339: call VecNorm(f,NORM_2,gnorm,ierr)
340: call VecNorm(x,NORM_2,xnorm,ierr)
341: call VecNorm(y,NORM_2,ynorm,ierr)
342: call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, &
343: & ierr)
344: flag = PETSC_FALSE
345: return
346: end