(git:374b731)
Loading...
Searching...
No Matches
arnoldi_geev.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief provides a unified interface to lapack geev routines
10!> \par History
11!> 2014.09 created [Florian Schiffmann]
12!> \author Florian Schiffmann
13! **************************************************************************************************
14
15MODULE arnoldi_geev
16 USE kinds, ONLY: real_4,&
17 real_8
18#include "../base/base_uses.f90"
19
20 IMPLICIT NONE
21
22 PRIVATE
23
24 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_geev'
25
26 PUBLIC :: arnoldi_general_local_diag, arnoldi_tridiag_local_diag, arnoldi_symm_local_diag
27
28 INTERFACE arnoldi_general_local_diag
29 MODULE PROCEDURE arnoldi_sgeev, arnoldi_dgeev, arnoldi_zgeev, arnoldi_cgeev
30 END INTERFACE
31
32 ! currently only specialzed for real matrices
33 INTERFACE arnoldi_tridiag_local_diag
34 MODULE PROCEDURE arnoldi_sstev, arnoldi_dstev, arnoldi_zgeev, arnoldi_cgeev
35 END INTERFACE
36
37 ! currently only specialzed for real matrices
38 INTERFACE arnoldi_symm_local_diag
39 MODULE PROCEDURE arnoldi_dsyevd, arnoldi_ssyevd, arnoldi_cheevd, arnoldi_zheevd
40 END INTERFACE
41
42CONTAINS
43
44! **************************************************************************************************
45!> \brief ...
46!> \param jobvr ...
47!> \param matrix ...
48!> \param ndim ...
49!> \param evals ...
50!> \param revec ...
51! **************************************************************************************************
52 SUBROUTINE arnoldi_zheevd(jobvr, matrix, ndim, evals, revec)
53 CHARACTER(1) :: jobvr
54 COMPLEX(real_8), DIMENSION(:, :) :: matrix
55 INTEGER :: ndim
56 COMPLEX(real_8), DIMENSION(:) :: evals
57 COMPLEX(real_8), DIMENSION(:, :) :: revec
58
59 INTEGER :: i, info, liwork, lrwork, lwork, &
60 iwork(3 + 5*ndim)
61 COMPLEX(real_8) :: work(2*ndim + ndim**2), &
62 tmp_array(ndim, ndim)
63 REAL(real_8) :: rwork(1 + 5*ndim + 2*ndim**2)
64
65 tmp_array(:, :) = matrix(:, :)
66 lwork = 2*ndim + ndim**2
67 lrwork = 1 + 5*ndim + 2*ndim**2
68 liwork = 3 + 5*ndim
69
70 CALL zheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info)
71
72 DO i = 1, ndim
73 revec(:, i) = tmp_array(:, i)
74 END DO
75
76 END SUBROUTINE arnoldi_zheevd
77
78! **************************************************************************************************
79!> \brief ...
80!> \param jobvr ...
81!> \param matrix ...
82!> \param ndim ...
83!> \param evals ...
84!> \param revec ...
85! **************************************************************************************************
86 SUBROUTINE arnoldi_cheevd(jobvr, matrix, ndim, evals, revec)
87 CHARACTER(1) :: jobvr
88 COMPLEX(real_4), DIMENSION(:, :) :: matrix
89 INTEGER :: ndim
90 COMPLEX(real_4), DIMENSION(:) :: evals
91 COMPLEX(real_4), DIMENSION(:, :) :: revec
92
93 INTEGER :: i, info, liwork, lrwork, lwork, &
94 iwork(3 + 5*ndim)
95 COMPLEX(real_4) :: work(2*ndim + ndim**2), &
96 tmp_array(ndim, ndim)
97 REAL(real_4) :: rwork(1 + 5*ndim + 2*ndim**2)
98
99 tmp_array(:, :) = matrix(:, :)
100 lwork = 2*ndim + ndim**2
101 lrwork = 1 + 5*ndim + 2*ndim**2
102 liwork = 3 + 5*ndim
103
104 CALL cheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info)
105
106 DO i = 1, ndim
107 revec(:, i) = tmp_array(:, i)
108 END DO
109
110 END SUBROUTINE arnoldi_cheevd
111
112! **************************************************************************************************
113!> \brief ...
114!> \param jobvr ...
115!> \param matrix ...
116!> \param ndim ...
117!> \param evals ...
118!> \param revec ...
119! **************************************************************************************************
120 SUBROUTINE arnoldi_dsyevd(jobvr, matrix, ndim, evals, revec)
121 CHARACTER(1) :: jobvr
122 REAL(real_8), DIMENSION(:, :) :: matrix
123 INTEGER :: ndim
124 COMPLEX(real_8), DIMENSION(:) :: evals
125 COMPLEX(real_8), DIMENSION(:, :) :: revec
126
127 INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
128 REAL(real_8) :: tmp_array(ndim, ndim), &
129 work(1 + 6*ndim + 2*ndim**2)
130 REAL(real_8), DIMENSION(ndim) :: eval
131
132 lwork = 1 + 6*ndim + 2*ndim**2
133 liwork = 3 + 5*ndim
134
135 tmp_array(:, :) = matrix(:, :)
136 CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
137
138 DO i = 1, ndim
139 revec(:, i) = cmplx(tmp_array(:, i), real(0.0, real_8), real_8)
140 evals(i) = cmplx(eval(i), 0.0, real_8)
141 END DO
142
143 END SUBROUTINE arnoldi_dsyevd
144
145! **************************************************************************************************
146!> \brief ...
147!> \param jobvr ...
148!> \param matrix ...
149!> \param ndim ...
150!> \param evals ...
151!> \param revec ...
152! **************************************************************************************************
153 SUBROUTINE arnoldi_ssyevd(jobvr, matrix, ndim, evals, revec)
154 CHARACTER(1) :: jobvr
155 REAL(real_4), DIMENSION(:, :) :: matrix
156 INTEGER :: ndim
157 COMPLEX(real_4), DIMENSION(:) :: evals
158 COMPLEX(real_4), DIMENSION(:, :) :: revec
159
160 INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
161 REAL(real_4) :: tmp_array(ndim, ndim), &
162 work(1 + 6*ndim + 2*ndim**2)
163 REAL(real_4), DIMENSION(ndim) :: eval
164
165 mark_used(jobvr) !the argument has to be here for the template to work
166 lwork = 1 + 6*ndim + 2*ndim**2
167 liwork = 3 + 5*ndim
168
169 tmp_array(:, :) = matrix(:, :)
170 CALL ssyevd("V", "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
171
172 DO i = 1, ndim
173 revec(:, i) = cmplx(tmp_array(:, i), real(0.0, real_4), real_4)
174 evals(i) = cmplx(eval(i), 0.0, real_4)
175 END DO
176
177 END SUBROUTINE arnoldi_ssyevd
178
179! **************************************************************************************************
180!> \brief ...
181!> \param jobvl ...
182!> \param jobvr ...
183!> \param matrix ...
184!> \param ndim ...
185!> \param evals ...
186!> \param revec ...
187!> \param levec ...
188! **************************************************************************************************
189 SUBROUTINE arnoldi_sstev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
190 CHARACTER(1) :: jobvl, jobvr
191 REAL(real_4), DIMENSION(:, :) :: matrix
192 INTEGER :: ndim
193 COMPLEX(real_4), DIMENSION(:) :: evals
194 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec
195
196 INTEGER :: i, info
197 REAL(real_4) :: work(20*ndim)
198 REAL(real_4), DIMENSION(ndim) :: diag, offdiag
199 REAL(real_4), DIMENSION(ndim, ndim) :: evec_r
200
201 mark_used(jobvl) !the argument has to be here for the template to work
202
203 levec(1, 1) = cmplx(0.0, 0.0, real_4)
204 info = 0
205 diag(ndim) = matrix(ndim, ndim)
206 DO i = 1, ndim - 1
207 diag(i) = matrix(i, i)
208 offdiag(i) = matrix(i + 1, i)
209 END DO
210
211 CALL sstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
212
213 DO i = 1, ndim
214 revec(:, i) = cmplx(evec_r(:, i), real(0.0, real_4), real_4)
215 evals(i) = cmplx(diag(i), 0.0, real_4)
216 END DO
217 END SUBROUTINE arnoldi_sstev
218
219! **************************************************************************************************
220!> \brief ...
221!> \param jobvl ...
222!> \param jobvr ...
223!> \param matrix ...
224!> \param ndim ...
225!> \param evals ...
226!> \param revec ...
227!> \param levec ...
228! **************************************************************************************************
229 SUBROUTINE arnoldi_dstev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
230 CHARACTER(1) :: jobvl, jobvr
231 REAL(real_8), DIMENSION(:, :) :: matrix
232 INTEGER :: ndim
233 COMPLEX(real_8), DIMENSION(:) :: evals
234 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec
235
236 INTEGER :: i, info
237 REAL(real_8) :: work(20*ndim)
238 REAL(real_8), DIMENSION(ndim) :: diag, offdiag
239 REAL(real_8), DIMENSION(ndim, ndim) :: evec_r
240
241 mark_used(jobvl) !the argument has to be here for the template to work
242
243 levec(1, 1) = cmplx(0.0, 0.0, real_8)
244 info = 0
245 diag(ndim) = matrix(ndim, ndim)
246 DO i = 1, ndim - 1
247 diag(i) = matrix(i, i)
248 offdiag(i) = matrix(i + 1, i)
249
250 END DO
251
252 CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
253
254 DO i = 1, ndim
255 revec(:, i) = cmplx(evec_r(:, i), real(0.0, real_8), real_8)
256 evals(i) = cmplx(diag(i), 0.0, real_8)
257 END DO
258 END SUBROUTINE arnoldi_dstev
259! **************************************************************************************************
260!> \brief ...
261!> \param jobvl ...
262!> \param jobvr ...
263!> \param matrix ...
264!> \param ndim ...
265!> \param evals ...
266!> \param revec ...
267!> \param levec ...
268! **************************************************************************************************
269 SUBROUTINE arnoldi_sgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
270 CHARACTER(1) :: jobvl, jobvr
271 REAL(real_4), DIMENSION(:, :) :: matrix
272 INTEGER :: ndim
273 COMPLEX(real_4), DIMENSION(:) :: evals
274 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec
275
276 INTEGER :: i, info, lwork
277 LOGICAL :: selects(ndim)
278 REAL(real_4) :: norm, tmp_array(ndim, ndim), &
279 work(20*ndim)
280 REAL(real_4), DIMENSION(ndim) :: eval1, eval2
281 REAL(real_4), DIMENSION(ndim, ndim) :: evec_l, evec_r
282
283 mark_used(jobvr) !the argument has to be here for the template to work
284 mark_used(jobvl) !the argument has to be here for the template to work
285
286 eval1 = real(0.0, real_4); eval2 = real(0.0, real_4)
287 tmp_array(:, :) = matrix(:, :)
288 ! ask lapack how much space it would like in the work vector, don't ask me why
289 lwork = -1
290 CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
291
292 lwork = min(20*ndim, int(work(1)))
293 CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
294 CALL strevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info)
295
296 ! compose the eigenvectors, lapacks way of storing them is a pain
297 ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec
298 ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized
299 i = 1
300 DO WHILE (i .LE. ndim)
301 IF (abs(eval2(i)) .LT. epsilon(real(0.0, real_4))) THEN
302 evec_r(:, i) = evec_r(:, i)/sqrt(dot_product(evec_r(:, i), evec_r(:, i)))
303 revec(:, i) = cmplx(evec_r(:, i), real(0.0, real_4), real_4)
304 levec(:, i) = cmplx(evec_l(:, i), real(0.0, real_4), real_4)
305 i = i + 1
306 ELSE IF (eval2(i) .GT. epsilon(real(0.0, real_4))) THEN
307 norm = sqrt(sum(evec_r(:, i)**2.0_real_4) + sum(evec_r(:, i + 1)**2.0_real_4))
308 revec(:, i) = cmplx(evec_r(:, i), evec_r(:, i + 1), real_4)/norm
309 revec(:, i + 1) = cmplx(evec_r(:, i), -evec_r(:, i + 1), real_4)/norm
310 levec(:, i) = cmplx(evec_l(:, i), evec_l(:, i + 1), real_4)
311 levec(:, i + 1) = cmplx(evec_l(:, i), -evec_l(:, i + 1), real_4)
312 i = i + 2
313 ELSE
314 cpabort('something went wrong while sorting the EV in arnoldi_geev')
315 END IF
316 END DO
317
318 ! this is to keep the interface consistent with complex geev
319 DO i = 1, ndim
320 evals(i) = cmplx(eval1(i), eval2(i), real_4)
321 END DO
322
323 END SUBROUTINE arnoldi_sgeev
324
325! **************************************************************************************************
326!> \brief ...
327!> \param jobvl ...
328!> \param jobvr ...
329!> \param matrix ...
330!> \param ndim ...
331!> \param evals ...
332!> \param revec ...
333!> \param levec ...
334! **************************************************************************************************
335 SUBROUTINE arnoldi_dgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
336 CHARACTER(1) :: jobvl, jobvr
337 REAL(real_8), DIMENSION(:, :) :: matrix
338 INTEGER :: ndim
339 COMPLEX(real_8), DIMENSION(:) :: evals
340 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec
341
342 INTEGER :: i, info, lwork
343 LOGICAL :: selects(ndim)
344 REAL(real_8) :: norm, tmp_array(ndim, ndim), &
345 work(20*ndim)
346 REAL(real_8), DIMENSION(ndim) :: eval1, eval2
347 REAL(real_8), DIMENSION(ndim, ndim) :: evec_l, evec_r
348
349 mark_used(jobvr) !the argument has to be here for the template to work
350 mark_used(jobvl) !the argument has to be here for the template to work
351
352 eval1 = real(0.0, real_8); eval2 = real(0.0, real_8)
353 tmp_array(:, :) = matrix(:, :)
354 ! ask lapack how much space it would like in the work vector, don't ask me why
355 lwork = -1
356 CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
357
358 lwork = min(20*ndim, int(work(1)))
359 CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
360 CALL dtrevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info)
361
362 ! compose the eigenvectors, lapacks way of storing them is a pain
363 ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec
364 ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized
365 i = 1
366 DO WHILE (i .LE. ndim)
367 IF (abs(eval2(i)) .LT. epsilon(real(0.0, real_8))) THEN
368 evec_r(:, i) = evec_r(:, i)/sqrt(dot_product(evec_r(:, i), evec_r(:, i)))
369 revec(:, i) = cmplx(evec_r(:, i), real(0.0, real_8), real_8)
370 levec(:, i) = cmplx(evec_l(:, i), real(0.0, real_8), real_8)
371 i = i + 1
372 ELSE IF (eval2(i) .GT. epsilon(real(0.0, real_8))) THEN
373 norm = sqrt(sum(evec_r(:, i)**2.0_real_8) + sum(evec_r(:, i + 1)**2.0_real_8))
374 revec(:, i) = cmplx(evec_r(:, i), evec_r(:, i + 1), real_8)/norm
375 revec(:, i + 1) = cmplx(evec_r(:, i), -evec_r(:, i + 1), real_8)/norm
376 levec(:, i) = cmplx(evec_l(:, i), evec_l(:, i + 1), real_8)
377 levec(:, i + 1) = cmplx(evec_l(:, i), -evec_l(:, i + 1), real_8)
378 i = i + 2
379 ELSE
380 cpabort('something went wrong while sorting the EV in arnoldi_geev')
381 END IF
382 END DO
383
384 ! this is to keep the interface consistent with complex geev
385 DO i = 1, ndim
386 evals(i) = cmplx(eval1(i), eval2(i), real_8)
387 END DO
388
389 END SUBROUTINE arnoldi_dgeev
390
391! **************************************************************************************************
392!> \brief ...
393!> \param jobvl ...
394!> \param jobvr ...
395!> \param matrix ...
396!> \param ndim ...
397!> \param evals ...
398!> \param revec ...
399!> \param levec ...
400! **************************************************************************************************
401 SUBROUTINE arnoldi_zgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
402 CHARACTER(1) :: jobvl, jobvr
403 COMPLEX(real_8), DIMENSION(:, :) :: matrix
404 INTEGER :: ndim
405 COMPLEX(real_8), DIMENSION(:) :: evals
406 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec
407
408 INTEGER :: info, lwork
409 COMPLEX(real_8) :: work(20*ndim), tmp_array(ndim, ndim)
410 REAL(real_8) :: work2(2*ndim)
411
412 evals = cmplx(0.0, 0.0, real_8)
413 ! ask lapack how much space it would like in the work vector, don't ask me why
414 lwork = -1
415 CALL zgeev(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
416 lwork = min(20*ndim, int(work(1)))
417
418 tmp_array(:, :) = matrix(:, :)
419 CALL zgeev(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
420
421 END SUBROUTINE arnoldi_zgeev
422
423! **************************************************************************************************
424!> \brief ...
425!> \param jobvl ...
426!> \param jobvr ...
427!> \param matrix ...
428!> \param ndim ...
429!> \param evals ...
430!> \param revec ...
431!> \param levec ...
432! **************************************************************************************************
433 SUBROUTINE arnoldi_cgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec)
434 CHARACTER(1) :: jobvl, jobvr
435 COMPLEX(real_4), DIMENSION(:, :) :: matrix
436 INTEGER :: ndim
437 COMPLEX(real_4), DIMENSION(:) :: evals
438 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec
439
440 INTEGER :: info, lwork
441 COMPLEX(real_4) :: work(20*ndim), tmp_array(ndim, ndim)
442 REAL(real_4) :: work2(2*ndim)
443
444 evals = cmplx(0.0, 0.0, real_4)
445 ! ask lapack how much space it would like in the work vector, don't ask me why
446 lwork = -1
447 CALL cgeev(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
448 lwork = min(20*ndim, int(work(1)))
449
450 tmp_array(:, :) = matrix(:, :)
451 CALL cgeev(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
452
453 END SUBROUTINE arnoldi_cgeev
454
455END MODULE arnoldi_geev
arnoldi_geev::arnoldi_general_local_diag
Definition arnoldi_geev.F:28
arnoldi_geev::arnoldi_symm_local_diag
Definition arnoldi_geev.F:38
arnoldi_geev::arnoldi_tridiag_local_diag
Definition arnoldi_geev.F:33
arnoldi_geev
provides a unified interface to lapack geev routines
Definition arnoldi_geev.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::real_4
integer, parameter, public real_4
Definition kinds.F:40
kinds::real_8
integer, parameter, public real_8
Definition kinds.F:41