Actual source code: ex1f.F

  1: !
  2: ! "$Id: ex1f.F,v 1.33 2001/08/07 03:04:16 balay Exp $";
  3: !
  4: !  Description: Uses the Newton method to solve a two-variable system.
  5: !
  6: !/*T
  7: !  Concepts: SNES^basic uniprocessor example
  8: !  Processors: 1
  9: !T*/
 10: !
 11: ! -----------------------------------------------------------------------

 13:       program main
 14:       implicit none

 16: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 17: !                    Include files
 18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19: !
 20: !  The following include statements are generally used in SNES Fortran
 21: !  programs:
 22: !     petsc.h       - base PETSc routines
 23: !     petscvec.h    - vectors
 24: !     petscmat.h    - matrices
 25: !     petscksp.h    - Krylov subspace methods
 26: !     petscpc.h     - preconditioners
 27: !     petscsnes.h   - SNES interface
 28: !  Other include statements may be needed if using additional PETSc
 29: !  routines in a Fortran program, e.g.,
 30: !     petscviewer.h - viewers
 31: !     petscis.h     - index sets
 32: !
 33:  #include include/finclude/petsc.h
 34:  #include include/finclude/petscvec.h
 35:  #include include/finclude/petscmat.h
 36:  #include include/finclude/petscksp.h
 37:  #include include/finclude/petscpc.h
 38:  #include include/finclude/petscsnes.h
 39: !
 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !                   Variable declarations
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !
 44: !  Variables:
 45: !     snes        - nonlinear solver
 46: !     ksp        - linear solver
 47: !     pc          - preconditioner context
 48: !     ksp         - Krylov subspace method context
 49: !     x, r        - solution, residual vectors
 50: !     J           - Jacobian matrix
 51: !     its         - iterations for convergence
 52: !
 53:       SNES     snes
 54:       PC       pc
 55:       KSP      ksp
 56:       Vec      x,r
 57:       Mat      J
 58:       integer  ierr,its,size,rank
 59:       PetscScalar   pfive
 60:       double precision   tol
 61:       PetscTruth  setls

 63: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 64: !  MUST be declared as external.

 66:       external FormFunction, FormJacobian, MyLineSearch

 68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69: !                   Macro definitions
 70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71: !
 72: !  Macros to make clearer the process of setting values in vectors and
 73: !  getting values from vectors.  These vectors are used in the routines
 74: !  FormFunction() and FormJacobian().
 75: !   - The element lx_a(ib) is element ib in the vector x
 76: !
 77: #define lx_a(ib) lx_v(lx_i + (ib))
 78: #define lf_a(ib) lf_v(lf_i + (ib))
 79: !
 80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81: !                 Beginning of program
 82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 84:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 85:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 86:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 87:       if (size .ne. 1) then
 88:          if (rank .eq. 0) then
 89:             write(6,*) 'This is a uniprocessor example only!'
 90:          endif
 91:          SETERRQ(1,' ',ierr)
 92:       endif

 94: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 95: !  Create nonlinear solver context
 96: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 98:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: !  Create matrix and vector data structures; set corresponding routines
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

104: !  Create vectors for solution and nonlinear function

106:       call VecCreateSeq(PETSC_COMM_SELF,2,x,ierr)
107:       call VecDuplicate(x,r,ierr)

109: !  Create Jacobian matrix data structure

111:       call MatCreate(PETSC_COMM_SELF,PETSC_DECIDE,PETSC_DECIDE,2,2,J,   &
112:      &               ierr)
113:       call MatSetFromOptions(J,ierr)

115: !  Set function evaluation routine and vector

117:       call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr)

119: !  Set Jacobian matrix data structure and Jacobian evaluation routine

121:       call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT,     &
122:      &     ierr)

124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: !  Customize nonlinear solver; set runtime options
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

128: !  Set linear solver defaults for this problem. By extracting the
129: !  KSP, KSP, and PC contexts from the SNES context, we can then
130: !  directly call any KSP, KSP, and PC routines to set various options.

132:       call SNESGetKSP(snes,ksp,ierr)
133:       call KSPGetPC(ksp,pc,ierr)
134:       call PCSetType(pc,PCNONE,ierr)
135:       tol = 1.e-4
136:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,     &
137:      &                      PETSC_DEFAULT_DOUBLE_PRECISION,20,ierr)

139: !  Set SNES/KSP/KSP/PC runtime options, e.g.,
140: !      -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
141: !  These options will override those specified above as long as
142: !  SNESSetFromOptions() is called _after_ any other customization
143: !  routines.


146:       call SNESSetFromOptions(snes,ierr)

148:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-setls',setls,ierr)

150:       if (setls .eq. PETSC_TRUE) then
151:         call SNESSetLineSearch(snes,MyLineSearch,                       &
152:      &                         PETSC_NULL_OBJECT,ierr)
153:       endif

155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: !  Evaluate initial guess; then solve nonlinear system
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

159: !  Note: The user should initialize the vector, x, with the initial guess
160: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
161: !  to employ an initial guess of zero, the user should explicitly set
162: !  this vector to zero by calling VecSet().

164:       pfive = 0.5
165:       call VecSet(pfive,x,ierr)
166:       call SNESSolve(snes,x,ierr)
167:       call SNESGetIterationNumber(snes,its,ierr);
168:       if (rank .eq. 0) then
169:          write(6,100) its
170:       endif
171:   100 format('Number of Newton iterations = ',i5)

