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: ! Beginning of program
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: PetscCallA(PetscInitialize(ierr))
51: PetscCallA(PetscLogNestedBegin(ierr))
52: threshold = 1.0
53: PetscCallA(PetscLogSetThreshold(threshold,oldthreshold,ierr))
54: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
55: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
56: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example')
58: i2 = 2
59: i20 = 20
60: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Create nonlinear solver context
62: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
64: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: ! Create matrix and vector data structures; set corresponding routines
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: ! Create vectors for solution and nonlinear function
72: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr))
73: PetscCallA(VecDuplicate(x,r,ierr))
75: ! Create Jacobian matrix data structure
77: PetscCallA(MatCreate(PETSC_COMM_SELF,J,ierr))
78: PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr))
79: PetscCallA(MatSetFromOptions(J,ierr))
80: PetscCallA(MatSetUp(J,ierr))
82: ! Set function evaluation routine and vector
84: PetscCallA(SNESSetFunction(snes,r,FormFunction,0,ierr))
86: ! Set Jacobian matrix data structure and Jacobian evaluation routine
88: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Customize nonlinear solver; set runtime options
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: ! Set linear solver defaults for this problem. By extracting the
95: ! KSP, KSP, and PC contexts from the SNES context, we can then
96: ! directly call any KSP, KSP, and PC routines to set various options.
98: PetscCallA(SNESGetKSP(snes,ksp,ierr))
99: PetscCallA(KSPGetPC(ksp,pc,ierr))
100: PetscCallA(PCSetType(pc,PCNONE,ierr))
101: tol = 1.e-4
102: PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,i20,ierr))
104: ! Set SNES/KSP/KSP/PC runtime options, e.g.,
105: ! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
106: ! These options will override those specified above as long as
107: ! SNESSetFromOptions() is called _after_ any other customization
108: ! routines.
110: PetscCallA(SNESSetFromOptions(snes,ierr))
112: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-setls',setls,ierr))
114: if (setls) then
115: PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
116: PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
117: PetscCallA(SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,0,ierr))
118: endif
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Evaluate initial guess; then solve nonlinear system
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! Note: The user should initialize the vector, x, with the initial guess
125: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
126: ! to employ an initial guess of zero, the user should explicitly set
127: ! this vector to zero by calling VecSet().
129: pfive = 0.5
130: PetscCallA(VecSet(x,pfive,ierr))
131: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
133: ! View solver converged reason; we could instead use the option -snes_converged_reason
134: PetscCallA(SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr))
136: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
137: if (rank .eq. 0) then
138: write(6,100) its
139: endif
140: 100 format('Number of SNES iterations = ',i5)
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: ! Free work space. All PETSc objects should be destroyed when they
144: ! are no longer needed.
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: PetscCallA(VecDestroy(x,ierr))
148: PetscCallA(VecDestroy(r,ierr))
149: PetscCallA(MatDestroy(J,ierr))
150: PetscCallA(SNESDestroy(snes,ierr))
151: #if defined(PETSC_USE_LOG)
152: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr))
153: PetscCallA(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr))
154: PetscCallA(PetscLogView(viewer,ierr))
155: PetscCallA(PetscViewerDestroy(viewer,ierr))
156: #endif
157: PetscCallA(PetscFinalize(ierr))
158: end
159: !
160: ! ------------------------------------------------------------------------
161: !
162: ! FormFunction - Evaluates nonlinear function, F(x).
163: !
164: ! Input Parameters:
165: ! snes - the SNES context
166: ! x - input vector
167: ! dummy - optional user-defined context (not used here)
168: !
169: ! Output Parameter:
170: ! f - function vector
171: !
172: subroutine FormFunction(snes,x,f,dummy,ierr)
173: use petscsnes
174: implicit none
176: SNES snes
177: Vec x,f
178: PetscErrorCode ierr
179: integer dummy(*)
181: ! Declarations for use with local arrays
182: PetscScalar,pointer :: lx_v(:),lf_v(:)
184: ! Get pointers to vector data.
185: ! - VecGetArrayF90() returns a pointer to the data array.
186: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
187: ! the array.
189: PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
190: PetscCall(VecGetArrayF90(f,lf_v,ierr))
192: ! Compute function
194: lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
195: lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0
197: ! Restore vectors
199: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
200: PetscCall(VecRestoreArrayF90(f,lf_v,ierr))
202: return
203: end
205: ! ---------------------------------------------------------------------
206: !
207: ! FormJacobian - Evaluates Jacobian matrix.
208: !
209: ! Input Parameters:
210: ! snes - the SNES context
211: ! x - input vector
212: ! dummy - optional user-defined context (not used here)
213: !
214: ! Output Parameters:
215: ! A - Jacobian matrix
216: ! B - optionally different preconditioning matrix
217: !
218: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
219: use petscsnes
220: implicit none
222: SNES snes
223: Vec X
224: Mat jac,B
225: PetscScalar A(4)
226: PetscErrorCode ierr
227: PetscInt idx(2),i2
228: integer dummy(*)
230: ! Declarations for use with local arrays
232: PetscScalar,pointer :: lx_v(:)
234: ! Get pointer to vector data
236: i2 = 2
237: PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
239: ! Compute Jacobian entries and insert into matrix.
240: ! - Since this is such a small problem, we set all entries for
241: ! the matrix at once.
242: ! - Note that MatSetValues() uses 0-based row and column numbers
243: ! in Fortran as well as in C (as set here in the array idx).
245: idx(1) = 0
246: idx(2) = 1
247: A(1) = 2.0*lx_v(1) + lx_v(2)
248: A(2) = lx_v(1)
249: A(3) = lx_v(2)
250: A(4) = lx_v(1) + 2.0*lx_v(2)
251: PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr))
253: ! Restore vector
255: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
257: ! Assemble matrix
259: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
260: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
261: if (B .ne. jac) then
262: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
263: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
264: endif
266: return
267: end
269: subroutine MyLineSearch(linesearch, lctx, ierr)
270: use petscsnes
271: implicit none
273: SNESLineSearch linesearch
274: SNES snes
275: integer lctx
276: Vec x, f,g, y, w
277: PetscReal ynorm,gnorm,xnorm
278: PetscErrorCode ierr
280: PetscScalar mone
282: mone = -1.0
283: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
284: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
285: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
286: PetscCall(VecAXPY(x,mone,y,ierr))
287: PetscCall(SNESComputeFunction(snes,x,f,ierr))
288: PetscCall(VecNorm(f,NORM_2,gnorm,ierr))
289: PetscCall(VecNorm(x,NORM_2,xnorm,ierr))
290: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
291: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr))
292: return
293: end
295: !/*TEST
296: !
297: ! test:
298: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
299: ! requires: !single
300: !
301: !TEST*/