-
Notifications
You must be signed in to change notification settings - Fork 1
/
eigen.f90
363 lines (309 loc) · 12.2 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
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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
!Eigen_solver is a F90 module which simplifies the eigen value sovling probelm for those
!first time users of lapack.
!author: Jiuzhou Tang ([email protected])
!history: 2017-8-8(first eition)
!Credits: this code is greatly inspired by the open source code(iQIST: An open
!source continuous-time quantum Monte Carlo impurity solver toolkit), with one
!of its author a close friend of mine.
!For detailed information of eigen solver, please refer to the lapack's user guide on line.
!this eigen module will follow matlab's eigen funtion style, i.e., [val,vec]=eigen(A), which means you can specify the output of eigen vectors as optional, which is necessary for many scientific applications.
!!!Very importmant notification:
! eigen values output are not ordered in procedure eigen_g_r/eigen_g_c
! eigen values output in Ascending order in procedure eigen_s_r/eigen_h_c
module eigen_solver
integer,parameter :: dp = kind(1.0d0) !double precision
interface eigen
module procedure eigen_g_r ! for general real matrix, by default, using double real
module procedure eigen_s_r ! for symmetric real matrix, by default,double complex
module procedure eigen_g_c ! for general complex matrix
module procedure eigen_h_c ! for hermite complex matrix
end interface
contains
subroutine eigen_g_r(A,val,vec,ierr)
implicit none
real(dp), intent(in) :: A(:,:) ! input matrix A
integer :: mat_dim ! matrix dimension, assuming square matrix
real(dp),intent(out) :: val(:) ! output array of eigen value
real(dp), allocatable, intent(out), optional :: vec(:,:) ! list of output arrays of eigen vector,
! stored in a matrix form. i.e., [[v1],[v2],...v[n]].
integer, intent(out), optional :: ierr ! error info, 0 for sucess, otherwise
! error, labeled with different integer.
integer :: dim(2) ! result of shape(A)
! lapack local variables
integer :: ndim !matrix dim used in lapack subroutine
integer :: istat ! status flag
integer :: info ! info of lapack computing subroutine,if info=0, the eigenvalues will
!output in descending order
integer :: lwork ! the length of the array work, lwork >= max(1,4*ndim)
real(dp), allocatable :: work(:) ! workspace array
real(dp), allocatable :: wr(:)
real(dp), allocatable :: wi(:) ! auxiliary real(dp) matrix: real and imaginary parts of eigenvalues
real(dp), allocatable :: vr(:,:)
real(dp), allocatable :: vl(:,:) ! auxiliary real(dp) matrix: left and right eigenvectors
real(dp), allocatable :: vec_tmp(:,:) ! a temporary place holder for input matrix and output eigen vec.
if (present(ierr)) ierr = 0 ! set ierr = 0 if present
dim = shape(A)
if (dim(1) .eq. dim(2)) then
mat_dim = dim(1)
else
if ( present(ierr)) ierr = 1
return
endif
! allocate memory
if(present(vec)) then
allocate(vec(mat_dim,mat_dim), stat=ierr)
else
allocate(vec_tmp(mat_dim,mat_dim), stat=ierr)
endif
ndim = mat_dim
! initialize lwork
lwork = 4 * ndim
allocate(work(lwork), stat=istat)
allocate(wr(ndim), stat=istat)
allocate(wi(ndim), stat=istat)
allocate(vr(ndim,ndim), stat=istat)
allocate(vl(ndim,ndim), stat=istat)
if ( istat /= 0 ) then
if ( present(ierr)) ierr = 2
print*, 'can not allocate enough memory'
return
endif
! call the computational subroutine: dgeev
if ( present(vec)) then
vec = A
call DGEEV('N', 'V', ndim, vec, ndim, wr, wi, vl, ndim, vr, ndim, work, lwork, info)
else
vec_tmp = A
call DGEEV('N', 'N', ndim, vec_tmp, ndim, wr, wi, vl, ndim, vr, ndim, work, lwork, info)
endif
! check the status
if ( info /= 0 ) then
if ( present(ierr)) ierr = 3
print*, 'error in lapack subroutine dgeev'
return
endif !
! copy eigenvalues and eigenvectors
val(1:ndim) = wr(1:ndim)
if(present(vec)) then
vec(1:ndim,1:ndim) = vr(1:ndim,1:ndim)
else
vec_tmp(1:ndim,1:ndim) = vr(1:ndim,1:ndim)
endif
! dealloate memory for workspace array
if ( allocated(work) ) deallocate(work)
if ( allocated(wr ) ) deallocate(wr )
if ( allocated(wi ) ) deallocate(wi )
if ( allocated(vr ) ) deallocate(vr )
if ( allocated(vl ) ) deallocate(vl )
if(.not. present(vec)) then
if ( allocated(vec_tmp ) ) deallocate(vec_tmp )
endif
end subroutine eigen_g_r
subroutine eigen_g_c(A,val,vec,ierr)
implicit none
complex(dp), intent(in) :: A(:,:) ! input matrix A
integer :: mat_dim ! matrix dimension, assuming square matrix
complex(dp),intent(out) :: val(:) ! output array of eigen value
complex(dp), allocatable, intent(out), optional :: vec(:,:) ! list of output arrays of eigen vector,
! stored in a matrix form. i.e., [[v1],[v2],...v[n]].
integer, intent(out), optional :: ierr ! error info, 0 for sucess, otherwise
! error, labeled with different integer.
integer :: dim(2) ! result of shape(A)
! lapack local variables
integer :: ndim !matrix dim used in lapack subroutine
integer :: istat ! status flag
integer :: info ! info of lapack computing subroutine,if info=0, the eigenvalues will
!output in descending order
integer :: lwork ! the length of the array work, lwork >= max(1,2*ndim)
complex(dp), allocatable :: work(:) ! workspace array
complex(dp), allocatable :: rwork(:)
complex(dp), allocatable :: vr(:,:)
complex(dp), allocatable :: vl(:,:) ! auxiliary matrix: left and right eigenvectors
complex(dp), allocatable :: vec_tmp(:,:) ! a temporary place holder for input matrix and output eigen vec.
if (present(ierr)) ierr = 0 ! set ierr = 0 if present
dim = shape(A)
if (dim(1) .eq. dim(2)) then
mat_dim = dim(1)
else
if ( present(ierr)) ierr = 1
return
endif
! allocate memory
if(present(vec)) then
allocate(vec(mat_dim,mat_dim), stat=ierr)
else
allocate(vec_tmp(mat_dim,mat_dim), stat=ierr)
endif
ndim = mat_dim
! initialize lwork
lwork = 2 * ndim
allocate(work(lwork), stat=istat)
allocate(rwork(ndim), stat=istat)
allocate(vr(ndim,ndim), stat=istat)
allocate(vl(ndim,ndim), stat=istat)
if ( istat /= 0 ) then
if ( present(ierr)) ierr = 2
print*, 'can not allocate enough memory'
return
endif
! call the computational subroutine: zgeev
if ( present(vec)) then
vec = A
call ZGEEV('N', 'V', ndim, vec, ndim,val, vl, ndim, vr, ndim, work, lwork,rwork, info)
else
vec_tmp = A
call ZGEEV('N', 'N', ndim, vec_tmp, ndim, val, vl, ndim, vr, ndim, work, lwork,rwork, info)
endif
! check the status
if ( info /= 0 ) then
if ( present(ierr)) ierr = 3
print*, 'error in lapack subroutine zgeev'
return
endif !
if(present(vec)) then
vec = vr
else
vec_tmp = vr
endif
! dealloate memory for workspace array
if ( allocated(work) ) deallocate(work)
if ( allocated(rwork ) ) deallocate(rwork )
if ( allocated(vr ) ) deallocate(vr )
if ( allocated(vl ) ) deallocate(vl )
if(.not. present(vec)) then
if ( allocated(vec_tmp ) ) deallocate(vec_tmp )
endif
end subroutine eigen_g_c
subroutine eigen_s_r(A,p,val,vec,ierr)
implicit none
real(dp), intent(in) :: A(:,:) ! input matrix A
character*1, intent(in) :: p ! using upper triangular or lower triangular of
!the symmetric matrix, 'U' for upper, 'L' for lower
integer :: mat_dim ! matrix dimension, assuming square matrix
real(dp),intent(out) :: val(:) ! output array of eigen value
real(dp), allocatable, intent(out), optional :: vec(:,:) ! list of output arrays of eigen vector,
! stored in a matrix form. i.e., [[v1],[v2],...v[n]].
integer, intent(out), optional :: ierr ! error info, 0 for sucess, otherwise
! error, labeled with different integer.
integer :: dim(2) ! result of shape(A)
! lapack local variables
integer :: ndim !matrix dim used in lapack subroutine
integer :: istat ! status flag
integer :: info ! info of lapack computing subroutine,if info=0, the eigenvalues will
!output in ascending order
integer :: lwork ! the length of the array work, lwork >= max(1,4*ndim)
real(dp), allocatable :: work(:) ! workspace array
real(dp), allocatable :: vec_tmp(:,:) ! a temporary place holder for input matrix and output eigen vec.
if (present(ierr)) ierr = 0 ! set ierr = 0 if present
dim = shape(A)
if (dim(1) .eq. dim(2)) then
mat_dim = dim(1)
else
if ( present(ierr)) ierr = 1
return
endif
! allocate memory
if(present(vec)) then
allocate(vec(mat_dim,mat_dim), stat=ierr)
else
allocate(vec_tmp(mat_dim,mat_dim), stat=ierr)
endif
ndim = mat_dim
! initialize lwork
lwork = 3 * ndim -1
allocate(work(lwork), stat=istat)
if ( istat /= 0 ) then
if ( present(ierr)) ierr = 2
print*, 'can not allocate enough memory'
return
endif
! call the computational subroutine: dgeev
if ( present(vec)) then
vec = A
call DSYEV('V', p, ndim, vec, ndim, val, work, lwork, info)
else
vec_tmp = A
call DSYEV('N', p, ndim, vec_tmp, ndim, val, work, lwork, info)
endif
! check the status
if ( info /= 0 ) then
if ( present(ierr)) ierr = 3
print*, 'error in lapack subroutine dsyev'
return
endif !
! dealloate memory for workspace array
if ( allocated(work) ) deallocate(work)
if(.not. present(vec)) then
if ( allocated(vec_tmp ) ) deallocate(vec_tmp )
endif
end subroutine eigen_s_r
subroutine eigen_h_c(A,p,val,vec,ierr)
implicit none
complex(dp), intent(in) :: A(:,:) ! input matrix A
character*1, intent(in) :: p ! using upper triangular or lower triangular of
!the symmetric matrix, 'U' for upper, 'L' for lower
integer :: mat_dim ! matrix dimension, assuming square matrix
complex(dp),intent(out) :: val(:) ! output array of eigen value
complex(dp), allocatable, intent(out), optional :: vec(:,:) ! list of output arrays of eigen vector,
! stored in a matrix form. i.e., [[v1],[v2],...v[n]].
integer, intent(out), optional :: ierr ! error info, 0 for sucess, otherwise
! error, labeled with different integer.
integer :: dim(2) ! result of shape(A)
! lapack local variables
integer :: ndim !matrix dim used in lapack subroutine
integer :: istat ! status flag
integer :: info ! info of lapack computing subroutine,if info=0, the eigenvalues will
!output in ascending order
integer :: lwork,lrwork ! the length of the array work
!! the length of the array work and rwork
! lwork >= max(1,2*ndim-1), lrwork >= max(1,3*ndim-2)
complex(dp), allocatable :: work(:) ! workspace array
complex(dp), allocatable :: rwork(:) ! workspace array
complex(dp), allocatable :: vec_tmp(:,:) ! a temporary place holder for input matrix and output eigen vec.
if (present(ierr)) ierr = 0 ! set ierr = 0 if present
dim = shape(A)
if (dim(1) .eq. dim(2)) then
mat_dim = dim(1)
else
if ( present(ierr)) ierr = 1
return
endif
! allocate memory
if(present(vec)) then
allocate(vec(mat_dim,mat_dim), stat=ierr)
else
allocate(vec_tmp(mat_dim,mat_dim), stat=ierr)
endif
ndim = mat_dim
! initialize lwork
lwork = 2 * ndim - 1
lrwork = 3 * ndim - 2
allocate(work(lwork), stat=istat)
allocate(rwork(lrwork), stat=istat)
if ( istat /= 0 ) then
if ( present(ierr)) ierr = 2
print*, 'can not allocate enough memory'
return
endif
! call the computational subroutine: dgeev
if ( present(vec)) then
vec = A
call ZHEEV('V', p, ndim, vec, ndim, val, work, lwork, rwork,info)
else
vec_tmp = A
call ZHEEV('N', p, ndim, vec_tmp, ndim, val, work, lwork, rwork,info)
endif
! check the status
if ( info /= 0 ) then
if ( present(ierr)) ierr = 3
print*, 'error in lapack subroutine zheev'
return
endif !
! dealloate memory for workspace array
if ( allocated(work) ) deallocate(work)
if ( allocated(rwork) ) deallocate(rwork)
if(.not. present(vec)) then
if ( allocated(vec_tmp ) ) deallocate(vec_tmp )
endif
end subroutine eigen_h_c
end module eigen_solver