Actual source code: ex1f.F

petsc-3.7.3 2016-07-24
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:       implicit none

 13: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 14: !                    Include files
 15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 16: !
 17: !  This program uses CPP for preprocessing, as indicated by the use of
 18: !  PETSc include files in the directory petsc/include/petsc/finclude.  This
 19: !  convention enables use of the CPP preprocessor, which allows the use
 20: !  of the #include statements that define PETSc objects and variables.
 21: !
 22: !  Use of the conventional Fortran include statements is also supported
 23: !  In this case, the PETsc include files are located in the directory
 24: !  petsc/include/foldinclude.
 25: !
 26: !  Since one must be very careful to include each file no more than once
 27: !  in a Fortran routine, application programmers must exlicitly list
 28: !  each file needed for the various PETSc components within their
 29: !  program (unlike the C/C++ interface).
 30: !
 31: !  See the Fortran section of the PETSc users manual for details.
 32: !
 33: !  The following include statements are required for KSP Fortran programs:
 34: !     petscsys.h       - base PETSc routines
 35: !     petscvec.h    - vectors
 36: !     petscmat.h    - matrices
 37: !     petscksp.h    - Krylov subspace methods
 38: !     petscpc.h     - preconditioners
 39: !  Other include statements may be needed if using additional PETSc
 40: !  routines in a Fortran program, e.g.,
 41: !     petscviewer.h - viewers
 42: !     petscis.h     - index sets
 43: !
 44: #include <petsc/finclude/petscsys.h>
 45: #include <petsc/finclude/petscvec.h>
 46: #include <petsc/finclude/petscmat.h>
 47: #include <petsc/finclude/petscksp.h>
 48: #include <petsc/finclude/petscpc.h>
 49: !
 50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51: !                   Variable declarations
 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53: !
 54: !  Variables:
 55: !     ksp     - linear solver context
 56: !     ksp      - Krylov subspace method context
 57: !     pc       - preconditioner context
 58: !     x, b, u  - approx solution, right-hand-side, exact solution vectors
 59: !     A        - matrix that defines linear system
 60: !     its      - iterations for convergence
 61: !     norm     - norm of error in solution
 62: !
 63:       Vec              x,b,u
 64:       Mat              A
 65:       KSP              ksp
 66:       PC               pc
 67:       PetscReal        norm,tol
 68:       PetscErrorCode   ierr
 69:       PetscInt i,n,col(3),its,i1,i2,i3
 70:       PetscBool  flg
 71:       PetscMPIInt size,rank
 72:       PetscScalar      none,one,value(3)

 74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75: !                 Beginning of program
 76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 78:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 79:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 80:       if (size .ne. 1) then
 81:          call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 82:          if (rank .eq. 0) then
 83:             write(6,*) 'This is a uniprocessor example only!'
 84:          endif
 85:          SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
 86:       endif
 87:       none = -1.0
 88:       one  = 1.0
 89:       n    = 10
 90:       i1 = 1
 91:       i2 = 2
 92:       i3 = 3
 93:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,         &
 94:      &                        '-n',n,flg,ierr)

 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97: !         Compute the matrix and right-hand-side vector that define
 98: !         the linear system, Ax = b.
 99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

104:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
105:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
106:       call MatSetFromOptions(A,ierr)
107:       call MatSetUp(A,ierr)

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

113:       value(1) = -1.0
114:       value(2) = 2.0
115:       value(3) = -1.0
116:       do 50 i=1,n-2
117:          col(1) = i-1
118:          col(2) = i
119:          col(3) = i+1
120:          call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
121:   50  continue
122:       i = n - 1
123:       col(1) = n - 2
124:       col(2) = n - 1
125:       call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
126:       i = 0
127:       col(1) = 0
128:       col(2) = 1
129:       value(1) = 2.0
130:       value(2) = -1.0
131:       call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
132:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
133:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

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

138:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
139:       call VecSetSizes(x,PETSC_DECIDE,n,ierr)
140:       call VecSetFromOptions(x,ierr)
141:       call VecDuplicate(x,b,ierr)
142:       call VecDuplicate(x,u,ierr)

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

146:       call VecSet(u,one,ierr)
147:       call MatMult(A,u,b,ierr)

149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !          Create the linear solver and set various options
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

153: !  Create linear solver context

155:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

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

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

162: !  Set linear solver defaults for this problem (optional).
163: !   - By extracting the KSP and PC contexts from the KSP context,
164: !     we can then directly directly call any KSP and PC routines
165: !     to set various options.
166: !   - The following four statements are optional; all of these
167: !     parameters could alternatively be specified at runtime via
168: !     KSPSetFromOptions();

170:       call KSPGetPC(ksp,pc,ierr)
171:       call PCSetType(pc,PCJACOBI,ierr)
172:       tol = .0000001
173:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                         &
174:      &     PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)

176: !  Set runtime options, e.g.,
177: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
178: !  These options will override those specified above as long as
179: !  KSPSetFromOptions() is called _after_ any other customization
180: !  routines.

182:       call KSPSetFromOptions(ksp,ierr)

184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: !                      Solve the linear system
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

192:       call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)

194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: !                      Check solution and clean up
196: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

198: !  Check the error

200:       call VecAXPY(x,none,u,ierr)
201:       call VecNorm(x,NORM_2,norm,ierr)
202:       call KSPGetIterationNumber(ksp,its,ierr)
203:       if (norm .gt. 1.e-12) then
204:         write(6,100) norm,its
205:       else
206:         write(6,200) its
207:       endif
208:  100  format('Norm of error ',e11.4,',  Iterations = ',i5)
209:  200  format('Norm of error < 1.e-12, Iterations = ',i5)

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

214:       call VecDestroy(x,ierr)
215:       call VecDestroy(u,ierr)
216:       call VecDestroy(b,ierr)
217:       call MatDestroy(A,ierr)
218:       call KSPDestroy(ksp,ierr)
219:       call PetscFinalize(ierr)

221:       end