Actual source code: ex1f.F90

petsc-3.11.2 2019-05-18
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
 38:       PetscScalar      none,one,value(3)
 39:       PetscLogStage    stages(2);

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

 45:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 46:       if (ierr .ne. 0) then
 47:         print*,'Unable to initialize PETSc'
 48:         stop
 49:       endif
 50:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 51:       if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,1,'This is a uniprocessor example only'); endif
 52:       none = -1.0
 53:       one  = 1.0
 54:       n    = 10
 55:       i1 = 1
 56:       i2 = 2
 57:       i3 = 3
 58:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)

 60:       call PetscLogStageRegister("MatVec Assembly",stages(1),ierr)
 61:       call PetscLogStageRegister("KSP Solve",stages(2),ierr)
 62:       call PetscLogStagePush(stages(1),ierr)
 63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64: !         Compute the matrix and right-hand-side vector that define
 65: !         the linear system, Ax = b.
 66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

 71:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 72:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
 73:       call MatSetFromOptions(A,ierr)
 74:       call MatSetUp(A,ierr)

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

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

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

105:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
106:       call VecSetSizes(x,PETSC_DECIDE,n,ierr)
107:       call VecSetFromOptions(x,ierr)
108:       call VecDuplicate(x,b,ierr)
109:       call VecDuplicate(x,u,ierr)

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

113:       call VecSet(u,one,ierr)
114:       call MatMult(A,u,b,ierr)
115:       call PetscLogStagePop(ierr)
116:       call PetscLogStagePush(stages(2),ierr)

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

122: !  Create linear solver context

124:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

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

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

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

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

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

151:       call KSPSetFromOptions(ksp,ierr)

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

157:       call KSPSolve(ksp,b,x,ierr)
158:       call PetscLogStagePop(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

193: !/*TEST
194: !
195: !     test:
196: !       args: -ksp_monitor_short
197: !
198: !TEST*/