Actual source code: ex1f.F
1: !
2: ! "$Id: ex1f.F,v 1.33 2001/08/07 03:04:16 balay Exp $";
3: !
4: ! Description: Uses the Newton method to solve a two-variable system.
5: !
6: !/*T
7: ! Concepts: SNES^basic uniprocessor example
8: ! Processors: 1
9: !T*/
10: !
11: ! -----------------------------------------------------------------------
13: program main
14: implicit none
16: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
17: ! Include files
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: !
20: ! The following include statements are generally used in SNES Fortran
21: ! programs:
22: ! petsc.h - base PETSc routines
23: ! petscvec.h - vectors
24: ! petscmat.h - matrices
25: ! petscksp.h - Krylov subspace methods
26: ! petscpc.h - preconditioners
27: ! petscsles.h - SLES interface
28: ! petscsnes.h - SNES interface
29: ! Other include statements may be needed if using additional PETSc
30: ! routines in a Fortran program, e.g.,
31: ! petscviewer.h - viewers
32: ! petscis.h - index sets
33: !
34: #include include/finclude/petsc.h
35: #include include/finclude/petscvec.h
36: #include include/finclude/petscmat.h
37: #include include/finclude/petscksp.h
38: #include include/finclude/petscpc.h
39: #include include/finclude/petscsles.h
40: #include include/finclude/petscsnes.h
41: !
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Variable declarations
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! Variables:
47: ! snes - nonlinear solver
48: ! sles - linear solver
49: ! pc - preconditioner context
50: ! ksp - Krylov subspace method context
51: ! x, r - solution, residual vectors
52: ! J - Jacobian matrix
53: ! its - iterations for convergence
54: !
55: SNES snes
56: SLES sles
57: PC pc
58: KSP ksp
59: Vec x,r
60: Mat J
61: integer ierr,its,size,rank
62: PetscScalar pfive
63: double precision tol
64: PetscTruth setls
66: ! Note: Any user-defined Fortran routines (such as FormJacobian)
67: ! MUST be declared as external.
69: external FormFunction, FormJacobian, MyLineSearch
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: ! Macro definitions
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: !
75: ! Macros to make clearer the process of setting values in vectors and
76: ! getting values from vectors. These vectors are used in the routines
77: ! FormFunction() and FormJacobian().
78: ! - The element lx_a(ib) is element ib in the vector x
79: !
80: #define lx_a(ib) lx_v(lx_i + (ib))
81: #define lf_a(ib) lf_v(lf_i + (ib))
82: !
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Beginning of program
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
88: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
89: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
90: if (size .ne. 1) then
91: if (rank .eq. 0) then
92: write(6,*) 'This is a uniprocessor example only!'
93: endif
94: SETERRQ(1,' ',ierr)
95: endif
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,2,x,ierr)
110: call VecDuplicate(x,r,ierr)
112: ! Create Jacobian matrix data structure
114: call MatCreate(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,2,2,J, &
115: & 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: ! SLES, KSP, and PC contexts from the SNES context, we can then
133: ! directly call any SLES, KSP, and PC routines to set various options.
135: call SNESGetSLES(snes,sles,ierr)
136: call SLESGetKSP(sles,ksp,ierr)
137: call SLESGetPC(sles,pc,ierr)
138: call PCSetType(pc,PCNONE,ierr)
139: tol = 1.e-4
140: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
141: & PETSC_DEFAULT_DOUBLE_PRECISION,20,ierr)
143: ! Set SNES/SLES/KSP/PC runtime options, e.g.,
144: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
145: ! These options will override those specified above as long as
146: ! SNESSetFromOptions() is called _after_ any other customization
147: ! routines.
150: call SNESSetFromOptions(snes,ierr)
152: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-setls',setls,ierr)
154: if (setls .eq. PETSC_TRUE) then
155: call SNESSetLineSearch(snes,MyLineSearch, &
156: & PETSC_NULL_OBJECT,ierr)
157: endif
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Evaluate initial guess; then solve nonlinear system
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: ! Note: The user should initialize the vector, x, with the initial guess
164: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
165: ! to employ an initial guess of zero, the user should explicitly set
166: ! this vector to zero by calling VecSet().
168: pfive = 0.5
169: call VecSet(pfive,x,ierr)
170: call SNESSolve(snes,x,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: ! FormFunction - Evaluates nonlinear function, F(x).
190: !
191: ! Input Parameters:
192: ! snes - the SNES context
193: ! x - input vector
194: ! dummy - optional user-defined context (not used here)
195: !
196: ! Output Parameter:
197: ! f - function vector
198: !
199: subroutine FormFunction(snes,x,f,dummy,ierr)
200: implicit none
202: #include include/finclude/petsc.h
203: #include include/finclude/petscvec.h
204: #include include/finclude/petscsnes.h
206: SNES snes
207: Vec x,f
208: integer ierr,dummy(*)
210: ! Declarations for use with local arrays
212: PetscScalar lx_v(1),lf_v(1)
213: PetscOffset lx_i,lf_i
215: ! Get pointers to vector data.
216: ! - For default PETSc vectors, VecGetArray() returns a pointer to
217: ! the data array. Otherwise, the routine is implementation dependent.
218: ! - You MUST call VecRestoreArray() when you no longer need access to
219: ! the array.
220: ! - Note that the Fortran interface to VecGetArray() differs from the
221: ! C version. See the Fortran chapter of the users manual for details.
223: call VecGetArray(x,lx_v,lx_i,ierr)
224: call VecGetArray(f,lf_v,lf_i,ierr)
226: ! Compute function
228: lf_a(1) = lx_a(1)*lx_a(1) &
229: & + lx_a(1)*lx_a(2) - 3.0
230: lf_a(2) = lx_a(1)*lx_a(2) &
231: & + lx_a(2)*lx_a(2) - 6.0
233: ! Restore vectors
235: call VecRestoreArray(x,lx_v,lx_i,ierr)
236: call VecRestoreArray(f,lf_v,lf_i,ierr)
238: return
239: end
241: ! ---------------------------------------------------------------------
242: !
243: ! FormJacobian - Evaluates Jacobian matrix.
244: !
245: ! Input Parameters:
246: ! snes - the SNES context
247: ! x - input vector
248: ! dummy - optional user-defined context (not used here)
249: !
250: ! Output Parameters:
251: ! A - Jacobian matrix
252: ! B - optionally different preconditioning matrix
253: ! flag - flag indicating matrix structure
254: !
255: subroutine FormJacobian(snes,X,jac,B,flag,dummy,ierr)
256: implicit none
258: #include include/finclude/petsc.h
259: #include include/finclude/petscvec.h
260: #include include/finclude/petscmat.h
261: #include include/finclude/petscpc.h
262: #include include/finclude/petscsnes.h
264: SNES snes
265: Vec X
266: Mat jac,B
267: MatStructure flag
268: PetscScalar A(4)
269: integer ierr,idx(2),dummy(*)
271: ! Declarations for use with local arrays
273: PetscScalar lx_v(1)
274: PetscOffset lx_i
276: ! Get pointer to vector data
278: call VecGetArray(x,lx_v,lx_i,ierr)
280: ! Compute Jacobian entries and insert into matrix.
281: ! - Since this is such a small problem, we set all entries for
282: ! the matrix at once.
283: ! - Note that MatSetValues() uses 0-based row and column numbers
284: ! in Fortran as well as in C (as set here in the array idx).
286: idx(1) = 0
287: idx(2) = 1
288: A(1) = 2.0*lx_a(1) + lx_a(2)
289: A(2) = lx_a(1)
290: A(3) = lx_a(2)
291: A(4) = lx_a(1) + 2.0*lx_a(2)
292: call MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES,ierr)
293: flag = SAME_NONZERO_PATTERN
295: ! Restore vector
297: call VecRestoreArray(x,lx_v,lx_i,ierr)
299: ! Assemble matrix
301: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
302: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
304: return
305: end
308: subroutine MyLineSearch(snes,lctx,x,f,g,y,w,fnorm,ynorm,gnorm, &
309: & flag,ierr)
310: #include include/finclude/petsc.h
311: #include include/finclude/petscvec.h
312: #include include/finclude/petscmat.h
313: #include include/finclude/petscksp.h
314: #include include/finclude/petscpc.h
315: #include include/finclude/petscsles.h
316: #include include/finclude/petscsnes.h
318: SNES snes
319: integer lctx
320: Vec x, f,g, y, w
321: double precision fnorm,ynorm,gnorm
322: integer flag,ierr
324: PetscScalar mone
326: mone = -1.0d0
327: flag = 0
328: call VecNorm(y,NORM_2,ynorm,ierr)
329: call VecAYPX(mone,x,y,ierr)
330: call SNESComputeFunction(snes,y,g,ierr)
331: call VecNorm(g,NORM_2,gnorm,ierr)
332: return
333: end