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