Actual source code: ex1f90.F90
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2009, Universidad 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 ex1f90 [-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
33: #include "finclude/petscdef.h"
34: #include finclude/slepcepsdef.h
35: use slepceps
37: implicit none
39: ! For usage without modules, uncomment the following lines and remove
40: ! the previous lines between 'program main' and 'implicit none'
41: !
42: !#include "finclude/petsc.h"
43: !#include "finclude/petscvec.h"
44: !#include "finclude/petscvec.h90"
45: !#include "finclude/petscmat.h"
46: !#include "finclude/petscmat.h90"
47: !#include "finclude/slepc.h"
48: !#include "finclude/slepceps.h"
49: !#include "finclude/slepceps.h90"
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: ! Declarations
53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: !
55: ! Variables:
56: ! A operator matrix
57: ! solver eigenproblem solver context
59: #if defined(PETSC_USE_FORTRAN_DATATYPES)
60: type(Mat) A
61: type(EPS) solver
62: #else
63: Mat A
64: EPS solver
65: #endif
66: EPSType tname
67: PetscReal tol, error
68: PetscScalar kr, ki
69: PetscInt n, i, Istart, Iend
70: PetscInt nev, maxit, its, nconv
71: PetscInt row(1), col(3)
72: PetscMPIInt rank
73: PetscErrorCode ierr
74: PetscTruth flg
75: PetscScalar value(3)
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: ! Beginning of program
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
82: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
83: n = 30
84: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
86: if (rank .eq. 0) then
87: write(*,100) n
88: endif
89: 100 format (/'1-D Laplacian Eigenproblem, n =',I4,' (Fortran)')
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Compute the operator matrix that defines the eigensystem, Ax=kx
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: call MatCreate(PETSC_COMM_WORLD,A,ierr)
96: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
97: call MatSetFromOptions(A,ierr)
99: call MatGetOwnershipRange(A,Istart,Iend,ierr)
100: if (Istart .eq. 0) then
101: row(1) = 0
102: col(1) = 0
103: col(2) = 1
104: value(1) = 2.0
105: value(2) = -1.0
106: call MatSetValues(A,1,row,2,col,value,INSERT_VALUES,ierr)
107: Istart = Istart+1
108: endif
109: if (Iend .eq. n) then
110: row(1) = n-1
111: col(1) = n-2
112: col(2) = n-1
113: value(1) = -1.0
114: value(2) = 2.0
115: call MatSetValues(A,1,row,2,col,value,INSERT_VALUES,ierr)
116: Iend = Iend-1
117: endif
118: value(1) = -1.0
119: value(2) = 2.0
120: value(3) = -1.0
121: do i=Istart,Iend-1
122: row(1) = i
123: col(1) = i-1
124: col(2) = i
125: col(3) = i+1
126: call MatSetValues(A,1,row,3,col,value,INSERT_VALUES,ierr)
127: enddo
129: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
130: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! Create the eigensolver and display info
134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! ** Create eigensolver context
137: call EPSCreate(PETSC_COMM_WORLD,solver,ierr)
139: ! ** Set operators. In this case, it is a standard eigenvalue problem
140: call EPSSetOperators(solver,A,PETSC_NULL_OBJECT,ierr)
141: call EPSSetProblemType(solver,EPS_HEP,ierr)
143: ! ** Set solver parameters at runtime
144: call EPSSetFromOptions(solver,ierr)
146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: ! Solve the eigensystem
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: call EPSSolve(solver,ierr)
151: call EPSGetIterationNumber(solver,its,ierr)
152: if (rank .eq. 0) then
153: write(*,110) its
154: endif
155: 110 format (/' Number of iterations of the method:',I4)
156:
157: ! ** Optional: Get some information from the solver and display it
158: call EPSGetType(solver,tname,ierr)
159: if (rank .eq. 0) then
160: write(*,120) tname
161: endif
162: 120 format (' Solution method: ',A)
163: call EPSGetDimensions(solver,nev,PETSC_NULL_INTEGER, &
164: PETSC_NULL_INTEGER,ierr)
165: if (rank .eq. 0) then
166: write(*,130) nev
167: endif
168: 130 format (' Number of requested eigenvalues:',I4)
169: call EPSGetTolerances(solver,tol,maxit,ierr)
170: if (rank .eq. 0) then
171: write(*,140) tol, maxit
172: endif
173: 140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)
175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: ! Display solution and clean up
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: ! ** Get number of converged eigenpairs
180: call EPSGetConverged(solver,nconv,ierr)
181: if (rank .eq. 0) then
182: write(*,150) nconv
183: endif
184: 150 format (' Number of converged eigenpairs:',I4/)
186: ! ** Display eigenvalues and relative errors
187: if (nconv.gt.0) then
188: if (rank .eq. 0) then
189: write(*,*) ' k ||Ax-kx||/||kx||'
190: write(*,*) ' ----------------- ------------------'
191: endif
192: do i=0,nconv-1
193: ! ** Get converged eigenpairs: i-th eigenvalue is stored in kr
194: ! ** (real part) and ki (imaginary part)
195: call EPSGetEigenpair(solver,i,kr,ki,PETSC_NULL,PETSC_NULL,ierr)
197: ! ** Compute the relative error associated to each eigenpair
198: call EPSComputeRelativeError(solver,i,error,ierr)
199: if (rank .eq. 0) then
200: write(*,160) PetscRealPart(kr), error
201: endif
202: 160 format (1P,' ',E12.4,' ',E12.4)
204: enddo
205: if (rank .eq. 0) then
206: write(*,*)
207: endif
208: endif
210: ! ** Free work space
211: call EPSDestroy(solver,ierr)
212: call MatDestroy(A,ierr)
214: call SlepcFinalize(ierr)
215: end