173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: !  Free work space.  All PETSc objects should be destroyed when they
175: !  are no longer needed.
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

178:       call VecDestroy(x,ierr)
179:       call VecDestroy(r,ierr)
180:       call MatDestroy(J,ierr)
181:       call SNESDestroy(snes,ierr)
182:       call PetscFinalize(ierr)
183:       end
184: !
185: ! ------------------------------------------------------------------------
186: !
187: !  FormFunction - Evaluates nonlinear function, F(x).
188: !
189: !  Input Parameters:
190: !  snes - the SNES context
191: !  x - input vector
192: !  dummy - optional user-defined context (not used here)
193: !
194: !  Output Parameter:
195: !  f - function vector
196: !
197:       subroutine FormFunction(snes,x,f,dummy,ierr)
198:       implicit none

200:  #include include/finclude/petsc.h
201:  #include include/finclude/petscvec.h
202:  #include include/finclude/petscsnes.h

204:       SNES     snes
205:       Vec      x,f
206:       integer  ierr,dummy(*)

208: !  Declarations for use with local arrays

210:       PetscScalar  lx_v(1),lf_v(1)
211:       PetscOffset  lx_i,lf_i

213: !  Get pointers to vector data.
214: !    - For default PETSc vectors, VecGetArray() returns a pointer to
215: !      the data array.  Otherwise, the routine is implementation dependent.
216: !    - You MUST call VecRestoreArray() when you no longer need access to
217: !      the array.
218: !    - Note that the Fortran interface to VecGetArray() differs from the
219: !      C version.  See the Fortran chapter of the users manual for details.

221:       call VecGetArray(x,lx_v,lx_i,ierr)
222:       call VecGetArray(f,lf_v,lf_i,ierr)

224: !  Compute function

226:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
227:      &          + lx_a(1)*lx_a(2) - 3.0
228:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
229:      &          + lx_a(2)*lx_a(2) - 6.0

231: !  Restore vectors

233:       call VecRestoreArray(x,lx_v,lx_i,ierr)
234:       call VecRestoreArray(f,lf_v,lf_i,ierr)

236:       return
237:       end

239: ! ---------------------------------------------------------------------
240: !
241: !  FormJacobian - Evaluates Jacobian matrix.
242: !
243: !  Input Parameters:
244: !  snes - the SNES context
245: !  x - input vector
246: !  dummy - optional user-defined context (not used here)
247: !
248: !  Output Parameters:
249: !  A - Jacobian matrix
250: !  B - optionally different preconditioning matrix
251: !  flag - flag indicating matrix structure
252: !
253:       subroutine FormJacobian(snes,X,jac,B,flag,dummy,ierr)
254:       implicit none

256:  #include include/finclude/petsc.h
257:  #include include/finclude/petscvec.h
258:  #include include/finclude/petscmat.h
259:  #include include/finclude/petscpc.h
260:  #include include/finclude/petscsnes.h

262:       SNES         snes
263:       Vec          X
264:       Mat          jac,B
265:       MatStructure flag
266:       PetscScalar  A(4)
267:       integer      ierr,idx(2),dummy(*)

269: !  Declarations for use with local arrays

271:       PetscScalar lx_v(1)
272:       PetscOffset lx_i

274: !  Get pointer to vector data

276:       call VecGetArray(x,lx_v,lx_i,ierr)

278: !  Compute Jacobian entries and insert into matrix.
279: !   - Since this is such a small problem, we set all entries for
280: !     the matrix at once.
281: !   - Note that MatSetValues() uses 0-based row and column numbers
282: !     in Fortran as well as in C (as set here in the array idx).

284:       idx(1) = 0
285:       idx(2) = 1
286:       A(1) = 2.0*lx_a(1) + lx_a(2)
287:       A(2) = lx_a(1)
288:       A(3) = lx_a(2)
289:       A(4) = lx_a(1) + 2.0*lx_a(2)
290:       call MatSetValues(jac,2,idx,2,idx,A,INSERT_VALUES,ierr)
291:       flag = SAME_NONZERO_PATTERN

293: !  Restore vector

295:       call VecRestoreArray(x,lx_v,lx_i,ierr)

297: !  Assemble matrix

299:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
300:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)

302:       return
303:       end


306:       subroutine MyLineSearch(snes,lctx,x,f,g,y,w,fnorm,ynorm,gnorm,     &
307:      &                        flag,ierr)
308:  #include include/finclude/petsc.h
309:  #include include/finclude/petscvec.h
310:  #include include/finclude/petscmat.h
311:  #include include/finclude/petscksp.h
312:  #include include/finclude/petscpc.h
313:  #include include/finclude/petscsnes.h

315:       SNES              snes
316:       integer           lctx
317:       Vec               x, f,g, y, w
318:       double  precision fnorm,ynorm,gnorm
319:       integer           flag,ierr

321:       PetscScalar       mone

323:       mone = -1.0d0
324:       flag = 0
325:       call VecNorm(y,NORM_2,ynorm,ierr)
326:       call VecAYPX(mone,x,y,ierr)
327:       call SNESComputeFunction(snes,y,g,ierr)
328:       call VecNorm(g,NORM_2,gnorm,ierr)
329:       return
330:       end