Actual source code: ex1f.F

  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:       implicit none

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

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

 67:       external FormFunction, FormJacobian, MyLineSearch

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

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

 95:       i2  = 2
 96:       i20 = 20
 97: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 98: !  Create nonlinear solver context
 99: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

101:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: !  Create matrix and vector data structures; set corresponding routines
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

107: !  Create vectors for solution and nonlinear function

109:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
110:       call VecDuplicate(x,r,ierr)

112: !  Create Jacobian matrix data structure

114:       call MatCreate(PETSC_COMM_SELF,J,ierr)
115:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
116:       call MatSetFromOptions(J,ierr)

118: !  Set function evaluation routine and vector

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

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

124:       call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT,     &
125:      &     ierr)

127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: !  Customize nonlinear solver; set runtime options
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

135:       call SNESGetKSP(snes,ksp,ierr)
136:       call KSPGetPC(ksp,pc,ierr)
137:       call PCSetType(pc,PCNONE,ierr)
138:       tol = 1.e-4
139:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,     &
140:      &                      PETSC_DEFAULT_DOUBLE_PRECISION,i20,ierr)

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


149:       call SNESSetFromOptions(snes,ierr)

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

153:       if (setls .eq. PETSC_TRUE) then
154:         call SNESLineSearchSet(snes,MyLineSearch,                       &
155:      &                         PETSC_NULL_OBJECT,ierr)
156:       endif

158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: !  Evaluate initial guess; then solve nonlinear system
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: !  Free work space.  All PETSc objects should be destroyed when they
178: !  are no longer needed.
179: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

181:       call VecDestroy(x,ierr)
182:       call VecDestroy(r,ierr)
183:       call MatDestroy(J,ierr)
184:       call SNESDestroy(snes,ierr)
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:       implicit none

203:  #include include/finclude/petsc.h
204:  #include include/finclude/petscvec.h
205:  #include include/finclude/petscsnes.h

207:       SNES     snes
208:       Vec      x,f
209:       PetscErrorCode ierr
210:       integer dummy(*)

212: !  Declarations for use with local arrays

214:       PetscScalar  lx_v(1),lf_v(1)
215:       PetscOffset  lx_i,lf_i

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

225:       call VecGetArray(x,lx_v,lx_i,ierr)
226:       call VecGetArray(f,lf_v,lf_i,ierr)

228: !  Compute function

230:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
231:      &          + lx_a(1)*lx_a(2) - 3.0
232:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
233:      &          + lx_a(2)*lx_a(2) - 6.0

235: !  Restore vectors

237:       call VecRestoreArray(x,lx_v,lx_i,ierr)
238:       call VecRestoreArray(f,lf_v,lf_i,ierr)

240:       return
241:       end

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

260:  #include include/finclude/petsc.h
261:  #include include/finclude/petscvec.h
262:  #include include/finclude/petscmat.h
263:  #include include/finclude/petscpc.h
264:  #include include/finclude/petscsnes.h

266:       SNES         snes
267:       Vec          X
268:       Mat          jac,B
269:       MatStructure flag
270:       PetscScalar  A(4)
271:       PetscErrorCode ierr
272:       PetscInt idx(2),i2
273:       integer dummy(*)

275: !  Declarations for use with local arrays

277:       PetscScalar lx_v(1)
278:       PetscOffset lx_i

280: !  Get pointer to vector data

282:       i2 = 2
283:       call VecGetArray(x,lx_v,lx_i,ierr)

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

291:       idx(1) = 0
292:       idx(2) = 1
293:       A(1) = 2.0*lx_a(1) + lx_a(2)
294:       A(2) = lx_a(1)
295:       A(3) = lx_a(2)
296:       A(4) = lx_a(1) + 2.0*lx_a(2)
297:       call MatSetValues(jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)
298:       flag = SAME_NONZERO_PATTERN

300: !  Restore vector

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

304: !  Assemble matrix

306:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
307:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)

309:       return
310:       end


313:       subroutine MyLineSearch(snes,lctx,x,f,g,y,w,fnorm,ynorm,gnorm,     &
314:      &                        flag,ierr)
315:  #include include/finclude/petsc.h
316:  #include include/finclude/petscvec.h
317:  #include include/finclude/petscmat.h
318:  #include include/finclude/petscksp.h
319:  #include include/finclude/petscpc.h
320:  #include include/finclude/petscsnes.h

322:       SNES              snes
323:       integer           lctx
324:       Vec               x, f,g, y, w
325:       double  precision fnorm,ynorm,gnorm
326:       PetscTruth           flag
327:       PetscErrorCode ierr

329:       PetscScalar       mone

331:       mone = -1.0d0
332:       flag = 0
333:       call VecNorm(y,NORM_2,ynorm,ierr)
334:       call VecAYPX(y,mone,x,ierr)
335:       call SNESComputeFunction(snes,y,g,ierr)
336:       call VecNorm(g,NORM_2,gnorm,ierr)
337:       return
338:       end