/
EigValSym.F95
160 lines (134 loc) · 4.57 KB
/
EigValSym.F95
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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
subroutine EigValSym(ain, n, eval, ul)
!-------------------------------------------------------------------------------
!
! This subroutine will return the eigenvalues of the symmetric square
! matrix Ain, ordered from greatest to least.
!
! Calling Parameters
!
! IN
! Ain Input symmetric matrix. By default, only the
! upper portion is used.
! n Order of the matrix Ain.
!
! OUT
! eval Vector of length n of the eigenvalues of Ain.
!
! OPTIONAL
! uplo Use the upper 'U' or lower 'L' portion of the
! input symmetric matrix.
!
! The eigenvalues and eigenvectors are determined by reducing the
! matrix to
!
! A = Z L Z = Q (S L S') Q'
!
! by the two operations:
!
! (1) The real symmetric square matrix is reduced to tridiagonal form
! A = Q T Q'
! where Q is orthogonal, and T is symmetric tridiagonal.
! (2) The tridiagonal matrix is reduced to
! T = S L S'
!
! The eigenvalues of A correspond to L (which is a diagonal).
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!-------------------------------------------------------------------------------
use ftypes
implicit none
real(dp), intent(in) :: ain(:,:)
integer, intent(in) :: n
real(dp), intent(out) :: eval(:)
character, intent(in), optional :: ul
integer, parameter :: nb = 80, nbl = 10
character :: uplo
real(dp) :: d(n), e(n), tau(n-1), work(nb*n), vl, vu, abstol, w(n)
real(dp), allocatable :: a(:,:), z(:,:)
integer :: lwork, info, il, iu, m, isuppz(2*n), liwork, iwork(nbl*n), &
i, astat(2)
#ifdef LAPACK_UNDERSCORE
#define dsytrd dsytrd_
#define dstegr dstegr_
#endif
external dsytrd, dstegr
if (size(ain(:,1)) < n .or. size(ain(1,:)) < n) then
print*, "Error --- EigValSym"
print*, "AIN must be dimensioned as (N, N) where N is ", n
print*, "Input array is dimensioned as ", size(ain(:,1)), size(ain(1,:))
stop
else if (size(eval) < n) then
print*, "Error --- EigValSym"
print*, "EVAL must be dimensioned as (N) where N is ", n
print*, "Input array is dimensioned as ", size(eval)
stop
end if
allocate (a(n,n), stat = astat(1))
allocate (z(n,n), stat = astat(2))
if (astat(1) /= 0 .or. astat(2) /= 0) then
print*, "Error --- EigValSym"
print*, "Problem allocating arrays A and Z", astat(1), astat(2)
stop
end if
lwork = nb*n
liwork = nbl*n
eval = 0.0_dp
a(1:n,1:n) = ain(1:n,1:n)
if (present(ul)) then
uplo = ul
else
uplo = "U"
end if
!---------------------------------------------------------------------------
!
! Factor A = Q T Q'
!
!---------------------------------------------------------------------------
call dsytrd(uplo, n, a, n, d, e(1:n-1), tau, work, lwork, info)
if (info /= 0) then
print*, "Error --- EigValSym"
print*, "Problem tri-diagonalizing input matrix"
stop
else
if ( work(1) > dble(lwork) ) then
print*, "Warning --- EigValSym"
print*, "Consider changing value of nb to ", work(1)/n, &
" and recompile the SHTOOLS archive."
end if
end if
!---------------------------------------------------------------------------
!
! Factor T = S L S'
!
!---------------------------------------------------------------------------
abstol = 0.0_dp
call dstegr('n','a', n, d, e, vl, vu, il, iu, abstol, m, w, &
z, n, isuppz, work, lwork, iwork, liwork, info)
if (info /= 0) then
print*, "Error --- EigValSym"
print*, "Problem determining eigenvalues and eigenvectors of " // &
"tridiagonal matrix."
if (info == 1) print*, "Internal error in DLARRE"
if (info == 2) print*, "Internal error in DLARRV"
stop
else
if (work(1) > dble(lwork) ) then
print*, "Warning --- EigValSym"
print*, "Consider changing value of nb to ", work(1)/n, &
" and recompile the SHTOOLS archive."
end if
if (iwork(1) > liwork ) then
print*, "Warning --- Eigsym"
print*, "Consider changing value of nb to ", iwork(1)/n, &
" and recompile the SHTOOLS archive."
end if
end if
! Reorder eigenvalues from greatest to least.
do i = 1, n
eval(i) = w(n+1-i)
end do
deallocate(a)
deallocate(z)
end subroutine EigValSym