Actual source code: ex1f.F90

  1: !
  2: !   Description: Solves a tridiagonal linear system with KSP.
  3: !
  4: ! -----------------------------------------------------------------------

  6:       program main
  7: #include <petsc/finclude/petscksp.h>
  8:       use petscksp
  9:       implicit none

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

 37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 38: !                 Beginning of program
 39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

 67:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 68:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
 69:       call MatSetFromOptions(A,ierr)
 70:       call MatSetUp(A,ierr)

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

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

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

101:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
102:       call VecSetSizes(x,PETSC_DECIDE,n,ierr)
103:       call VecSetFromOptions(x,ierr)
104:       call VecDuplicate(x,b,ierr)
105:       call VecDuplicate(x,u,ierr)

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

109:       call VecSet(u,one,ierr)
110:       call MatMult(A,u,b,ierr)
111:       call PetscLogStagePop(ierr)
112:       call PetscLogStagePush(stages(2),ierr)

114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: !          Create the linear solver and set various options
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

118: !  Create linear solver context

120:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

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

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

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

135:       call KSPGetPC(ksp,pc,ierr)
136:       call PCSetType(pc,PCJACOBI,ierr)
137:       tol = .0000001
138:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                         &
139:      &     PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)

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

147:       call KSPSetFromOptions(ksp,ierr)

149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !                      Solve the linear system
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

153:       call KSPSolve(ksp,b,x,ierr)
154:       call PetscLogStagePop(ierr)

156: !  View solver converged reason; we could instead use the option -ksp_converged_reason
157:       call KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)

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

161:       call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)

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

167: !  Check the error

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

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

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

190:       end

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