Actual source code: ex1f.F
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
5: !
6: ! This file is part of SLEPc.
7: ! SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: !
10: ! Program usage: mpiexec -n <np> ./ex1f [-help] [-n <n>] [all SLEPc options]
11: !
12: ! Description: Simple example that solves an eigensystem with the EPS object.
13: ! The standard symmetric eigenvalue problem to be solved corresponds to the
14: ! Laplacian operator in 1 dimension.
15: !
16: ! The command line options are:
17: ! -n <n>, where <n> = number of grid points = matrix size
18: !
19: ! ----------------------------------------------------------------------
20: !
21: program main
22: #include <slepc/finclude/slepceps.h>
23: use slepceps
24: implicit none
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! Declarations
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: !
30: ! Variables:
31: ! A operator matrix
32: ! eps eigenproblem solver context
34: Mat A
35: EPS eps
36: EPSType tname
37: PetscReal tol, error
38: PetscScalar kr, ki
39: Vec xr, xi
40: PetscInt n, i, Istart, Iend
41: PetscInt nev, maxit, its, nconv
42: PetscInt col(3)
43: PetscInt i1,i2,i3
44: PetscMPIInt rank
45: PetscErrorCode ierr
46: PetscBool flg
47: PetscScalar value(3)
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Beginning of program
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
54: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
55: n = 30
56: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
57: & '-n',n,flg,ierr)
59: if (rank .eq. 0) then
60: write(*,100) n
61: endif
62: 100 format (/'1-D Laplacian Eigenproblem, n =',I3,' (Fortran)')
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: ! Compute the operator matrix that defines the eigensystem, Ax=kx
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: call MatCreate(PETSC_COMM_WORLD,A,ierr)
69: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
70: call MatSetFromOptions(A,ierr)
71: call MatSetUp(A,ierr)
73: i1 = 1
74: i2 = 2
75: i3 = 3
76: call MatGetOwnershipRange(A,Istart,Iend,ierr)
77: if (Istart .eq. 0) then
78: i = 0
79: col(1) = 0
80: col(2) = 1
81: value(1) = 2.0
82: value(2) = -1.0
83: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
84: Istart = Istart+1
85: endif
86: if (Iend .eq. n) then
87: i = n-1
88: col(1) = n-2
89: col(2) = n-1
90: value(1) = -1.0
91: value(2) = 2.0
92: call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)
93: Iend = Iend-1
94: endif
95: value(1) = -1.0
96: value(2) = 2.0
97: value(3) = -1.0
98: do i=Istart,Iend-1
99: col(1) = i-1
100: col(2) = i
101: col(3) = i+1
102: call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)
103: enddo
105: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
106: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
108: call MatCreateVecs(A,xr,xi,ierr)
110: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: ! Create the eigensolver and display info
112: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! ** Create eigensolver context
115: call EPSCreate(PETSC_COMM_WORLD,eps,ierr)
117: ! ** Set operators. In this case, it is a standard eigenvalue problem
118: call EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr)
119: call EPSSetProblemType(eps,EPS_HEP,ierr)
121: ! ** Set solver parameters at runtime
122: call EPSSetFromOptions(eps,ierr)
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Solve the eigensystem
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: call EPSSolve(eps,ierr)
129: call EPSGetIterationNumber(eps,its,ierr)
130: if (rank .eq. 0) then
131: write(*,110) its
132: endif
133: 110 format (/' Number of iterations of the method:',I4)
135: ! ** Optional: Get some information from the solver and display it
136: call EPSGetType(eps,tname,ierr)
137: if (rank .eq. 0) then
138: write(*,120) tname
139: endif
140: 120 format (' Solution method: ',A)
141: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER, &
142: & PETSC_NULL_INTEGER,ierr)
143: if (rank .eq. 0) then
144: write(*,130) nev
145: endif
146: 130 format (' Number of requested eigenvalues:',I2)
147: call EPSGetTolerances(eps,tol,maxit,ierr)
148: if (rank .eq. 0) then
149: write(*,140) tol, maxit
150: endif
151: 140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Display solution and clean up
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! ** Get number of converged eigenpairs
158: call EPSGetConverged(eps,nconv,ierr)
159: if (rank .eq. 0) then
160: write(*,150) nconv
161: endif
162: 150 format (' Number of converged eigenpairs:',I2/)
164: ! ** Display eigenvalues and relative errors
165: if (nconv.gt.0) then
166: if (rank .eq. 0) then
167: write(*,*) ' k ||Ax-kx||/||kx||'
168: write(*,*) ' ----------------- ------------------'
169: endif
170: do i=0,nconv-1
171: ! ** Get converged eigenpairs: i-th eigenvalue is stored in kr
172: ! ** (real part) and ki (imaginary part)
173: call EPSGetEigenpair(eps,i,kr,ki,xr,xi,ierr)
175: ! ** Compute the relative error associated to each eigenpair
176: call EPSComputeError(eps,i,EPS_ERROR_RELATIVE,error,ierr)
177: if (rank .eq. 0) then
178: write(*,160) PetscRealPart(kr), error
179: endif
180: 160 format (1P,' ',E12.4,' ',E12.4)
182: enddo
183: if (rank .eq. 0) then
184: write(*,*)
185: endif
186: endif
188: ! ** Free work space
189: call EPSDestroy(eps,ierr)
190: call MatDestroy(A,ierr)
191: call VecDestroy(xr,ierr)
192: call VecDestroy(xi,ierr)
194: call SlepcFinalize(ierr)
195: end