-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy patheigen.f90
126 lines (91 loc) · 3.59 KB
/
eigen.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
! complex version on cpu
! a subroutine to calculate eigenvector and eigenvalue
subroutine eigensystem_c(JOBZ,UPLO,N,A,W)
use para, only : Dp
implicit none
! JOBZ (input) CHARACTER*1
! = 'N': Compute eigenvalues only;
! = 'V': Compute eigenvalues and eigenvectors.
!
character*1, intent(in) :: JOBZ
! UPLO (input) CHARACTER*1
! = 'U': Upper triangle of A is stored;
! = 'L': Lower triangle of A is stored.
character*1, intent(in) :: UPLO
! N (input) INTEGER
! The order of the matrix A. N >= 0.
integer, intent(in) :: N
! A (input/output) COMPLEX*16 array, dimension (LDA, N)
! On entry, the Hermitian matrix A. If UPLO = 'U', the
! leading N-by-N upper triangular part of A contains the
! upper triangular part of the matrix A. If UPLO = 'L',
! the leading N-by-N lower triangular part of A contains
! the lower triangular part of the matrix A.
! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
! orthonormal eigenvectors of the matrix A.
! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
! or the upper triangle (if UPLO='U') of A, including the
! diagonal, is destroyed.
complex(Dp),intent(inout) :: A(N,N)
! W (output) DOUBLE PRECISION array, dimension (N)
! If INFO = 0, the eigenvalues in ascending order.
real(Dp), intent(inout) :: W(N)
integer :: info
integer :: lwork
real(Dp),allocatable :: rwork(:)
complex(Dp),allocatable :: work(:)
lwork= 16*N
allocate(rwork(lwork))
allocate( work(lwork))
W= 0.0d0
info= 0
call zheev( JOBZ, UPLO, N, A, N, &
W, work, lwork, rwork, info )
return
end subroutine
! complex version on cpu
! a subroutine to calculate eigenvector and eigenvalue
subroutine eigensystem_zge(JOBZ, N, A, W)
use para, only : Dp
implicit none
! JOBZ (input) CHARACTER*1
! = 'N': Compute eigenvalues only;
! = 'V': Compute eigenvalues and eigenvectors.
!
character*1, intent(in) :: JOBZ
! N (input) INTEGER
! The order of the matrix A. N >= 0.
integer, intent(in) :: N
! A (input/output) COMPLEX*16 array, dimension (LDA, N)
! On entry, the Hermitian matrix A. If UPLO = 'U', the
! leading N-by-N upper triangular part of A contains the
! upper triangular part of the matrix A. If UPLO = 'L',
! the leading N-by-N lower triangular part of A contains
! the lower triangular part of the matrix A.
! On exit, if JOBZ = 'V', then if INFO = 0, A contains the
! orthonormal eigenvectors of the matrix A.
! If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
! or the upper triangle (if UPLO='U') of A, including the
! diagonal, is destroyed.
complex(Dp),intent(inout) :: A(N,N)
! W (output) DOUBLE PRECISION array, dimension (N)
! If INFO = 0, the eigenvalues in ascending order.
complex(Dp), intent(inout) :: W(N)
integer :: info
integer :: lwork
real(Dp),allocatable :: rwork(:)
complex(Dp),allocatable :: work(:)
complex(Dp),allocatable :: VL(:, :)
complex(Dp),allocatable :: VR(:, :)
lwork= 16*N
allocate(rwork(lwork))
allocate( work(lwork))
allocate( VL(N, N))
allocate( VR(N, N))
info=0
W=0.0d0
call zgeev(JOBZ, 'N', N, A, N, &
W, VL, N, VR, N, work, lwork, rwork, info )
A= VR
return
end subroutine