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: PetscCallA(PetscInitialize(ierr))
63: PetscCallA(PetscLogNestedBegin(ierr))
64: threshold = 1.0
65: PetscCallA(PetscLogSetThreshold(threshold,oldthreshold,ierr))
66: ! dummy test of logging a reduction
67: #if defined(PETSC_USE_LOG)
68: PetscAReduce()
69: #endif
70: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
71: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
72: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif
74: i2 = 2
75: i20 = 20
76: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! Create nonlinear solver context
78: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
80: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! Create matrix and vector data structures; set corresponding routines
84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: ! Create vectors for solution and nonlinear function
88: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr))
89: PetscCallA(VecDuplicate(x,r,ierr))
91: ! Create Jacobian matrix data structure
93: PetscCallA(MatCreate(PETSC_COMM_SELF,J,ierr))
94: PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr))
95: PetscCallA(MatSetFromOptions(J,ierr))
96: PetscCallA(MatSetUp(J,ierr))
98: ! Set function evaluation routine and vector
100: PetscCallA(SNESSetFunction(snes,r,FormFunction,0,ierr))
102: ! Set Jacobian matrix data structure and Jacobian evaluation routine
104: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Customize nonlinear solver; set runtime options
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: ! Set linear solver defaults for this problem. By extracting the
111: ! KSP, KSP, and PC contexts from the SNES context, we can then
112: ! directly call any KSP, KSP, and PC routines to set various options.
114: PetscCallA(SNESGetKSP(snes,ksp,ierr))
115: PetscCallA(KSPGetPC(ksp,pc,ierr))
116: PetscCallA(PCSetType(pc,PCNONE,ierr))
117: tol = 1.e-4
118: PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,i20,ierr))
120: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
121: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
122: ! These options will override those specified above as long as
123: ! SNESSetFromOptions() is called _after_ any other customization
124: ! routines.
126: PetscCallA(SNESSetFromOptions(snes,ierr))
128: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-setls',setls,ierr))
130: if (setls) then
131: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
132: PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
133: PetscCallA(SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,0,ierr))
134: endif
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Evaluate initial guess; then solve nonlinear system
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: ! Note: The user should initialize the vector, x, with the initial guess
141: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
142: ! to employ an initial guess of zero, the user should explicitly set
143: ! this vector to zero by calling VecSet().
145: pfive = 0.5
146: PetscCallA(VecSet(x,pfive,ierr))
147: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
149: ! View solver converged reason; we could instead use the option -snes_converged_reason
150: PetscCallA(SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr))
152: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
153: if (rank .eq. 0) then
154: write(6,100) its
155: endif
156: 100 format('Number of SNES iterations = ',i5)
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Free work space. All PETSc objects should be destroyed when they
160: ! are no longer needed.
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: PetscCallA(VecDestroy(x,ierr))
164: PetscCallA(VecDestroy(r,ierr))
165: PetscCallA(MatDestroy(J,ierr))
166: PetscCallA(SNESDestroy(snes,ierr))
167: #if defined(PETSC_USE_LOG)
168: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr))
169: PetscCallA(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr))
170: PetscCallA(PetscLogView(viewer,ierr))
171: PetscCallA(PetscViewerDestroy(viewer,ierr))
172: #endif
173: PetscCallA(PetscFinalize(ierr))
174: end
175: !
176: ! ------------------------------------------------------------------------
177: !
178: ! FormFunction - Evaluates nonlinear function, F(x).
179: !
180: ! Input Parameters:
181: ! snes - the SNES context
182: ! x - input vector
183: ! dummy - optional user-defined context (not used here)
184: !
185: ! Output Parameter:
186: ! f - function vector
187: !
188: subroutine FormFunction(snes,x,f,dummy,ierr)
189: use petscsnes
190: implicit none
192: SNES snes
193: Vec x,f
194: PetscErrorCode ierr
195: integer dummy(*)
197: ! Declarations for use with local arrays
199: PetscScalar lx_v(2),lf_v(2)
200: PetscOffset lx_i,lf_i
202: ! Get pointers to vector data.
203: ! - For default PETSc vectors, VecGetArray() returns a pointer to
204: ! the data array. Otherwise, the routine is implementation dependent.
205: ! - You MUST call VecRestoreArray() when you no longer need access to
206: ! the array.
207: ! - Note that the Fortran interface to VecGetArray() differs from the
208: ! C version. See the Fortran chapter of the users manual for details.
210: PetscCall(VecGetArrayRead(x,lx_v,lx_i,ierr))
211: PetscCall(VecGetArray(f,lf_v,lf_i,ierr))
213: ! Compute function
215: lf_a(1) = lx_a(1)*lx_a(1) + lx_a(1)*lx_a(2) - 3.0
216: lf_a(2) = lx_a(1)*lx_a(2) + lx_a(2)*lx_a(2) - 6.0
218: ! Restore vectors
220: PetscCall(VecRestoreArrayRead(x,lx_v,lx_i,ierr))
221: PetscCall(VecRestoreArray(f,lf_v,lf_i,ierr))
223: return
224: end
226: ! ---------------------------------------------------------------------
227: !
228: ! FormJacobian - Evaluates Jacobian matrix.
229: !
230: ! Input Parameters:
231: ! snes - the SNES context
232: ! x - input vector
233: ! dummy - optional user-defined context (not used here)
234: !
235: ! Output Parameters:
236: ! A - Jacobian matrix
237: ! B - optionally different preconditioning matrix
238: !
239: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
240: use petscsnes
241: implicit none
243: SNES snes
244: Vec X
245: Mat jac,B
246: PetscScalar A(4)
247: PetscErrorCode ierr
248: PetscInt idx(2),i2
249: integer dummy(*)
251: ! Declarations for use with local arrays
253: PetscScalar lx_v(2)
254: PetscOffset lx_i
256: ! Get pointer to vector data
258: i2 = 2
259: PetscCall(VecGetArrayRead(x,lx_v,lx_i,ierr))
261: ! Compute Jacobian entries and insert into matrix.
262: ! - Since this is such a small problem, we set all entries for
263: ! the matrix at once.
264: ! - Note that MatSetValues() uses 0-based row and column numbers
265: ! in Fortran as well as in C (as set here in the array idx).
267: idx(1) = 0
268: idx(2) = 1
269: A(1) = 2.0*lx_a(1) + lx_a(2)
270: A(2) = lx_a(1)
271: A(3) = lx_a(2)
272: A(4) = lx_a(1) + 2.0*lx_a(2)
273: PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr))
275: ! Restore vector
277: PetscCall(VecRestoreArrayRead(x,lx_v,lx_i,ierr))
279: ! Assemble matrix
281: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
282: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
283: if (B .ne. jac) then
284: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
285: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
286: endif
288: return
289: end
291: subroutine MyLineSearch(linesearch, lctx, ierr)
292: use petscsnes
293: implicit none
295: SNESLineSearch linesearch
296: SNES snes
297: integer lctx
298: Vec x, f,g, y, w
299: PetscReal ynorm,gnorm,xnorm
300: PetscErrorCode ierr
302: PetscScalar mone
304: mone = -1.0
305: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
306: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
307: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
308: PetscCall(VecAXPY(x,mone,y,ierr))
309: PetscCall(SNESComputeFunction(snes,x,f,ierr))
310: PetscCall(VecNorm(f,NORM_2,gnorm,ierr))
311: PetscCall(VecNorm(x,NORM_2,xnorm,ierr))
312: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
313: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr))
314: return
315: end
317: !/*TEST
318: !
319: ! test:
320: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
321: ! requires: !single
322: !
323: !TEST*/