Actual source code: ex1f.F

  1: !
  2: !  Program usage: mpirun -np n ex1f [-help] [-n <n>] [all SLEPc options]
  3: !
  4: !  Description: Simple example that solves an eigensystem with the EPS object.
  5: !  The standard symmetric eigenvalue problem to be solved corresponds to the
  6: !  Laplacian operator in 1 dimension.
  7: !
  8: !  The command line options are:
  9: !    -n <n>, where <n> = number of grid points = matrix size
 10: !
 11: ! ----------------------------------------------------------------------
 12: !
 13:       program main
 14:       implicit none

 16: #include "finclude/petsc.h"
 17: #include "finclude/petscvec.h"
 18: #include "finclude/petscmat.h"
 19:  #include finclude/slepc.h
 20:  #include finclude/slepceps.h

 22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 23: !     Declarations
 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 25: !
 26: !  Variables:
 27: !     A     operator matrix
 28: !     eps   eigenproblem solver context

 30:       Mat          A
 31:       EPS          eps
 32:       EPSType      type
 33:       PetscReal    tol, error
 34:       PetscScalar  kr, ki
 35:       integer      rank, n, nev, ierr, maxit, i, its, nconv
 36:       integer      col(3), Istart, Iend
 37:       PetscTruth   flg
 38:       PetscScalar  value(3)

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

 44:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 45:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 46:       n = 30
 47:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)

 49:       if (rank .eq. 0) then
 50:         write(*,100) n
 51:       endif
 52:  100  format (/'1-D Laplacian Eigenproblem, n =',I3,' (Fortran)')

 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55: !     Compute the operator matrix that defines the eigensystem, Ax=kx
 56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 58:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 59:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
 60:       call MatSetFromOptions(A,ierr)

 62:       call MatGetOwnershipRange(A,Istart,Iend,ierr)
 63:       if (Istart .eq. 0) then
 64:         i = 0
 65:         col(1) = 0
 66:         col(2) = 1
 67:         value(1) =  2.0
 68:         value(2) = -1.0
 69:         call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
 70:         Istart = Istart+1
 71:       endif
 72:       if (Iend .eq. n) then
 73:         i = n-1
 74:         col(1) = n-2
 75:         col(2) = n-1
 76:         value(1) = -1.0
 77:         value(2) =  2.0
 78:         call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
 79:         Iend = Iend-1
 80:       endif
 81:       value(1) = -1.0
 82:       value(2) =  2.0
 83:       value(3) = -1.0
 84:       do i=Istart,Iend-1
 85:         col(1) = i-1
 86:         col(2) = i
 87:         col(3) = i+1
 88:         call MatSetValues(A,1,i,3,col,value,INSERT_VALUES,ierr)
 89:       enddo

 91:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 92:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95: !     Create the eigensolver and display info
 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 98: !     ** Create eigensolver context
 99:       call EPSCreate(PETSC_COMM_WORLD,eps,ierr)

101: !     ** Set operators. In this case, it is a standard eigenvalue problem
102:       call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
103:       call EPSSetProblemType(eps,EPS_HEP,ierr)

105: !     ** Set solver parameters at runtime
106:       call EPSSetFromOptions(eps,ierr)

108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: !     Solve the eigensystem
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

112:       call EPSSolve(eps,ierr)
113:       call EPSGetIterationNumber(eps,its,ierr)
114:       if (rank .eq. 0) then
115:         write(*,110) its
116:       endif
117:  110  format (/' Number of iterations of the method:',I4)
118: 
119: !     ** Optional: Get some information from the solver and display it
120:       call EPSGetType(eps,type,ierr)
121:       if (rank .eq. 0) then
122:         write(*,120) type
123:       endif
124:  120  format (' Solution method: ',A)
125:       call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,ierr)
126:       if (rank .eq. 0) then
127:         write(*,130) nev
128:       endif
129:  130  format (' Number of requested eigenvalues:',I2)
130:       call EPSGetTolerances(eps,tol,maxit,ierr)
131:       if (rank .eq. 0) then
132:         write(*,140) tol, maxit
133:       endif
134:  140  format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)

136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: !     Display solution and clean up
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

140: !     ** Get number of converged eigenpairs
141:       call EPSGetConverged(eps,nconv,ierr)
142:       if (rank .eq. 0) then
143:         write(*,150) nconv
144:       endif
145:  150  format (' Number of converged eigenpairs:',I2/)

147: !     ** Display eigenvalues and relative errors
148:       if (nconv.gt.0) then
149:         if (rank .eq. 0) then
150:           write(*,*) '         k          ||Ax-kx||/||kx||'
151:           write(*,*) ' ----------------- ------------------'
152:         endif
153:         do i=0,nconv-1
154: !         ** Get converged eigenpairs: i-th eigenvalue is stored in kr
155: !         ** (real part) and ki (imaginary part)
156:           call EPSGetEigenpair(eps,i,kr,ki,PETSC_NULL,PETSC_NULL,ierr)

158: !         ** Compute the relative error associated to each eigenpair
159:           call EPSComputeRelativeError(eps,i,error,ierr)
160:           if (rank .eq. 0) then
161:             write(*,160) PetscRealPart(kr), error
162:           endif
163:  160      format (1P,'   ',E12.4,'       ',E12.4)

165:         enddo
166:         if (rank .eq. 0) then
167:           write(*,*)
168:         endif
169:       endif

171: !     ** Free work space
172:       call EPSDestroy(eps,ierr)
173:       call MatDestroy(A,ierr)

175:       call SlepcFinalize(ierr)
176:       end