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:       PetscCallA(PetscInitialize(ierr))
 42:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 43:       if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
 44:       none = -1.0
 45:       one  = 1.0
 46:       n    = 10
 47:       i1 = 1
 48:       i2 = 2
 49:       i3 = 3
 50:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))

 52:       PetscCallA(PetscLogStageRegister("MatVec Assembly",stages(1),ierr))
 53:       PetscCallA(PetscLogStageRegister("KSP Solve",stages(2),ierr))
 54:       PetscCallA(PetscLogStagePush(stages(1),ierr))
 55: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56: !         Compute the matrix and right-hand-side vector that define
 57: !         the linear system, Ax = b.
 58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

 63:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
 64:       PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
 65:       PetscCallA(MatSetFromOptions(A,ierr))
 66:       PetscCallA(MatSetUp(A,ierr))

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

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

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

 97:       PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
 98:       PetscCallA(VecSetSizes(x,PETSC_DECIDE,n,ierr))
 99:       PetscCallA(VecSetFromOptions(x,ierr))
100:       PetscCallA(VecDuplicate(x,b,ierr))
101:       PetscCallA(VecDuplicate(x,u,ierr))

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

105:       PetscCallA(VecSet(u,one,ierr))
106:       PetscCallA(MatMult(A,u,b,ierr))
107:       PetscCallA(PetscLogStagePop(ierr))
108:       PetscCallA(PetscLogStagePush(stages(2),ierr))

110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: !          Create the linear solver and set various options
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

114: !  Create linear solver context

116:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))

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

121:       PetscCallA(KSPSetOperators(ksp,A,A,ierr))

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

131:       PetscCallA(KSPGetPC(ksp,pc,ierr))
132:       PetscCallA(PCSetType(pc,PCJACOBI,ierr))
133:       tol = .0000001
134:       PetscCallA(KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr))

136: !  Set runtime options, e.g.,
137: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
138: !  These options will override those specified above as long as
139: !  KSPSetFromOptions() is called _after_ any other customization
140: !  routines.

142:       PetscCallA(KSPSetFromOptions(ksp,ierr))

144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: !                      Solve the linear system
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

148:       PetscCallA(KSPSolve(ksp,b,x,ierr))
149:       PetscCallA(PetscLogStagePop(ierr))

151: !  View solver converged reason; we could instead use the option -ksp_converged_reason
152:       PetscCallA(KSPConvergedReasonView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr))

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

156:       PetscCallA(KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr))

158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: !                      Check solution and clean up
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

162: !  Check the error

164:       PetscCallA(VecAXPY(x,none,u,ierr))
165:       PetscCallA(VecNorm(x,NORM_2,norm,ierr))
166:       PetscCallA(KSPGetIterationNumber(ksp,its,ierr))
167:       if (norm .gt. 1.e-12) then
168:         write(6,100) norm,its
169:       else
170:         write(6,200) its
171:       endif
172:  100  format('Norm of error ',e11.4,',  Iterations = ',i5)
173:  200  format('Norm of error < 1.e-12, Iterations = ',i5)

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

178:       PetscCallA(VecDestroy(x,ierr))
179:       PetscCallA(VecDestroy(u,ierr))
180:       PetscCallA(VecDestroy(b,ierr))
181:       PetscCallA(MatDestroy(A,ierr))
182:       PetscCallA(KSPDestroy(ksp,ierr))
183:       PetscCallA(PetscFinalize(ierr))

185:       end

187: !/*TEST
188: !
189: !     test:
190: !       args: -ksp_monitor_short
191: !
192: !TEST*/