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 .eq. 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: endif
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 .eq. 0) then
145: write(6,100) its
146: endif
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 preconditioning matrix
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 .ne. jac) then
271: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
272: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
273: endif
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*/