Actual source code: ex1f.F90

  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*/

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

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

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

 48:       external FormFunction, FormJacobian, MyLineSearch

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

 66:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 67:       if (ierr .ne. 0) then
 68:         print*,'Unable to initialize PETSc'
 69:         stop
 70:       endif
 71:       call PetscLogNestedBegin(ierr);CHKERRA(ierr)
 72:       threshold = 1.0
 73:       call PetscLogSetThreshold(threshold,oldthreshold,ierr)
 74: ! dummy test of logging a reduction
 75: #if defined(PETSC_USE_LOG)
 76:       PetscAReduce()
 77: #endif
 78:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 79:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 80:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif

 82:       i2  = 2
 83:       i20 = 20
 84: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 85: !  Create nonlinear solver context
 86: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 88:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

 90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 91: !  Create matrix and vector data structures; set corresponding routines
 92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 94: !  Create vectors for solution and nonlinear function

 96:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
 97:       call VecDuplicate(x,r,ierr)

 99: !  Create Jacobian matrix data structure

101:       call MatCreate(PETSC_COMM_SELF,J,ierr)
102:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
103:       call MatSetFromOptions(J,ierr)
104:       call MatSetUp(J,ierr)

106: !  Set function evaluation routine and vector

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

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

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

114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: !  Customize nonlinear solver; set runtime options
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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


136:       call SNESSetFromOptions(snes,ierr)

138:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
139:      &                         '-setls',setls,ierr)

141:       if (setls) then
142:         call SNESGetLineSearch(snes, linesearch, ierr)
143:         call SNESLineSearchSetType(linesearch, 'shell', ierr)
144:         call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,   &
145:      &                                      0, ierr)
146:       endif

148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: !  Evaluate initial guess; then solve nonlinear system
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

157:       pfive = 0.5
158:       call VecSet(x,pfive,ierr)
159:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)

161: !  View solver converged reason; we could instead use the option -snes_converged_reason
162:       call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr)

164:       call SNESGetIterationNumber(snes,its,ierr);
165:       if (rank .eq. 0) then
166:          write(6,100) its
167:       endif
168:   100 format('Number of SNES iterations = ',i5)

170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: !  Free work space.  All PETSc objects should be destroyed when they
172: !  are no longer needed.
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

175:       call VecDestroy(x,ierr)
176:       call VecDestroy(r,ierr)
177:       call MatDestroy(J,ierr)
178:       call SNESDestroy(snes,ierr)
179: #if defined(PETSC_USE_LOG)
180:       call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)
181:       call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)
182:       call PetscLogView(viewer,ierr)
183:       call PetscViewerDestroy(viewer,ierr)
184: #endif
185:       call PetscFinalize(ierr)
186:       end
187: !
188: ! ------------------------------------------------------------------------
189: !
190: !  FormFunction - Evaluates nonlinear function, F(x).
191: !
192: !  Input Parameters:
193: !  snes - the SNES context
194: !  x - input vector
195: !  dummy - optional user-defined context (not used here)
196: !
197: !  Output Parameter:
198: !  f - function vector
199: !
200:       subroutine FormFunction(snes,x,f,dummy,ierr)
201:       use petscsnes
202:       implicit none

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

209: !  Declarations for use with local arrays

211:       PetscScalar  lx_v(2),lf_v(2)
212:       PetscOffset  lx_i,lf_i

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

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

225: !  Compute function

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

232: !  Restore vectors

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

237:       return
238:       end

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

258:       SNES         snes
259:       Vec          X
260:       Mat          jac,B
261:       PetscScalar  A(4)
262:       PetscErrorCode ierr
263:       PetscInt idx(2),i2
264:       integer dummy(*)

266: !  Declarations for use with local arrays

268:       PetscScalar lx_v(2)
269:       PetscOffset lx_i

271: !  Get pointer to vector data

273:       i2 = 2
274:       call VecGetArrayRead(x,lx_v,lx_i,ierr)

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

282:       idx(1) = 0
283:       idx(2) = 1
284:       A(1) = 2.0*lx_a(1) + lx_a(2)
285:       A(2) = lx_a(1)
286:       A(3) = lx_a(2)
287:       A(4) = lx_a(1) + 2.0*lx_a(2)
288:       call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)

290: !  Restore vector

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

294: !  Assemble matrix

296:       call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
297:       call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
298:       if (B .ne. jac) then
299:         call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
300:         call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
301:       endif

303:       return
304:       end


307:       subroutine MyLineSearch(linesearch, lctx, ierr)
308:       use petscsnes
309:       implicit none

311:       SNESLineSearch    linesearch
312:       SNES              snes
313:       integer           lctx
314:       Vec               x, f,g, y, w
315:       PetscReal         ynorm,gnorm,xnorm
316:       PetscBool         flag
317:       PetscErrorCode    ierr

319:       PetscScalar       mone

321:       mone = -1.0
322:       call SNESLineSearchGetSNES(linesearch, snes, ierr)
323:       call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
324:       call VecNorm(y,NORM_2,ynorm,ierr)
325:       call VecAXPY(x,mone,y,ierr)
326:       call SNESComputeFunction(snes,x,f,ierr)
327:       call VecNorm(f,NORM_2,gnorm,ierr)
328:       call VecNorm(x,NORM_2,xnorm,ierr)
329:       call VecNorm(y,NORM_2,ynorm,ierr)
330:       call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,      &
331:      & ierr)
332:       flag = PETSC_FALSE
333:       return
334:       end

336: !/*TEST
337: !
338: !   test:
339: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
340: !      requires: !single
341: !
342: !TEST*/