Actual source code: ex1f.F

petsc-master 2016-09-27
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:       implicit none

 15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 16: !                    Include files
 17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 18: !
 19: !  The following include statements are generally used in SNES Fortran
 20: !  programs:
 21: !     petscsys.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 <petsc/finclude/petscsys.h>
 33:  #include <petsc/finclude/petscvec.h>
 34:  #include <petsc/finclude/petscmat.h>
 35:  #include <petsc/finclude/petscksp.h>
 36:  #include <petsc/finclude/petscpc.h>
 37:  #include <petsc/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:       SNESLineSearch linesearch
 58:       PetscErrorCode  ierr
 59:       PetscInt its,i2,i20
 60:       PetscMPIInt size,rank
 61:       PetscScalar   pfive
 62:       PetscReal   tol
 63:       PetscBool   setls

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

 68:       external FormFunction, FormJacobian, MyLineSearch

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

 86:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 87:       if (ierr .ne. 0) then
 88:         print*,'Unable to initialize PETSc'
 89:         stop
 90:       endif
 91:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 92:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 93:       if (size .ne. 1) then
 94:          if (rank .eq. 0) then
 95:             write(6,*) 'This is a uniprocessor example only!'
 96:          endif
 97:          SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
 98:       endif

100:       i2  = 2
101:       i20 = 20
102: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
103: !  Create nonlinear solver context
104: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

106:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: !  Create matrix and vector data structures; set corresponding routines
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

112: !  Create vectors for solution and nonlinear function

114:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
115:       call VecDuplicate(x,r,ierr)

117: !  Create Jacobian matrix data structure

119:       call MatCreate(PETSC_COMM_SELF,J,ierr)
120:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
121:       call MatSetFromOptions(J,ierr)
122:       call MatSetUp(J,ierr)

124: !  Set function evaluation routine and vector

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

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

130:       call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT,     &
131:      &     ierr)

133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: !  Customize nonlinear solver; set runtime options
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

141:       call SNESGetKSP(snes,ksp,ierr)
142:       call KSPGetPC(ksp,pc,ierr)
143:       call PCSetType(pc,PCNONE,ierr)
144:       tol = 1.e-4
145:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                  &
146:      &                      PETSC_DEFAULT_REAL,i20,ierr)

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


155:       call SNESSetFromOptions(snes,ierr)

157:       call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,   &
158:      &                         '-setls',setls,ierr)

160:       if (setls) then
161:         call SNESGetLineSearch(snes, linesearch, ierr)
162:         call SNESLineSearchSetType(linesearch, 'shell', ierr)
163:         call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,   &
164:      &PETSC_NULL_OBJECT, ierr)
165:       endif

167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: !  Evaluate initial guess; then solve nonlinear system
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

176:       pfive = 0.5
177:       call VecSet(x,pfive,ierr)
178:       call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
179:       call SNESGetIterationNumber(snes,its,ierr);
180:       if (rank .eq. 0) then
181:          write(6,100) its
182:       endif
183:   100 format('Number of SNES iterations = ',i5)

185: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: !  Free work space.  All PETSc objects should be destroyed when they
187: !  are no longer needed.
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

190:       call VecDestroy(x,ierr)
191:       call VecDestroy(r,ierr)
192:       call MatDestroy(J,ierr)
193:       call SNESDestroy(snes,ierr)
194:       call PetscFinalize(ierr)
195:       end
196: !
197: ! ------------------------------------------------------------------------
198: !
199: !  FormFunction - Evaluates nonlinear function, F(x).
200: !
201: !  Input Parameters:
202: !  snes - the SNES context
203: !  x - input vector
204: !  dummy - optional user-defined context (not used here)
205: !
206: !  Output Parameter:
207: !  f - function vector
208: !
209:       subroutine FormFunction(snes,x,f,dummy,ierr)
210:       implicit none

212:  #include <petsc/finclude/petscsys.h>
213:  #include <petsc/finclude/petscvec.h>
214:  #include <petsc/finclude/petscsnes.h>

216:       SNES     snes
217:       Vec      x,f
218:       PetscErrorCode ierr
219:       integer dummy(*)

