Actual source code: ex1f.F
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: ! petsc.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 include/finclude/petsc.h
33: #include include/finclude/petscvec.h
34: #include include/finclude/petscmat.h
35: #include include/finclude/petscksp.h
36: #include include/finclude/petscpc.h
37: #include include/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: PetscErrorCode ierr
58: PetscInt its,i2,i20
59: PetscMPIInt size,rank
60: PetscScalar pfive
61: double precision tol
62: PetscTruth setls
64: ! Note: Any user-defined Fortran routines (such as FormJacobian)
65: ! MUST be declared as external.
67: external FormFunction, FormJacobian, MyLineSearch
69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: ! Macro definitions
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: !
73: ! Macros to make clearer the process of setting values in vectors and
74: ! getting values from vectors. These vectors are used in the routines
75: ! FormFunction() and FormJacobian().
76: ! - The element lx_a(ib) is element ib in the vector x
77: !
78: #define lx_a(ib) lx_v(lx_i + (ib))
79: #define lf_a(ib) lf_v(lf_i + (ib))
80: !
81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: ! Beginning of program
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
86: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
87: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
88: if (size .ne. 1) then
89: if (rank .eq. 0) then
90: write(6,*) 'This is a uniprocessor example only!'
91: endif
92: SETERRQ(1,' ',ierr)
93: endif
95: i2 = 2
96: i20 = 20
97: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
98: ! Create nonlinear solver context
99: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
101: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! Create matrix and vector data structures; set corresponding routines
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Create vectors for solution and nonlinear function
109: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
110: call VecDuplicate(x,r,ierr)
112: ! Create Jacobian matrix data structure
114: call MatCreate(PETSC_COMM_SELF,J,ierr)
115: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
116: call MatSetFromOptions(J,ierr)
118: ! Set function evaluation routine and vector
120: call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr)
122: ! Set Jacobian matrix data structure and Jacobian evaluation routine
124: call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT, &
125: & ierr)
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! Customize nonlinear solver; set runtime options
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: ! Set linear solver defaults for this problem. By extracting the
132: ! KSP, KSP, and PC contexts from the SNES context, we can then
133: ! directly call any KSP, KSP, and PC routines to set various options.
135: call SNESGetKSP(snes,ksp,ierr)
136: call KSPGetPC(ksp,pc,ierr)
137: call PCSetType(pc,PCNONE,ierr)
138: tol = 1.e-4
139: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
140: & PETSC_DEFAULT_DOUBLE_PRECISION,i20,ierr)
142: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
143: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
144: ! These options will override those specified above as long as
145: ! SNESSetFromOptions() is called _after_ any other customization
146: ! routines.
149: call SNESSetFromOptions(snes,ierr)
151: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-setls',setls,ierr)
153: if (setls .eq. PETSC_TRUE) then
154: call SNESLineSearchSet(snes,MyLineSearch, &
155: & PETSC_NULL_OBJECT,ierr)
156: endif
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Evaluate initial guess; then solve nonlinear system
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! Note: The user should initialize the vector, x, with the initial guess
163: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
164: ! to employ an initial guess of zero, the user should explicitly set
165: ! this vector to zero by calling VecSet().
167: pfive = 0.5
168: call VecSet(x,pfive,ierr)
169: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
170: call SNESGetIterationNumber(snes,its,ierr);
171: if (rank .eq. 0) then
172: write(6,100) its
173: endif
174: 100 format('Number of Newton iterations = ',i5)
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: ! Free work space. All PETSc objects should be destroyed when they
178: ! are no longer needed.
179: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: call VecDestroy(x,ierr)
182: call VecDestroy(r,ierr)
183: call MatDestroy(J,ierr)
184: call SNESDestroy(snes,ierr)
185: call PetscFinalize(ierr)
186: end
187: !
188: ! ------------------------------------------------------------------------
189: !
190: ! FormFunction - Evaluates nonlinear function, F(x).
191: !
192: ! Input Parameters:
193: ! snes - the SNES context
194: ! x - input vector
195: ! dummy - optional user-defined context (not used here)
196: !
197: ! Output Parameter:
198: ! f - function vector
199: !
200: subroutine FormFunction(snes,x,f,dummy,ierr)
201: implicit none
203: #include include/finclude/petsc.h
204: #include include/finclude/petscvec.h
205: #include include/finclude/petscsnes.h
207: SNES snes
208: Vec x,f
209: PetscErrorCode ierr
210: integer dummy(*)
212: ! Declarations for use with local arrays
214: PetscScalar lx_v(1),lf_v(1)
215: PetscOffset lx_i,lf_i
217: ! Get pointers to vector data.
218: ! - For default PETSc vectors, VecGetArray() returns a pointer to
219: ! the data array. Otherwise, the routine is implementation dependent.
220: ! - You MUST call VecRestoreArray() when you no longer need access to
221: ! the array.
222: ! - Note that the Fortran interface to VecGetArray() differs from the
223: ! C version. See the Fortran chapter of the users manual for details.
225: call VecGetArray(x,lx_v,lx_i,ierr)
226: call VecGetArray(f,lf_v,lf_i,ierr)
228: ! Compute function
230: lf_a(1) = lx_a(1)*lx_a(1) &
231: & + lx_a(1)*lx_a(2) - 3.0
232: lf_a(2) = lx_a(1)*lx_a(2) &
233: & + lx_a(2)*lx_a(2) - 6.0
235: ! Restore vectors
237: call VecRestoreArray(x,lx_v,lx_i,ierr)
238: call VecRestoreArray(f,lf_v,lf_i,ierr)
240: return
241: end
243: ! ---------------------------------------------------------------------
244: !
245: ! FormJacobian - Evaluates Jacobian matrix.
246: !
247: ! Input Parameters:
248: ! snes - the SNES context
249: ! x - input vector
250: ! dummy - optional user-defined context (not used here)
251: !
252: ! Output Parameters:
253: ! A - Jacobian matrix
254: ! B - optionally different preconditioning matrix
255: ! flag - flag indicating matrix structure
256: !
257: subroutine FormJacobian(snes,X,jac,B,flag,dummy,ierr)
258: implicit none
260: #include include/finclude/petsc.h
261: #include include/finclude/petscvec.h
262: #include include/finclude/petscmat.h
263: #include include/finclude/petscpc.h
264: #include include/finclude/petscsnes.h
266: SNES snes
267: Vec X
268: Mat jac,B
269: MatStructure flag
270: PetscScalar A(4)
271: PetscErrorCode ierr
272: PetscInt idx(2),i2
273: integer dummy(*)
275: ! Declarations for use with local arrays
277: PetscScalar lx_v(1)
278: PetscOffset lx_i
280: ! Get pointer to vector data
282: i2 = 2
283: call VecGetArray(x,lx_v,lx_i,ierr)
285: ! Compute Jacobian entries and insert into matrix.
286: ! - Since this is such a small problem, we set all entries for
287: ! the matrix at once.
288: ! - Note that MatSetValues() uses 0-based row and column numbers
289: ! in Fortran as well as in C (as set here in the array idx).
291: idx(1) = 0
292: idx(2) = 1
293: A(1) = 2.0*lx_a(1) + lx_a(2)
294: A(2) = lx_a(1)
295: A(3) = lx_a(2)
296: A(4) = lx_a(1) + 2.0*lx_a(2)
297: call MatSetValues(jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
298: flag = SAME_NONZERO_PATTERN
300: ! Restore vector
302: call VecRestoreArray(x,lx_v,lx_i,ierr)
304: ! Assemble matrix
306: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
307: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
309: return
310: end
313: subroutine MyLineSearch(snes,lctx,x,f,g,y,w,fnorm,ynorm,gnorm, &
314: & flag,ierr)
315: #include include/finclude/petsc.h
316: #include include/finclude/petscvec.h
317: #include include/finclude/petscmat.h
318: #include include/finclude/petscksp.h
319: #include include/finclude/petscpc.h
320: #include include/finclude/petscsnes.h
322: SNES snes
323: integer lctx
324: Vec x, f,g, y, w
325: double precision fnorm,ynorm,gnorm
326: PetscTruth flag
327: PetscErrorCode ierr
329: PetscScalar mone
331: mone = -1.0d0
332: flag = 0
333: call VecNorm(y,NORM_2,ynorm,ierr)
334: call VecAYPX(y,mone,x,ierr)
335: call SNESComputeFunction(snes,y,g,ierr)
336: call VecNorm(g,NORM_2,gnorm,ierr)
337: return
338: end