Actual source code: ex1f.F
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain
4: !
5: ! This file is part of SLEPc.
6: !
7: ! SLEPc is free software: you can redistribute it and/or modify it under the
8: ! terms of version 3 of the GNU Lesser General Public License as published by
9: ! the Free Software Foundation.
10: !
11: ! SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
12: ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13: ! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14: ! more details.
15: !
16: ! You should have received a copy of the GNU Lesser General Public License
17: ! along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: !
20: ! Program usage: mpirun -np n ex1f [-help] [-n <n>] [all SLEPc options]
21: !
22: ! Description: Simple example that solves an eigensystem with the EPS object.
23: ! The standard symmetric eigenvalue problem to be solved corresponds to the
24: ! Laplacian operator in 1 dimension.
25: !
26: ! The command line options are:
27: ! -n <n>, where <n> = number of grid points = matrix size
28: !
29: ! ----------------------------------------------------------------------
30: !
31: program main
32: implicit none
34: #include <finclude/petscsys.h>
35: #include <finclude/petscvec.h>
36: #include <finclude/petscmat.h>
37: #include <finclude/slepcsys.h>
38: #include <finclude/slepceps.h>
40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: ! Declarations
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: !
44: ! Variables:
45: ! A operator matrix
46: ! eps eigenproblem solver context
48: Mat A
49: EPS eps
50: EPSType tname
51: PetscReal tol, error
52: PetscScalar kr, ki
53: PetscInt n, i, Istart, Iend
54: PetscInt nev, maxit, its, nconv
55: PetscInt col(3)
56: PetscInt i1,i2,i3
57: PetscMPIInt rank
58: PetscErrorCode ierr
59: PetscBool flg
60: PetscScalar value(3)
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: ! Beginning of program
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
67: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
68: n = 30
69: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
71: if (rank .eq. 0) then
72: write(*,100) n
73: endif
74: 100 format (/'1-D Laplacian Eigenproblem, n =',I3,' (Fortran)')
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! Compute the operator matrix that defines the eigensystem, Ax=kx
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: call MatCreate(PETSC_COMM_WORLD,A,ierr)
81: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
82: call MatSetFromOptions(A,ierr)
83: call MatSetUp(A,ierr)
85: i1 = 1
86: i2 = 2
87: i3 = 3
88: call MatGetOwnershipRange(A,Istart,Iend,ierr)
89: if (Istart .eq. 0) then
90: i = 0
91: col(1) = 0
92: col(2) = 1
93: value(1) = 2.0
94: value(2) = -1.0
95: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
96: Istart = Istart+1
97: endif
98: if (Iend .eq. n) then
99: i = n-1
100: col(1) = n-2
101: col(2) = n-1
102: value(1) = -1.0
103: value(2) = 2.0
104: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
105: Iend = Iend-1
106: endif
107: value(1) = -1.0
108: value(2) = 2.0
109: value(3) = -1.0
110: do i=Istart,Iend-1
111: col(1) = i-1
112: col(2) = i
113: col(3) = i+1
114: call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
115: enddo
117: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
118: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Create the eigensolver and display info
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! ** Create eigensolver context
125: call EPSCreate(PETSC_COMM_WORLD,eps,ierr)
127: ! ** Set operators. In this case, it is a standard eigenvalue problem
128: call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
129: call EPSSetProblemType(eps,EPS_HEP,ierr)
131: ! ** Set solver parameters at runtime
132: call EPSSetFromOptions(eps,ierr)
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: ! Solve the eigensystem
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: call EPSSolve(eps,ierr)
139: call EPSGetIterationNumber(eps,its,ierr)
140: if (rank .eq. 0) then
141: write(*,110) its
142: endif
143: 110 format (/' Number of iterations of the method:',I4)
145: ! ** Optional: Get some information from the solver and display it
146: call EPSGetType(eps,tname,ierr)
147: if (rank .eq. 0) then
148: write(*,120) tname
149: endif
150: 120 format (' Solution method: ',A)
151: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER, &
152: & PETSC_NULL_INTEGER,ierr)
153: if (rank .eq. 0) then
154: write(*,130) nev
155: endif
156: 130 format (' Number of requested eigenvalues:',I2)
157: call EPSGetTolerances(eps,tol,maxit,ierr)
158: if (rank .eq. 0) then
159: write(*,140) tol, maxit
160: endif
161: 140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Display solution and clean up
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: ! ** Get number of converged eigenpairs
168: call EPSGetConverged(eps,nconv,ierr)
169: if (rank .eq. 0) then
170: write(*,150) nconv
171: endif
172: 150 format (' Number of converged eigenpairs:',I2/)
174: ! ** Display eigenvalues and relative errors
175: if (nconv.gt.0) then
176: if (rank .eq. 0) then
177: write(*,*) ' k ||Ax-kx||/||kx||'
178: write(*,*) ' ----------------- ------------------'
179: endif
180: do i=0,nconv-1
181: ! ** Get converged eigenpairs: i-th eigenvalue is stored in kr
182: ! ** (real part) and ki (imaginary part)
183: call EPSGetEigenpair(eps,i,kr,ki,PETSC_NULL_OBJECT, &
184: & PETSC_NULL_OBJECT,ierr)
186: ! ** Compute the relative error associated to each eigenpair
187: call EPSComputeRelativeError(eps,i,error,ierr)
188: if (rank .eq. 0) then
189: write(*,160) PetscRealPart(kr), error
190: endif
191: 160 format (1P,' ',E12.4,' ',E12.4)
193: enddo
194: if (rank .eq. 0) then
195: write(*,*)
196: endif
197: endif
199: ! ** Free work space
200: call EPSDestroy(eps,ierr)
201: call MatDestroy(A,ierr)
203: call SlepcFinalize(ierr)
204: end