Actual source code: ex1f.F

petsc-master 2017-03-23
Report Typos and Errors
  1: !
  2: !
  3: !  Description: Uses the Newton method to solve a two-variable system.
  4: !
  5: !/*T
  6: !  Concepts: SNES^basic uniprocessor example
  7: !  Processors: 1
  8: !T*/
  9: !
 10: ! -----------------------------------------------------------------------

 12:       program main
 13:  #include <petsc/finclude/petsc.h>
 14:       use petsc
 15:       implicit none


 18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19: !                   Variable declarations
 20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 21: !
 22: !  Variables:
 23: !     snes        - nonlinear solver
 24: !     ksp        - linear solver
 25: !     pc          - preconditioner context
 26: !     ksp         - Krylov subspace method context
 27: !     x, r        - solution, residual vectors
 28: !     J           - Jacobian matrix
 29: !     its         - iterations for convergence
 30: !
 31:       SNES     snes
 32:       PC       pc
 33:       KSP      ksp
 34:       Vec      x,r
 35:       Mat      J
 36:       SNESLineSearch linesearch
 37:       PetscErrorCode  ierr
 38:       PetscInt its,i2,i20
 39:       PetscMPIInt size,rank
 40:       PetscScalar   pfive
 41:       PetscReal   tol
 42:       PetscBool   setls

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

 47:       external FormFunction, FormJacobian, MyLineSearch

 49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50: !                   Macro definitions
 51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52: !
 53: !  Macros to make clearer the process of setting values in vectors and
 54: !  getting values from vectors.  These vectors are used in the routines
 55: !  FormFunction() and FormJacobian().
 56: !   - The element lx_a(ib) is element ib in the vector x
 57: !
 58: #define lx_a(ib) lx_v(lx_i + (ib))
 59: #define lf_a(ib) lf_v(lf_i + (ib))
 60: !
 61: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62: !                 Beginning of program
 63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 65:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 66:       if (ierr .ne. 0) then
 67:         print*,'Unable to initialize PETSc'
 68:         stop
 69:       endif
 70:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 71:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 72:       if (size .ne. 1) then
 73:          if (rank .eq. 0) then
 74:             write(6,*) 'This is a uniprocessor example only!'
 75:          endif
 76:          SETERRQ(PETSC_COMM_SELF,1,'')
 77:       endif

 79:       i2  = 2
 80:       i20 = 20
 81: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 82: !  Create nonlinear solver context
 83: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 85:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

 87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88: !  Create matrix and vector data structures; set corresponding routines
 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 91: !  Create vectors for solution and nonlinear function

 93:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
 94:       call VecDuplicate(x,r,ierr)

 96: !  Create Jacobian matrix data structure

 98:       call MatCreate(PETSC_COMM_SELF,J,ierr)
 99:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
100:       call MatSetFromOptions(J,ierr)
101:       call MatSetUp(J,ierr)

103: !  Set function evaluation routine and vector

105:       call SNESSetFunction(snes,r,FormFunction,0,ierr)

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

109:       call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)

111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: !  Customize nonlinear solver; set runtime options
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

119:       call SNESGetKSP(snes,ksp,ierr)
120:       call KSPGetPC(ksp,pc,ierr)
121:       call PCSetType(pc,PCNONE,ierr)
122:       tol = 1.e-4
123:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                  &
124:      &                      PETSC_DEFAULT_REAL,i20,ierr)

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


133:       call SNESSetFromOptions(snes,ierr)

135:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
136:      &                         '-setls',setls,ierr)

138:       if (setls) then
139:         call SNESGetLineSearch(snes, linesearch, ierr)
140:         call SNESLineSearchSetType(linesearch, 'shell', ierr)
141:         call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,   &
142:      &                                      0, ierr)
143:       endif

145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: !  Evaluate initial guess; then solve nonlinear system
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

154:       pfive = 0.5
155:       call VecSet(x,pfive,ierr)
156:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
157:       call SNESGetIterationNumber(snes,its,ierr);
158:       if (rank .eq. 0) then
159:          write(6,100) its
160:       endif
161:   100 format('Number of SNES iterations = ',i5)

163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: !  Free work space.  All PETSc objects should be destroyed when they
165: !  are no longer needed.
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

