Actual source code: ex1f.F

petsc-master 2016-12-08
Report Typos and Errors
  1: !
  2: !   Description: Solves a tridiagonal linear system with KSP.
  3: !
  4: !/*T
  5: !   Concepts: KSP^solving a system of linear equations
  6: !   Processors: 1
  7: !T*/
  8: ! -----------------------------------------------------------------------

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

 15: !
 16: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 17: !                   Variable declarations
 18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19: !
 20: !  Variables:
 21: !     ksp     - linear solver context
 22: !     ksp      - Krylov subspace method context
 23: !     pc       - preconditioner context
 24: !     x, b, u  - approx solution, right-hand-side, exact solution vectors
 25: !     A        - matrix that defines linear system
 26: !     its      - iterations for convergence
 27: !     norm     - norm of error in solution
 28: !
 29:       Vec              x,b,u
 30:       Mat              A
 31:       KSP              ksp
 32:       PC               pc
 33:       PetscReal        norm,tol
 34:       PetscErrorCode   ierr
 35:       PetscInt i,n,col(3),its,i1,i2,i3
 36:       PetscBool  flg
 37:       PetscMPIInt size,rank
 38:       PetscScalar      none,one,value(3)

 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !                 Beginning of program
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 44:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 45:       if (ierr .ne. 0) then
 46:         print*,'Unable to initialize PETSc'
 47:         stop
 48:       endif
 49:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 50:       if (size .ne. 1) then
 51:          call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 52:          if (rank .eq. 0) then
 53:             write(6,*) 'This is a uniprocessor example only!'
 54:          endif
 55:          SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
 56:       endif
 57:       none = -1.0
 58:       one  = 1.0
 59:       n    = 10
 60:       i1 = 1
 61:       i2 = 2
 62:       i3 = 3
 63:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,         &
 64:      &                        '-n',n,flg,ierr)

 66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67: !         Compute the matrix and right-hand-side vector that define
 68: !         the linear system, Ax = b.
 69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 71: !  Create matrix.  When using MatCreate(), the matrix format can
 72: !  be specified at runtime.

 74:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 75:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
 76:       call MatSetFromOptions(A,ierr)
 77:       call MatSetUp(A,ierr)

 79: !  Assemble matrix.
 80: !   - Note that MatSetValues() uses 0-based row and column numbers
 81: !     in Fortran as well as in C (as set here in the array "col").

 83:       value(1) = -1.0
 84:       value(2) = 2.0
 85:       value(3) = -1.0
 86:       do 50 i=1,n-2
 87:          col(1) = i-1
 88:          col(2) = i
 89:          col(3) = i+1
 90:          call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
 91:   50  continue
 92:       i = n - 1
 93:       col(1) = n - 2
 94:       col(2) = n - 1
 95:       call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
 96:       i = 0
 97:       col(1) = 0
 98:       col(2) = 1
 99:       value(1) = 2.0
100:       value(2) = -1.0
101:       call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
102:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
103:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

105: !  Create vectors.  Note that we form 1 vector from scratch and
106: !  then duplicate as needed.

108:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
109:       call VecSetSizes(x,PETSC_DECIDE,n,ierr)
110:       call VecSetFromOptions(x,ierr)
111:       call VecDuplicate(x,b,ierr)
112:       call VecDuplicate(x,u,ierr)

114: !  Set exact solution; then compute right-hand-side vector.

116:       call VecSet(u,one,ierr)
117:       call MatMult(A,u,b,ierr)

119: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: !          Create the linear solver and set various options
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

123: !  Create linear solver context

125:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

127: !  Set operators. Here the matrix that defines the linear system
128: !  also serves as the preconditioning matrix.

130:       call KSPSetOperators(ksp,A,A,ierr)

132: !  Set linear solver defaults for this problem (optional).
133: !   - By extracting the KSP and PC contexts from the KSP context,
134: !     we can then directly directly call any KSP and PC routines
135: !     to set various options.
136: !   - The following four statements are optional; all of these
137: !     parameters could alternatively be specified at runtime via
138: !     KSPSetFromOptions();

140:       call KSPGetPC(ksp,pc,ierr)
141:       call PCSetType(pc,PCJACOBI,ierr)
142:       tol = .0000001
143:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                         &
144:      &     PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)

146: !  Set runtime options, e.g.,
147: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
148: !  These options will override those specified above as long as
149: !  KSPSetFromOptions() is called _after_ any other customization
150: !  routines.

152:       call KSPSetFromOptions(ksp,ierr)

154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: !                      Solve the linear system
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

158:       call KSPSolve(ksp,b,x,ierr)

160: !  View solver info; we could instead use the option -ksp_view

162:       call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)

164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: !                      Check solution and clean up
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

168: !  Check the error

170:       call VecAXPY(x,none,u,ierr)
171:       call VecNorm(x,NORM_2,norm,ierr)
172:       call KSPGetIterationNumber(ksp,its,ierr)
173:       if (norm .gt. 1.e-12) then
174:         write(6,100) norm,its
175:       else
176:         write(6,200) its
177:       endif
178:  100  format('Norm of error ',e11.4,',  Iterations = ',i5)
179:  200  format('Norm of error < 1.e-12, Iterations = ',i5)

181: !  Free work space.  All PETSc objects should be destroyed when they
182: !  are no longer needed.

184:       call VecDestroy(x,ierr)
185:       call VecDestroy(u,ierr)
186:       call VecDestroy(b,ierr)
187:       call MatDestroy(A,ierr)
188:       call KSPDestroy(ksp,ierr)
189:       call PetscFinalize(ierr)

191:       end