Actual source code: ex1f.F90
1: !
2: ! Description: Uses the Newton method to solve a two-variable system.
3: !
5: program main
6: #include <petsc/finclude/petsc.h>
7: use petsc
8: implicit none
10: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
11: ! Variable declarations
12: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
13: !
14: ! Variables:
15: ! snes - nonlinear solver
16: ! ksp - linear solver
17: ! pc - preconditioner context
18: ! ksp - Krylov subspace method context
19: ! x, r - solution, residual vectors
20: ! J - Jacobian matrix
21: ! its - iterations for convergence
22: !
23: SNES snes
24: PC pc
25: KSP ksp
26: Vec x, r
27: Mat J
28: SNESLineSearch linesearch
29: PetscErrorCode ierr
30: PetscInt its, i2, i20
31: PetscMPIInt size, rank
32: PetscScalar pfive
33: PetscReal tol
34: PetscBool setls
35: PetscReal, pointer :: rhistory(:)
36: PetscInt, pointer :: itshistory(:)
37: PetscInt nhistory
38: #if defined(PETSC_USE_LOG)
39: PetscViewer viewer
40: #endif
41: double precision threshold, oldthreshold
43: ! Note: Any user-defined Fortran routines (such as FormJacobian)
44: ! MUST be declared as external.
46: external FormFunction, FormJacobian, MyLineSearch
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: ! Beginning of program
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: PetscCallA(PetscInitialize(ierr))
53: PetscCallA(PetscLogNestedBegin(ierr))
54: threshold = 1.0
55: PetscCallA(PetscLogSetThreshold(threshold, oldthreshold, ierr))
56: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
57: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
58: PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'Uniprocessor example')
60: i2 = 2
61: i20 = 20
62: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
63: ! Create nonlinear solver context
64: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
66: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
68: PetscCallA(SNESSetConvergenceHistory(snes, PETSC_NULL_REAL_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_DECIDE, PETSC_FALSE, ierr))
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: ! Create matrix and vector data structures; set corresponding routines
72: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: ! Create vectors for solution and nonlinear function
76: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i2, x, ierr))
77: PetscCallA(VecDuplicate(x, r, ierr))
79: ! Create Jacobian matrix data structure
81: PetscCallA(MatCreate(PETSC_COMM_SELF, J, ierr))
82: PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, i2, i2, ierr))
83: PetscCallA(MatSetFromOptions(J, ierr))
84: PetscCallA(MatSetUp(J, ierr))
86: ! Set function evaluation routine and vector
88: PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))
90: ! Set Jacobian matrix data structure and Jacobian evaluation routine
92: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Customize nonlinear solver; set runtime options
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: ! Set linear solver defaults for this problem. By extracting the
99: ! KSP, KSP, and PC contexts from the SNES context, we can then
100: ! directly call any KSP, KSP, and PC routines to set various options.
102: PetscCallA(SNESGetKSP(snes, ksp, ierr))
103: PetscCallA(KSPGetPC(ksp, pc, ierr))
104: PetscCallA(PCSetType(pc, PCNONE, ierr))
105: tol = 1.e-4
106: PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, i20, ierr))
108: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
109: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
110: ! These options will override those specified above as long as
111: ! SNESSetFromOptions() is called _after_ any other customization
112: ! routines.
114: PetscCallA(SNESSetFromOptions(snes, ierr))
116: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-setls', setls, ierr))
118: if (setls) then
119: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
120: PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
121: PetscCallA(SNESLineSearchShellSetApply(linesearch, MyLineSearch, 0, ierr))
122: end if
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Evaluate initial guess; then solve nonlinear system
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! Note: The user should initialize the vector, x, with the initial guess
129: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
130: ! to employ an initial guess of zero, the user should explicitly set
131: ! this vector to zero by calling VecSet().
133: pfive = 0.5
134: PetscCallA(VecSet(x, pfive, ierr))
135: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
137: PetscCallA(SNESGetConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
138: PetscCallA(SNESRestoreConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
140: ! View solver converged reason; we could instead use the option -snes_converged_reason
141: PetscCallA(SNESConvergedReasonView(snes, PETSC_VIEWER_STDOUT_WORLD, ierr))
143: PetscCallA(SNESGetIterationNumber(snes, its, ierr))
144: if (rank == 0) then
145: write (6, 100) its
146: end if
147: 100 format('Number of SNES iterations = ', i5)
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Free work space. All PETSc objects should be destroyed when they
151: ! are no longer needed.
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: PetscCallA(VecDestroy(x, ierr))
155: PetscCallA(VecDestroy(r, ierr))
156: PetscCallA(MatDestroy(J, ierr))
157: PetscCallA(SNESDestroy(snes, ierr))
158: #if defined(PETSC_USE_LOG)
159: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD, 'filename.xml', viewer, ierr))
160: PetscCallA(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_XML, ierr))
161: PetscCallA(PetscLogView(viewer, ierr))
162: PetscCallA(PetscViewerDestroy(viewer, ierr))
163: #endif
164: PetscCallA(PetscFinalize(ierr))
165: end
166: !
167: ! ------------------------------------------------------------------------
168: !
169: ! FormFunction - Evaluates nonlinear function, F(x).
170: !
171: ! Input Parameters:
172: ! snes - the SNES context
173: ! x - input vector
174: ! dummy - optional user-defined context (not used here)
175: !
176: ! Output Parameter:
177: ! f - function vector
178: !
179: subroutine FormFunction(snes, x, f, dummy, ierr)
180: use petscvec
181: use petscsnesdef
182: implicit none
184: SNES snes
185: Vec x, f
186: PetscErrorCode ierr
187: integer dummy(*)
189: ! Declarations for use with local arrays
190: PetscScalar, pointer :: lx_v(:), lf_v(:)
192: ! Get pointers to vector data.
193: ! - VecGetArray() returns a pointer to the data array.
194: ! - You MUST call VecRestoreArray() when you no longer need access to
195: ! the array.
197: PetscCall(VecGetArrayRead(x, lx_v, ierr))
198: PetscCall(VecGetArray(f, lf_v, ierr))
200: ! Compute function
202: lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
203: lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0
205: ! Restore vectors
207: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
208: PetscCall(VecRestoreArray(f, lf_v, ierr))
210: end
212: ! ---------------------------------------------------------------------
213: !
214: ! FormJacobian - Evaluates Jacobian matrix.
215: !
216: ! Input Parameters:
217: ! snes - the SNES context
218: ! x - input vector
219: ! dummy - optional user-defined context (not used here)
220: !
221: ! Output Parameters:
222: ! A - Jacobian matrix
223: ! B - optionally different matrix used to construct the preconditioner
224: !
225: subroutine FormJacobian(snes, X, jac, B, dummy, ierr)
226: use petscvec
227: use petscmat
228: use petscsnesdef
229: implicit none
231: SNES snes
232: Vec X
233: Mat jac, B
234: PetscScalar A(4)
235: PetscErrorCode ierr
236: PetscInt idx(2), i2
237: integer dummy(*)
239: ! Declarations for use with local arrays
241: PetscScalar, pointer :: lx_v(:)
243: ! Get pointer to vector data
245: i2 = 2
246: PetscCall(VecGetArrayRead(x, lx_v, ierr))
248: ! Compute Jacobian entries and insert into matrix.
249: ! - Since this is such a small problem, we set all entries for
250: ! the matrix at once.
251: ! - Note that MatSetValues() uses 0-based row and column numbers
252: ! in Fortran as well as in C (as set here in the array idx).
254: idx(1) = 0
255: idx(2) = 1
256: A(1) = 2.0*lx_v(1) + lx_v(2)
257: A(2) = lx_v(1)
258: A(3) = lx_v(2)
259: A(4) = lx_v(1) + 2.0*lx_v(2)
260: PetscCall(MatSetValues(B, i2, idx, i2, idx, A, INSERT_VALUES, ierr))
262: ! Restore vector
264: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
266: ! Assemble matrix
268: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
269: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
270: if (B /= jac) then
271: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
272: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
273: end if
275: end
277: subroutine MyLineSearch(linesearch, lctx, ierr)
278: use petscsnes
279: use petscvec
280: implicit none
282: SNESLineSearch linesearch
283: SNES snes
284: integer lctx
285: Vec x, f, g, y, w
286: PetscReal ynorm, gnorm, xnorm
287: PetscErrorCode ierr
289: PetscScalar mone
291: mone = -1.0
292: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
293: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
294: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
295: PetscCall(VecAXPY(x, mone, y, ierr))
296: PetscCall(SNESComputeFunction(snes, x, f, ierr))
297: PetscCall(VecNorm(f, NORM_2, gnorm, ierr))
298: PetscCall(VecNorm(x, NORM_2, xnorm, ierr))
299: PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
300: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, ierr))
301: end
303: !/*TEST
304: !
305: ! test:
306: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
307: ! requires: !single
308: !
309: !TEST*/