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: ! dummy test of logging a reduction
55: #if defined(PETSC_USE_LOG)
56: ierr = PetscAReduce()
57: #endif
58: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
59: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
60: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif
62: i2 = 2
63: i20 = 20
64: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
65: ! Create nonlinear solver context
66: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
68: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,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_DEFAULT_REAL,PETSC_DEFAULT_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(SNESLineSearchShellSetUserFunc(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: ! View solver converged reason; we could instead use the option -snes_converged_reason
138: PetscCallA(SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr))
140: PetscCallA(SNESGetIterationNumber(snes,its,ierr))
141: if (rank .eq. 0) then
142: write(6,100) its
143: endif
144: 100 format('Number of SNES iterations = ',i5)
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Free work space. All PETSc objects should be destroyed when they
148: ! are no longer needed.
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: PetscCallA(VecDestroy(x,ierr))
152: PetscCallA(VecDestroy(r,ierr))
153: PetscCallA(MatDestroy(J,ierr))
154: PetscCallA(SNESDestroy(snes,ierr))
155: #if defined(PETSC_USE_LOG)
156: PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr))
157: PetscCallA(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr))
158: PetscCallA(PetscLogView(viewer,ierr))
159: PetscCallA(PetscViewerDestroy(viewer,ierr))
160: #endif
161: PetscCallA(PetscFinalize(ierr))
162: end
163: !
164: ! ------------------------------------------------------------------------
165: !
166: ! FormFunction - Evaluates nonlinear function, F(x).
167: !
168: ! Input Parameters:
169: ! snes - the SNES context
170: ! x - input vector
171: ! dummy - optional user-defined context (not used here)
172: !
173: ! Output Parameter:
174: ! f - function vector
175: !
176: subroutine FormFunction(snes,x,f,dummy,ierr)
177: use petscsnes
178: implicit none
180: SNES snes
181: Vec x,f
182: PetscErrorCode ierr
183: integer dummy(*)
185: ! Declarations for use with local arrays
186: PetscScalar,pointer :: lx_v(:),lf_v(:)
188: ! Get pointers to vector data.
189: ! - VecGetArrayF90() returns a pointer to the data array.
190: ! - You MUST call VecRestoreArrayF90() when you no longer need access to
191: ! the array.
193: PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
194: PetscCall(VecGetArrayF90(f,lf_v,ierr))
196: ! Compute function
198: lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
199: lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0
201: ! Restore vectors
203: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
204: PetscCall(VecRestoreArrayF90(f,lf_v,ierr))
206: return
207: end
209: ! ---------------------------------------------------------------------
210: !
211: ! FormJacobian - Evaluates Jacobian matrix.
212: !
213: ! Input Parameters:
214: ! snes - the SNES context
215: ! x - input vector
216: ! dummy - optional user-defined context (not used here)
217: !
218: ! Output Parameters:
219: ! A - Jacobian matrix
220: ! B - optionally different preconditioning matrix
221: !
222: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
223: use petscsnes
224: implicit none
226: SNES snes
227: Vec X
228: Mat jac,B
229: PetscScalar A(4)
230: PetscErrorCode ierr
231: PetscInt idx(2),i2
232: integer dummy(*)
234: ! Declarations for use with local arrays
236: PetscScalar,pointer :: lx_v(:)
238: ! Get pointer to vector data
240: i2 = 2
241: PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
243: ! Compute Jacobian entries and insert into matrix.
244: ! - Since this is such a small problem, we set all entries for
245: ! the matrix at once.
246: ! - Note that MatSetValues() uses 0-based row and column numbers
247: ! in Fortran as well as in C (as set here in the array idx).
249: idx(1) = 0
250: idx(2) = 1
251: A(1) = 2.0*lx_v(1) + lx_v(2)
252: A(2) = lx_v(1)
253: A(3) = lx_v(2)
254: A(4) = lx_v(1) + 2.0*lx_v(2)
255: PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr))
257: ! Restore vector
259: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
261: ! Assemble matrix
263: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
264: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
265: if (B .ne. jac) then
266: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
267: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
268: endif
270: return
271: end
273: subroutine MyLineSearch(linesearch, lctx, ierr)
274: use petscsnes
275: implicit none
277: SNESLineSearch linesearch
278: SNES snes
279: integer lctx
280: Vec x, f,g, y, w
281: PetscReal ynorm,gnorm,xnorm
282: PetscErrorCode ierr
284: PetscScalar mone
286: mone = -1.0
287: PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
288: PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
289: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
290: PetscCall(VecAXPY(x,mone,y,ierr))
291: PetscCall(SNESComputeFunction(snes,x,f,ierr))
292: PetscCall(VecNorm(f,NORM_2,gnorm,ierr))
293: PetscCall(VecNorm(x,NORM_2,xnorm,ierr))
294: PetscCall(VecNorm(y,NORM_2,ynorm,ierr))
295: PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr))
296: return
297: end
299: !/*TEST
300: !
301: ! test:
302: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
303: ! requires: !single
304: !
305: !TEST*/