221: !  Declarations for use with local arrays

223:       PetscScalar  lx_v(2),lf_v(2)
224:       PetscOffset  lx_i,lf_i

226: !  Get pointers to vector data.
227: !    - For default PETSc vectors, VecGetArray() returns a pointer to
228: !      the data array.  Otherwise, the routine is implementation dependent.
229: !    - You MUST call VecRestoreArray() when you no longer need access to
230: !      the array.
231: !    - Note that the Fortran interface to VecGetArray() differs from the
232: !      C version.  See the Fortran chapter of the users manual for details.

234:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
235:       call VecGetArray(f,lf_v,lf_i,ierr)

237: !  Compute function

239:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
240:      &          + lx_a(1)*lx_a(2) - 3.0
241:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
242:      &          + lx_a(2)*lx_a(2) - 6.0

244: !  Restore vectors

246:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
247:       call VecRestoreArray(f,lf_v,lf_i,ierr)

249:       return
250:       end

252: ! ---------------------------------------------------------------------
253: !
254: !  FormJacobian - Evaluates Jacobian matrix.
255: !
256: !  Input Parameters:
257: !  snes - the SNES context
258: !  x - input vector
259: !  dummy - optional user-defined context (not used here)
260: !
261: !  Output Parameters:
262: !  A - Jacobian matrix
263: !  B - optionally different preconditioning matrix
264: !  flag - flag indicating matrix structure
265: !
266:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
267:       implicit none

269:  #include <petsc/finclude/petscsys.h>
270:  #include <petsc/finclude/petscvec.h>
271:  #include <petsc/finclude/petscmat.h>
272:  #include <petsc/finclude/petscpc.h>
273:  #include <petsc/finclude/petscsnes.h>

275:       SNES         snes
276:       Vec          X
277:       Mat          jac,B
278:       PetscScalar  A(4)
279:       PetscErrorCode ierr
280:       PetscInt idx(2),i2
281:       integer dummy(*)

283: !  Declarations for use with local arrays

285:       PetscScalar lx_v(2)
286:       PetscOffset lx_i

288: !  Get pointer to vector data

290:       i2 = 2
291:       call VecGetArrayRead(x,lx_v,lx_i,ierr)

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

299:       idx(1) = 0
300:       idx(2) = 1
301:       A(1) = 2.0*lx_a(1) + lx_a(2)
302:       A(2) = lx_a(1)
303:       A(3) = lx_a(2)
304:       A(4) = lx_a(1) + 2.0*lx_a(2)
305:       call MatSetValues(jac,i2,idx,i2,idx,A,INSERT_VALUES,ierr)

307: !  Restore vector

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

311: !  Assemble matrix

313:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
314:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)

316:       return
317:       end


320:       subroutine MyLineSearch(linesearch, lctx, ierr)
321:  #include <petsc/finclude/petscsys.h>
322:  #include <petsc/finclude/petscvec.h>
323:  #include <petsc/finclude/petscmat.h>
324:  #include <petsc/finclude/petscksp.h>
325:  #include <petsc/finclude/petscpc.h>
326:  #include <petsc/finclude/petscsnes.h>

328:       SNES              snes
329:       integer           lctx
330:       Vec               x, f,g, y, w
331:       PetscReal         ynorm,gnorm,xnorm
332:       PetscBool         flag
333:       PetscErrorCode    ierr

335:       PetscScalar       mone

337:       mone = -1.0
338:       call SNESLineSearchGetSNES(linesearch, snes, ierr)
339:       call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
340:       call VecNorm(y,NORM_2,ynorm,ierr)
341:       call VecAXPY(x,mone,y,ierr)
342:       call SNESComputeFunction(snes,x,f,ierr)
343:       call VecNorm(f,NORM_2,gnorm,ierr)
344:       call VecNorm(x,NORM_2,xnorm,ierr)
345:       call VecNorm(y,NORM_2,ynorm,ierr)
346:       call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,      &
347:      & ierr)
348:       flag = PETSC_FALSE
349:       return
350:       end