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: ! petscsnes.h - SNES interface
28: ! Other include statements may be needed if using additional PETSc
29: ! routines in a Fortran program, e.g.,
30: ! petscviewer.h - viewers
31: ! petscis.h - index sets
32: !
33: #include include/finclude/petsc.h
34: #include include/finclude/petscvec.h
35: #include include/finclude/petscmat.h
36: #include include/finclude/petscksp.h
37: #include include/finclude/petscpc.h
38: #include include/finclude/petscsnes.h
39: !
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: ! Variable declarations
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: !
44: ! Variables:
45: ! snes - nonlinear solver
46: ! ksp - linear solver
47: ! pc - preconditioner context
48: ! ksp - Krylov subspace method context
49: ! x, r - solution, residual vectors
50: ! J - Jacobian matrix
51: ! its - iterations for convergence
52: !
53: SNES snes
54: PC pc
55: KSP ksp
56: Vec x,r
57: Mat J
58: integer ierr,its,size,rank
59: PetscScalar pfive
60: double precision tol
61: PetscTruth setls
63: ! Note: Any user-defined Fortran routines (such as FormJacobian)
64: ! MUST be declared as external.
66: external FormFunction, FormJacobian, MyLineSearch
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: ! Macro definitions
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: !
72: ! Macros to make clearer the process of setting values in vectors and
73: ! getting values from vectors. These vectors are used in the routines
74: ! FormFunction() and FormJacobian().
75: ! - The element lx_a(ib) is element ib in the vector x
76: !
77: #define lx_a(ib) lx_v(lx_i + (ib))
78: #define lf_a(ib) lf_v(lf_i + (ib))
79: !
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: ! Beginning of program
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
85: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
86: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
87: if (size .ne. 1) then
88: if (rank .eq. 0) then
89: write(6,*) 'This is a uniprocessor example only!'
90: endif
91: SETERRQ(1,' ',ierr)
92: endif
94: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Create nonlinear solver context
96: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
98: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Create matrix and vector data structures; set corresponding routines
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! Create vectors for solution and nonlinear function
106: call VecCreateSeq(PETSC_COMM_SELF,2,x,ierr)
107: call VecDuplicate(x,r,ierr)
109: ! Create Jacobian matrix data structure
111: call MatCreate(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,2,2,J, &
112: & ierr)
113: call MatSetFromOptions(J,ierr)
115: ! Set function evaluation routine and vector
117: call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr)
119: ! Set Jacobian matrix data structure and Jacobian evaluation routine
121: call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT, &
122: & ierr)
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Customize nonlinear solver; set runtime options
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! Set linear solver defaults for this problem. By extracting the
129: ! KSP, KSP, and PC contexts from the SNES context, we can then
130: ! directly call any KSP, KSP, and PC routines to set various options.
132: call SNESGetKSP(snes,ksp,ierr)
133: call KSPGetPC(ksp,pc,ierr)
134: call PCSetType(pc,PCNONE,ierr)
135: tol = 1.e-4
136: call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, &
137: & PETSC_DEFAULT_DOUBLE_PRECISION,20,ierr)
139: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
140: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
141: ! These options will override those specified above as long as
142: ! SNESSetFromOptions() is called _after_ any other customization
143: ! routines.
146: call SNESSetFromOptions(snes,ierr)
148: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-setls',setls,ierr)
150: if (setls .eq. PETSC_TRUE) then
151: call SNESSetLineSearch(snes,MyLineSearch, &
152: & PETSC_NULL_OBJECT,ierr)
153: endif
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: ! Evaluate initial guess; then solve nonlinear system
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Note: The user should initialize the vector, x, with the initial guess
160: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
161: ! to employ an initial guess of zero, the user should explicitly set
162: ! this vector to zero by calling VecSet().
164: pfive = 0.5
165: call VecSet(pfive,x,ierr)
166: call SNESSolve(snes,x,ierr)
167: call SNESGetIterationNumber(snes,its,ierr);
168: if (rank .eq. 0) then
169: write(6,100) its
170: endif
171: 100 format('Number of Newton iterations = ',i5)
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: ! Free work space. All PETSc objects should be destroyed when they
175: ! are no longer needed.
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: call VecDestroy(x,ierr)
179: call VecDestroy(r,ierr)
180: call MatDestroy(J,ierr)
181: call SNESDestroy(snes,ierr)
182: call PetscFinalize(ierr)
183: end
184: !
185: ! ------------------------------------------------------------------------
186: !
187: ! FormFunction - Evaluates nonlinear function, F(x).
188: !
189: ! Input Parameters:
190: ! snes - the SNES context
191: ! x - input vector
192: ! dummy - optional user-defined context (not used here)
193: !
194: ! Output Parameter:
195: ! f - function vector
196: !
197: subroutine FormFunction(snes,x,f,dummy,ierr)
198: implicit none
200: #include include/finclude/petsc.h
201: #include include/finclude/petscvec.h
202: #include include/finclude/petscsnes.h
204: SNES snes
205: Vec x,f
206: integer ierr,dummy(*)
208: ! Declarations for use with local arrays
210: PetscScalar lx_v(1),lf_v(1)
211: PetscOffset lx_i,lf_i
213: ! Get pointers to vector data.
214: ! - For default PETSc vectors, VecGetArray() returns a pointer to
215: ! the data array. Otherwise, the routine is implementation dependent.
216: ! - You MUST call VecRestoreArray() when you no longer need access to
217: ! the array.
218: ! - Note that the Fortran interface to VecGetArray() differs from the
219: ! C version. See the Fortran chapter of the users manual for details.
221: call VecGetArray(x,lx_v,lx_i,ierr)
222: call VecGetArray(f,lf_v,lf_i,ierr)
224: ! Compute function
226: lf_a(1) = lx_a(1)*lx_a(1) &
227: & + lx_a(1)*lx_a(2) - 3.0
228: lf_a(2) = lx_a(1)*lx_a(2) &
229: & + lx_a(2)*lx_a(2) - 6.0
231: ! Restore vectors
233: call VecRestoreArray(x,lx_v,lx_i,ierr)
234: call VecRestoreArray(f,lf_v,lf_i,ierr)
236: return
237: end
239: ! ---------------------------------------------------------------------
240: !
241: ! FormJacobian - Evaluates Jacobian matrix.
242: !
243: ! Input Parameters:
244: ! snes - the SNES context
245: ! x - input vector
246: ! dummy - optional user-defined context (not used here)
247: !
248: ! Output Parameters:
249: ! A - Jacobian matrix
250: ! B - optionally different preconditioning matrix
251: ! flag - flag indicating matrix structure
252: !
253: subroutine FormJacobian(snes,X,jac,B,flag,dummy,ierr)
254: implicit none
256: #include include/finclude/petsc.h
257: #include include/finclude/petscvec.h
258: #include include/finclude/petscmat.h
259: #include include/finclude/petscpc.h
260: #include include/finclude/petscsnes.h
262: SNES snes
263: Vec X
264: Mat jac,B
265: MatStructure flag
266: PetscScalar A(4)
267: integer ierr,idx(2),dummy(*)
269: ! Declarations for use with local arrays
271: PetscScalar lx_v(1)
272: PetscOffset lx_i
274: ! Get pointer to vector data
276: call VecGetArray(x,lx_v,lx_i,ierr)
278: ! Compute Jacobian entries and insert into matrix.
279: ! - Since this is such a small problem, we set all entries for
280: ! the matrix at once.
281: ! - Note that MatSetValues() uses 0-based row and column numbers
282: ! in Fortran as well as in C (as set here in the array idx).
284: idx(1) = 0
285: idx(2) = 1
286: A(1) = 2.0*lx_a(1) + lx_a(2)
287: A(2) = lx_a(1)
288: A(3) = lx_a(2)
289: A(4) = lx_a(1) + 2.0*lx_a(2)
290: call MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES,ierr)
291: flag = SAME_NONZERO_PATTERN
293: ! Restore vector
295: call VecRestoreArray(x,lx_v,lx_i,ierr)
297: ! Assemble matrix
299: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
300: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
302: return
303: end
306: subroutine MyLineSearch(snes,lctx,x,f,g,y,w,fnorm,ynorm,gnorm, &
307: & flag,ierr)
308: #include include/finclude/petsc.h
309: #include include/finclude/petscvec.h
310: #include include/finclude/petscmat.h
311: #include include/finclude/petscksp.h
312: #include include/finclude/petscpc.h
313: #include include/finclude/petscsnes.h
315: SNES snes
316: integer lctx
317: Vec x, f,g, y, w
318: double precision fnorm,ynorm,gnorm
319: integer flag,ierr
321: PetscScalar mone
323: mone = -1.0d0
324: flag = 0
325: call VecNorm(y,NORM_2,ynorm,ierr)
326: call VecAYPX(mone,x,y,ierr)
327: call SNESComputeFunction(snes,y,g,ierr)
328: call VecNorm(g,NORM_2,gnorm,ierr)
329: return
330: end