(git:34ef472)
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 
15 MODULE 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 
42 CONTAINS
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 
455 END MODULE arnoldi_geev
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