Actual source code: ex12f.F
1: !
2: ! "$Id: ex12.F,v 1.60 2001/08/07 03:04:13 balay Exp $";
3: !
4: ! This example demonstrates basic use of the SNES Fortran interface.
5: !
6: ! Note: The program ex10.f is the same as this example, except that it
7: ! uses the Fortran .f suffix rather than the .F suffix.
8: !
9: ! In this example the application context is a Fortran integer array:
10: ! ctx(1) = da - distributed array
11: ! 2 = F - global vector where the function is stored
12: ! 3 = xl - local work vector
13: ! 4 = rank - processor rank
14: ! 5 = size - number of processors
15: ! 6 = N - system size
16: !
17: ! Note: Any user-defined Fortran routines (such as FormJacobian)
18: ! MUST be declared as external.
19: !
20: !
21: ! Macros to make setting/getting values into vector clearer.
22: ! 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))
26: program main
27: implicit none
29: #include include/finclude/petsc.h
30: #include include/finclude/petscvec.h
31: #include include/finclude/petscda.h
32: #include include/finclude/petscmat.h
33: #include include/finclude/petscsnes.h
35: PetscFortranAddr ctx(6)
36: integer rank,size,ierr,N,start,end,nn,i,ii,its,flg
37: SNES snes
38: Mat J
39: Vec x,r,u
40: PetscScalar xp,FF,UU,h
41: character*(10) matrixname
42: external FormJacobian,FormFunction
44: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
45: N = 10
46: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',N,flg,ierr)
47: h = 1.d0/(N-1.d0)
48: ctx(6) = N
50: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
51: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
52: ctx(4) = rank
53: ctx(5) = size
55: ! Set up data structures
56: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,N,1,1, &
57: & PETSC_NULL_INTEGER,ctx(1),ierr)
59: call DACreateGlobalVector(ctx(1),x,ierr)
60: call DACreateLocalVector(ctx(1),ctx(3),ierr)
62: call PetscObjectSetName(x,'Approximate Solution',ierr)
63: call VecDuplicate(x,r,ierr)
64: call VecDuplicate(x,ctx(2),ierr)
65: call VecDuplicate(x,U,ierr)
66: call PetscObjectSetName(U,'Exact Solution',ierr)
68: call MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N, &
69: & N,3,PETSC_NULL_INTEGER,0,PETSC_NULL_INTEGER,J,ierr)
71: call MatGetType(J,matrixname,ierr)
73: ! Store right-hand-side of PDE and exact solution
74: call VecGetOwnershipRange(x,start,end,ierr)
75: xp = h*start
76: nn = end - start
77: ii = start
78: do 10, i=0,nn-1
79: FF = 6.0*xp + (xp+1.e-12)**6.e0
80: UU = xp*xp*xp
81: call VecSetValues(ctx(2),1,ii,FF,INSERT_VALUES,ierr)
82: call VecSetValues(U,1,ii,UU,INSERT_VALUES,ierr)
83: xp = xp + h
84: ii = ii + 1
85: 10 continue
86: call VecAssemblyBegin(ctx(2),ierr)
87: call VecAssemblyEnd(ctx(2),ierr)
88: call VecAssemblyBegin(U,ierr)
89: call VecAssemblyEnd(U,ierr)
91: ! Create nonlinear solver
92: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
94: ! Set various routines and options
95: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
96: call SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr)
97: call SNESSetFromOptions(snes,ierr)
99: ! Solve nonlinear system
100: call FormInitialGuess(snes,x,ierr)
101: call SNESSolve(snes,x,ierr)
102: call SNESGetIterationNumber(snes,its,ierr);
104: ! Write results if first processor
105: if (ctx(4) .eq. 0) then
106: write(6,100) its
107: endif
108: 100 format('Number of Newton iterations = ',i5)
110: ! Free work space. All PETSc objects should be destroyed when they
111: ! are no longer needed.
112: call VecDestroy(x,ierr)
113: call VecDestroy(ctx(3),ierr)
114: call VecDestroy(r,ierr)
115: call VecDestroy(U,ierr)
116: call VecDestroy(ctx(2),ierr)
117: call MatDestroy(J,ierr)
118: call SNESDestroy(snes,ierr)
119: call DADestroy(ctx(1),ierr)
120: call PetscFinalize(ierr)
121: end
124: ! -------------------- Evaluate Function F(x) ---------------------
126: subroutine FormFunction(snes,x,f,ctx,ierr)
127: implicit none
128: SNES snes
129: Vec x,f
130: PetscFortranAddr ctx(*)
131: integer rank,size,i,s,n,ierr
132: PetscOffset ixx,iff,iF2
133: PetscScalar h,d,vf2(1),vxx(1),vff(1)
134: #include include/finclude/petsc.h
135: #include include/finclude/petscvec.h
136: #include include/finclude/petscda.h
137: #include include/finclude/petscmat.h
138: #include include/finclude/petscsnes.h
141: rank = ctx(4)
142: size = ctx(5)
143: h = 1.d0/(ctx(6) - 1.d0)
144: call DAGlobalToLocalBegin(ctx(1),x,INSERT_VALUES,ctx(3),ierr)
145: call DAGlobalToLocalEnd(ctx(1),x,INSERT_VALUES,ctx(3),ierr)
147: call VecGetLocalSize(ctx(3),n,ierr)
148: if (n .gt. 1000) then
149: print*, 'Local work array not big enough'
150: call MPI_Abort(PETSC_COMM_WORLD,0,ierr)
151: endif
153: !
154: ! This sets the index ixx so that vxx(ixx+1) is the first local
155: ! element in the vector indicated by ctx(3).
156: !
157: call VecGetArray(ctx(3),vxx,ixx,ierr)
158: call VecGetArray(f,vff,iff,ierr)
159: call VecGetArray(ctx(2),vF2,iF2,ierr)
161: d = h*h
163: !
164: ! Note that the array vxx() was obtained from a ghosted local vector
165: ! ctx(3) while the array vff() was obtained from the non-ghosted parallel
166: ! vector F. This is why there is a need for shift variable s. Since vff()
167: ! does not have locations for the ghost variables we need to index in it
168: ! slightly different then indexing into vxx(). For example on processor
169: ! 1 (the second processor)
170: !
171: ! xx(1) xx(2) xx(3) .....
172: ! ^^^^^^^ ^^^^^ ^^^^^
173: ! ghost value 1st local value 2nd local value
174: !
175: ! ff(1) ff(2)
176: ! ^^^^^^^ ^^^^^^^
177: ! 1st local value 2nd local value
178: !
179: if (rank .eq. 0) then
180: s = 0
181: ff(1) = xx(1)
182: else
183: s = 1
184: endif
186: do 10 i=1,n-2
187: ff(i-s+1) = d*(xx(i) - 2.d0*xx(i+1) &
188: & + xx(i+2)) + xx(i+1)*xx(i+1) &
189: & - F2(i-s+1)
190: 10 continue
192: if (rank .eq. size-1) then
193: ff(n-s) = xx(n) - 1.d0
194: endif
196: call VecRestoreArray(f,vff,iff,ierr)
197: call VecRestoreArray(ctx(3),vxx,ixx,ierr)
198: call VecRestoreArray(ctx(2),vF2,iF2,ierr)
199: return
200: end
202: ! -------------------- Form initial approximation -----------------
204: subroutine FormInitialGuess(snes,x,ierr)
205: implicit none
206: #include include/finclude/petsc.h
207: #include include/finclude/petscvec.h
208: #include include/finclude/petscsnes.h
209: integer ierr
210: Vec x
211: SNES snes
212: PetscScalar five
214: five = 5.d-1
215: call VecSet(five,x,ierr)
216: return
217: end
219: ! -------------------- Evaluate Jacobian --------------------
221: subroutine FormJacobian(snes,x,jac,B,flag,ctx,ierr)
222: implicit none
223: #include include/finclude/petsc.h
224: #include include/finclude/petscvec.h
225: #include include/finclude/petscda.h
226: #include include/finclude/petscmat.h
227: #include include/finclude/petscsnes.h
228: SNES snes
229: Vec x
230: Mat jac,B
231: PetscFortranAddr ctx(*)
232: integer flag,ii,istart
233: PetscOffset ixx
234: integer iend,i,j,n,rank,size,end,start,ierr
235: PetscScalar d,A,h,vxx(1)
238: h = 1.d0/(ctx(6) - 1.d0)
239: d = h*h
240: rank = ctx(4)
241: size = ctx(5)
243: call VecGetArray(x,vxx,ixx,ierr)
244: call VecGetOwnershipRange(x,start,end,ierr)
245: n = end - start
247: if (rank .eq. 0) then
248: A = 1.0
249: call MatSetValues(jac,1,start,1,start,A,INSERT_VALUES,ierr)
250: istart = 1
251: else
252: istart = 0
253: endif
254: if (rank .eq. size-1) then
255: i = ctx(6)-1
256: A = 1.0
257: call MatSetValues(jac,1,i,1,i,A,INSERT_VALUES,ierr)
258: iend = n-1
259: else
260: iend = n
261: endif
262: do 10 i=istart,iend-1
263: ii = i + start
264: j = start + i - 1
265: call MatSetValues(jac,1,ii,1,j,d,INSERT_VALUES,ierr)
266: j = start + i + 1
267: call MatSetValues(jac,1,ii,1,j,d,INSERT_VALUES,ierr)
268: A = -2.0*d + 2.0*xx(i+1)
269: call MatSetValues(jac,1,ii,1,ii,A,INSERT_VALUES,ierr)
270: 10 continue
271: call VecRestoreArray(x,vxx,ixx,ierr)
272: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
273: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
274: return
275: end