168:       call VecDestroy(x,ierr)
169:       call VecDestroy(r,ierr)
170:       call MatDestroy(J,ierr)
171:       call SNESDestroy(snes,ierr)
172:       call PetscFinalize(ierr)
173:       end
174: !
175: ! ------------------------------------------------------------------------
176: !
177: !  FormFunction - Evaluates nonlinear function, F(x).
178: !
179: !  Input Parameters:
180: !  snes - the SNES context
181: !  x - input vector
182: !  dummy - optional user-defined context (not used here)
183: !
184: !  Output Parameter:
185: !  f - function vector
186: !
187:       subroutine FormFunction(snes,x,f,dummy,ierr)
188:       use petscsnes
189:       implicit none

191:       SNES     snes
192:       Vec      x,f
193:       PetscErrorCode ierr
194:       integer dummy(*)

196: !  Declarations for use with local arrays

198:       PetscScalar  lx_v(2),lf_v(2)
199:       PetscOffset  lx_i,lf_i

201: !  Get pointers to vector data.
202: !    - For default PETSc vectors, VecGetArray() returns a pointer to
203: !      the data array.  Otherwise, the routine is implementation dependent.
204: !    - You MUST call VecRestoreArray() when you no longer need access to
205: !      the array.
206: !    - Note that the Fortran interface to VecGetArray() differs from the
207: !      C version.  See the Fortran chapter of the users manual for details.

209:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
210:       call VecGetArray(f,lf_v,lf_i,ierr)

212: !  Compute function

214:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
215:      &          + lx_a(1)*lx_a(2) - 3.0
216:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
217:      &          + lx_a(2)*lx_a(2) - 6.0

219: !  Restore vectors

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

224:       return
225:       end

227: ! ---------------------------------------------------------------------
228: !
229: !  FormJacobian - Evaluates Jacobian matrix.
230: !
231: !  Input Parameters:
232: !  snes - the SNES context
233: !  x - input vector
234: !  dummy - optional user-defined context (not used here)
235: !
236: !  Output Parameters:
237: !  A - Jacobian matrix
238: !  B - optionally different preconditioning matrix
239: !  flag - flag indicating matrix structure
240: !
241:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
242:       use petscsnes
243:       implicit none

245:       SNES         snes
246:       Vec          X
247:       Mat          jac,B
248:       PetscScalar  A(4)
249:       PetscErrorCode ierr
250:       PetscInt idx(2),i2
251:       integer dummy(*)

253: !  Declarations for use with local arrays

255:       PetscScalar lx_v(2)
256:       PetscOffset lx_i

258: !  Get pointer to vector data

260:       i2 = 2
261:       call VecGetArrayRead(x,lx_v,lx_i,ierr)

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

269:       idx(1) = 0
270:       idx(2) = 1
271:       A(1) = 2.0*lx_a(1) + lx_a(2)
272:       A(2) = lx_a(1)
273:       A(3) = lx_a(2)
274:       A(4) = lx_a(1) + 2.0*lx_a(2)
275:       call MatSetValues(jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)

277: !  Restore vector

279:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

281: !  Assemble matrix

283:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
284:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)

286:       return
287:       end


290:       subroutine MyLineSearch(linesearch, lctx, ierr)
291:       use petscsnes
292:       implicit none

294:       SNESLineSearch    linesearch
295:       SNES              snes
296:       integer           lctx
297:       Vec               x, f,g, y, w
298:       PetscReal         ynorm,gnorm,xnorm
299:       PetscBool         flag
300:       PetscErrorCode    ierr

302:       PetscScalar       mone

304:       mone = -1.0
305:       call SNESLineSearchGetSNES(linesearch, snes, ierr)
306:       call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
307:       call VecNorm(y,NORM_2,ynorm,ierr)
308:       call VecAXPY(x,mone,y,ierr)
309:       call SNESComputeFunction(snes,x,f,ierr)
310:       call VecNorm(f,NORM_2,gnorm,ierr)
311:       call VecNorm(x,NORM_2,xnorm,ierr)
312:       call VecNorm(y,NORM_2,ynorm,ierr)
313:       call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,      &
314:      & ierr)
315:       flag = PETSC_FALSE
316:       return
317:       end