Actual source code: ex14f.F
1: !
2: !
3: ! This example demonstrates use of the SNES Fortran interface.
4: !
5: ! Note: The program is modified from ex12f.F
6: ! Used for testing MUMPS interface, and control when the Jacobian is rebuilt
7: !
8: ! In this example the application context is a Fortran integer array:
9: ! ctx(1) = da - distributed array
10: ! 2 = F - global vector where the function is stored
11: ! 3 = xl - local work vector
12: ! 4 = rank - processor rank
13: ! 5 = size - number of processors
14: ! 6 = N - system size
15: !
16: ! Note: Any user-defined Fortran routines (such as FormJacobian)
17: ! MUST be declared as external.
18: !
19: !
20: ! Macros to make setting/getting values into vector clearer.
21: ! The element xx(ib) is the ibth element in the vector indicated by ctx(3)
23: #define xx(ib) vxx(ixx + (ib))
24: #define ff(ib) vff(iff + (ib))
25: #define F2(ib) vF2(iF2 + (ib))
27: module Petscmod
28: implicit none
29: #include <finclude/petscsys.h>
30: #include <finclude/petscvec.h>
31: #include <finclude/petscvec.h90>
32: #include <finclude/petscmat.h>
33: #include <finclude/petscmat.h90>
34: #include <finclude/petscviewer.h>
35: #include <finclude/petscksp.h>
36: #include <finclude/petscpc.h>
37: #include <finclude/petscsnes.h>
38: #include <finclude/petscis.h>
39: #include <finclude/petscdmda.h>
40: end module Petscmod
42: module Snesmonitormod
43: use Petscmod
44: implicit none
45: save
46: type monctx
47: PetscInt :: its,lag
48: end type monctx
49: end module Snesmonitormod
51: ! ---------------------------------------------------------------------
52: ! ---------------------------------------------------------------------
53: ! Subroutine FormMonitor
54: ! This function lets up keep track of the SNES progress at each step
55: ! In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag'
56: !
57: ! Input Parameters:
58: ! snes - SNES nonlinear solver context
59: ! its - current nonlinear iteration, starting from a call of SNESSolve()
60: ! norm - 2-norm of current residual (may be approximate)
61: ! snesmonitor - monctx designed module (included in Snesmonitormod)
62: ! ---------------------------------------------------------------------
63: subroutine FormMonitor(snes,its,norm,snesmonitor,ierr)
65: use Snesmonitormod
66: implicit none
68: SNES :: snes
69: PetscInt :: its
70: PetscScalar :: norm
71: type(monctx) :: snesmonitor
72: PetscErrorCode :: ierr
74: c write(6,*) " "
75: c write(6,*) " its ",its,snesmonitor%its,"lag",
76: c & snesmonitor%lag
77: c call flush(6)
78: if (mod(snesmonitor%its,snesmonitor%lag).eq.0)
79: & then
80: call SNESSetLagJacobian(snes,1,ierr) ! build jacobian
81: else
82: call SNESSetLagJacobian(snes,-1,ierr) ! do NOT build jacobian
83: endif
84: snesmonitor%its = snesmonitor%its + 1
85: end subroutine FormMonitor
87: ! ---------------------------------------------------------------------
88: ! ---------------------------------------------------------------------
89: program main
91: use Petscmod
92: use Snesmonitormod
94: implicit none
96: type(monctx) :: snesmonitor
97: PetscFortranAddr ctx(6)
98: PetscMPIInt rank,size
99: PetscErrorCode ierr
100: PetscInt N,start,end,nn,i
101: PetscInt ii,its,i1,i0,i3
102: PetscBool flg
103: SNES snes
104: PC pc
105: KSP ksp
106: Mat J,F
107: Vec x,r,u
108: PetscScalar xp,FF,UU,h
109: character*(10) matrixname
110: external FormJacobian,FormFunction
111: external FormMonitor
113: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
114: i1 = 1
115: i0 = 0
116: i3 = 3
117: N = 10
118: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',N,flg,ierr)
119: h = 1.d0/(N-1.d0)
120: ctx(6) = N
122: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
123: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
124: ctx(4) = rank
125: ctx(5) = size
127: ! Set up data structures
128: call DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,N,i1,i1, &
129: & PETSC_NULL_INTEGER,ctx(1),ierr)
131: call DMCreateGlobalVector(ctx(1),x,ierr)
132: call DMCreateLocalVector(ctx(1),ctx(3),ierr)
134: call PetscObjectSetName(x,'Approximate Solution',ierr)
135: call VecDuplicate(x,r,ierr)
136: call VecDuplicate(x,ctx(2),ierr)
137: call VecDuplicate(x,U,ierr)
138: call PetscObjectSetName(U,'Exact Solution',ierr)
140: call MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N, &
141: & N,i3,PETSC_NULL_INTEGER,i0,PETSC_NULL_INTEGER,J,ierr)
143: call MatGetType(J,matrixname,ierr)
145: ! Store right-hand-side of PDE and exact solution
146: call VecGetOwnershipRange(x,start,end,ierr)
147: xp = h*start
148: nn = end - start
149: ii = start
150: do 10, i=0,nn-1
151: FF = 6.0*xp + (xp+1.e-12)**6.e0
152: UU = xp*xp*xp
153: call VecSetValues(ctx(2),i1,ii,FF,INSERT_VALUES,ierr)
154: call VecSetValues(U,i1,ii,UU,INSERT_VALUES,ierr)
155: xp = xp + h
156: ii = ii + 1
157: 10 continue
158: call VecAssemblyBegin(ctx(2),ierr)
159: call VecAssemblyEnd(ctx(2),ierr)
160: call VecAssemblyBegin(U,ierr)
161: call VecAssemblyEnd(U,ierr)
163: ! Create nonlinear solver
164: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
166: ! Set various routines and options
167: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
168: call SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr)
169: call SNESSetLagJacobian(snes,3,ierr)
171: ! Set linear solver defaults for this problem.
172: call SNESGetKSP(snes,ksp,ierr)
173: call KSPGetPC(ksp,pc,ierr)
174: #ifdef PETSC_HAVE_MUMPS
175: call PCSetType(pc,PCLU,ierr)
176: call PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS,ierr)
177: #endif
179: snesmonitor%its = 0
180: call SNESGetLagJacobian(snes,snesmonitor%lag,ierr)
181: call SNESMonitorSet(snes,FormMonitor,snesmonitor,
182: & PETSC_NULL_FUNCTION,ierr)
183: call SNESSetFromOptions(snes,ierr)
184: call FormInitialGuess(snes,x,ierr)
186: ! Setup nonlinear solver, and call SNESSolve() for first few iterations, which calls SNESSetUp() :-(
187: !--------------------------------------------------------------------------------------------------
188: call SNESSetTolerances(snes,PETSC_DEFAULT_DOUBLE_PRECISION,
189: & PETSC_DEFAULT_DOUBLE_PRECISION,
190: & PETSC_DEFAULT_DOUBLE_PRECISION,
191: & 3,PETSC_DEFAULT_INTEGER,ierr)
192: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
194: #ifdef PETSC_HAVE_MUMPS
195: ! Get PCFactor to set MUMPS matrix ordering option
196: call PCFactorGetMatrix(pc,F,ierr)
197: call MatMumpsSetIcntl(F,7,2,ierr) ! must be called before next SNESSetUp? -- not being used!!!
198: #endif
200: ! Solve nonlinear system
201: call SNESSetTolerances(snes,PETSC_DEFAULT_DOUBLE_PRECISION,
202: & PETSC_DEFAULT_DOUBLE_PRECISION,
203: & PETSC_DEFAULT_DOUBLE_PRECISION,
204: & 1000,PETSC_DEFAULT_INTEGER,ierr)
206: ! Call SNESSolve() for next few iterations
207: !--------------------------------------------------
208: snesmonitor%its = snesmonitor%its - 1 !do not count the 0-th iteration
209: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
210: call SNESGetIterationNumber(snes,its,ierr);
212: ! Write results if first processor
213: if (ctx(4) .eq. 0) then
214: write(6,100) its
215: endif
216: 100 format('Number of Newton iterations = ',i5)
218: ! Free work space. All PETSc objects should be destroyed when they
219: ! are no longer needed.
220: call VecDestroy(x,ierr)
221: call VecDestroy(ctx(3),ierr)
222: call VecDestroy(r,ierr)
223: call VecDestroy(U,ierr)
224: call VecDestroy(ctx(2),ierr)
225: call MatDestroy(J,ierr)
226: call SNESDestroy(snes,ierr)
227: call DMDestroy(ctx(1),ierr)
228: call PetscFinalize(ierr)
229: end
232: ! -------------------- Evaluate Function F(x) ---------------------
234: subroutine FormFunction(snes,x,f,ctx,ierr)
235: implicit none
236: SNES snes
237: Vec x,f
238: PetscFortranAddr ctx(*)
239: PetscMPIInt rank,size
240: PetscInt i,s,n
241: PetscErrorCode ierr
242: PetscOffset ixx,iff,iF2
243: PetscScalar h,d,vf2(1),vxx(1),vff(1)
244: #include <finclude/petscsys.h>
245: #include <finclude/petscvec.h>
246: #include <finclude/petscdmda.h>
247: #include <finclude/petscmat.h>
248: #include <finclude/petscsnes.h>
251: rank = ctx(4)
252: size = ctx(5)
253: h = 1.d0/(ctx(6) - 1.d0)
254: call DMGlobalToLocalBegin(ctx(1),x,INSERT_VALUES,ctx(3),ierr)
255: call DMGlobalToLocalEnd(ctx(1),x,INSERT_VALUES,ctx(3),ierr)
257: call VecGetLocalSize(ctx(3),n,ierr)
258: if (n .gt. 1000) then
259: print*, 'Local work array not big enough'
260: call MPI_Abort(PETSC_COMM_WORLD,0,ierr)
261: endif
263: !
264: ! This sets the index ixx so that vxx(ixx+1) is the first local
265: ! element in the vector indicated by ctx(3).
266: !
267: call VecGetArray(ctx(3),vxx,ixx,ierr)
268: call VecGetArray(f,vff,iff,ierr)
269: call VecGetArray(ctx(2),vF2,iF2,ierr)
271: d = h*h
273: !
274: ! Note that the array vxx() was obtained from a ghosted local vector
275: ! ctx(3) while the array vff() was obtained from the non-ghosted parallel
276: ! vector F. This is why there is a need for shift variable s. Since vff()
277: ! does not have locations for the ghost variables we need to index in it
278: ! slightly different then indexing into vxx(). For example on processor
279: ! 1 (the second processor)
280: !
281: ! xx(1) xx(2) xx(3) .....
282: ! ^^^^^^^ ^^^^^ ^^^^^
283: ! ghost value 1st local value 2nd local value
284: !
285: ! ff(1) ff(2)
286: ! ^^^^^^^ ^^^^^^^
287: ! 1st local value 2nd local value
288: !
289: if (rank .eq. 0) then
290: s = 0
291: ff(1) = xx(1)
292: else
293: s = 1
294: endif
296: do 10 i=1,n-2
297: ff(i-s+1) = d*(xx(i) - 2.d0*xx(i+1) &
298: & + xx(i+2)) + xx(i+1)*xx(i+1) &
299: & - F2(i-s+1)
300: 10 continue
302: if (rank .eq. size-1) then
303: ff(n-s) = xx(n) - 1.d0
304: endif
306: call VecRestoreArray(f,vff,iff,ierr)
307: call VecRestoreArray(ctx(3),vxx,ixx,ierr)
308: call VecRestoreArray(ctx(2),vF2,iF2,ierr)
309: return
310: end
312: ! -------------------- Form initial approximation -----------------
314: subroutine FormInitialGuess(snes,x,ierr)
315: implicit none
316: #include <finclude/petscsys.h>
317: #include <finclude/petscvec.h>
318: #include <finclude/petscsnes.h>
319: PetscErrorCode ierr
320: Vec x
321: SNES snes
322: PetscScalar five
324: five = 5.d-1
325: call VecSet(x,five,ierr)
326: return
327: end
329: ! -------------------- Evaluate Jacobian --------------------
331: subroutine FormJacobian(snes,x,jac,B,flag,ctx,ierr)
333: use Petscmod
334: implicit none
336: SNES snes
337: Vec x
338: Mat jac,B
339: PetscFortranAddr ctx(*)
340: PetscBool flag
341: PetscInt ii,istart,iend
342: PetscInt i,j,n,end,start,i1
343: PetscErrorCode ierr
344: PetscMPIInt rank,size
345: PetscOffset ixx
346: PetscScalar d,A,h,vxx(1)
348: rank = ctx(4)
349: size = ctx(5)
350: if (rank .eq. 0) then
351: write(6,*)" Jacobian is built ..."
352: call flush(6)
353: endif
354: i1 = 1
355: h = 1.d0/(ctx(6) - 1.d0)
356: d = h*h
357: rank = ctx(4)
358: size = ctx(5)
360: call VecGetArray(x,vxx,ixx,ierr)
361: call VecGetOwnershipRange(x,start,end,ierr)
362: n = end - start
364: if (rank .eq. 0) then
365: A = 1.0
366: call MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr)
367: istart = 1
368: else
369: istart = 0
370: endif
371: if (rank .eq. size-1) then
372: i = ctx(6)-1
373: A = 1.0
374: call MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr)
375: iend = n-1
376: else
377: iend = n
378: endif
379: do 10 i=istart,iend-1
380: ii = i + start
381: j = start + i - 1
382: call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
383: j = start + i + 1
384: call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
385: A = -2.0*d + 2.0*xx(i+1)
386: call MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr)
387: 10 continue
388: call VecRestoreArray(x,vxx,ixx,ierr)
389: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
390: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
391: return
392: end