Actual source code: ex1f.F

  1: !
  2: !    "$Id: ex1f.F,v 1.30 2001/08/07 03:04:00 balay Exp $";
  3: !
  4: !   Description: Solves a tridiagonal linear system with KSP.
  5: !
  6: !/*T
  7: !   Concepts: KSP^solving a system of linear equations
  8: !   Processors: 1
  9: !T*/
 10: ! -----------------------------------------------------------------------

 12:       program main
 13:       implicit none

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

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

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

 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92: !         Compute the matrix and right-hand-side vector that define
 93: !         the linear system, Ax = b.
 94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

 99:       call MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,A,   &
100:      &               ierr)
101:       call MatSetFromOptions(A,ierr)

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

107:       value(1) = -1.0
108:       value(2) = 2.0
109:       value(3) = -1.0
110:       do 50 i=1,n-2
111:          col(1) = i-1
112:          col(2) = i
113:          col(3) = i+1
114:          call MatSetValues(A,1,i,3,col,value,INSERT_VALUES,ierr)
115:   50  continue
116:       i = n - 1
117:       col(1) = n - 2
118:       col(2) = n - 1
119:       call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
120:       i = 0
121:       col(1) = 0
122:       col(2) = 1
123:       value(1) = 2.0
124:       value(2) = -1.0
125:       call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
126:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
127:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

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

132:       call VecCreate(PETSC_COMM_WORLD,x,ierr)
133:       call VecSetSizes(x,PETSC_DECIDE,n,ierr)
134:       call VecSetFromOptions(x,ierr)
135:       call VecDuplicate(x,b,ierr)
136:       call VecDuplicate(x,u,ierr)

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

140:       call VecSet(one,u,ierr)
141:       call MatMult(A,u,b,ierr)

143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: !          Create the linear solver and set various options
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

147: !  Create linear solver context

149:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

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

154:       call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)

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

164:       call KSPGetPC(ksp,pc,ierr)
165:       call PCSetType(pc,PCJACOBI,ierr)
166:       tol = 1.d-7
167:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,     &
168:      &     PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)

170: !  Set runtime options, e.g.,
171: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
172: !  These options will override those specified above as long as
173: !  KSPSetFromOptions() is called _after_ any other customization
174: !  routines.

176:       call KSPSetFromOptions(ksp,ierr)

178: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: !                      Solve the linear system
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

182:       call KSPSetRhs(ksp,b,ierr)
183:       call KSPSetSolution(ksp,x,ierr)
184:       call KSPSolve(ksp,ierr)

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

188:       call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)

190: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: !                      Check solution and clean up
192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

194: !  Check the error

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

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

210:       call VecDestroy(x,ierr)
211:       call VecDestroy(u,ierr)
212:       call VecDestroy(b,ierr)
213:       call MatDestroy(A,ierr)
214:       call KSPDestroy(ksp,ierr)
215:       call PetscFinalize(ierr)

217:       end