Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : !> 2023.12 Removed support for single-precision [Ole Schuett]
13 : !> 2024.12 Removed support for complex input matrices [Ole Schuett]
14 : !> \author Florian Schiffmann
15 : ! **************************************************************************************************
16 : MODULE arnoldi_geev
17 : #if defined (__HAS_IEEE_EXCEPTIONS)
18 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
19 : ieee_set_halting_mode, &
20 : IEEE_ALL
21 : #endif
22 : USE kinds, ONLY: dp
23 : #include "../base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_geev'
30 :
31 : PUBLIC :: arnoldi_general_local_diag, arnoldi_tridiag_local_diag, arnoldi_symm_local_diag
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief ...
37 : !> \param jobvr ...
38 : !> \param matrix ...
39 : !> \param ndim ...
40 : !> \param evals ...
41 : !> \param revec ...
42 : ! **************************************************************************************************
43 5063 : SUBROUTINE arnoldi_symm_local_diag(jobvr, matrix, ndim, evals, revec)
44 : CHARACTER(1) :: jobvr
45 : REAL(dp), DIMENSION(:, :) :: matrix
46 : INTEGER :: ndim
47 : COMPLEX(dp), DIMENSION(:) :: evals
48 : COMPLEX(dp), DIMENSION(:, :) :: revec
49 :
50 10126 : INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
51 10126 : REAL(dp) :: tmp_array(ndim, ndim), &
52 10126 : work(1 + 6*ndim + 2*ndim**2)
53 10126 : REAL(dp), DIMENSION(ndim) :: eval
54 :
55 5063 : lwork = 1 + 6*ndim + 2*ndim**2
56 5063 : liwork = 3 + 5*ndim
57 :
58 1526179 : tmp_array(:, :) = matrix(:, :)
59 5063 : CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
60 :
61 87300 : DO i = 1, ndim
62 1521116 : revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, dp), dp)
63 87300 : evals(i) = CMPLX(eval(i), 0.0, dp)
64 : END DO
65 :
66 5063 : END SUBROUTINE arnoldi_symm_local_diag
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param jobvl ...
71 : !> \param jobvr ...
72 : !> \param matrix ...
73 : !> \param ndim ...
74 : !> \param evals ...
75 : !> \param revec ...
76 : !> \param levec ...
77 : ! **************************************************************************************************
78 9487 : SUBROUTINE arnoldi_tridiag_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
79 : CHARACTER(1) :: jobvl, jobvr
80 : REAL(dp), DIMENSION(:, :) :: matrix
81 : INTEGER :: ndim
82 : COMPLEX(dp), DIMENSION(:) :: evals
83 : COMPLEX(dp), DIMENSION(:, :) :: revec, levec
84 : #if defined (__HAS_IEEE_EXCEPTIONS)
85 : LOGICAL, DIMENSION(5) :: halt
86 : #endif
87 : INTEGER :: i, info
88 18974 : REAL(dp) :: work(20*ndim)
89 18974 : REAL(dp), DIMENSION(ndim) :: diag, offdiag
90 18974 : REAL(dp), DIMENSION(ndim, ndim) :: evec_r
91 :
92 : MARK_USED(jobvl) !the argument has to be here for the template to work
93 :
94 9487 : levec(1, 1) = CMPLX(0.0, 0.0, dp)
95 9487 : info = 0
96 9487 : diag(ndim) = matrix(ndim, ndim)
97 98427 : DO i = 1, ndim - 1
98 88940 : diag(i) = matrix(i, i)
99 98427 : offdiag(i) = matrix(i + 1, i)
100 :
101 : END DO
102 :
103 : #if defined (__HAS_IEEE_EXCEPTIONS)
104 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
105 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
106 : #endif
107 :
108 9487 : CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
109 :
110 : #if defined (__HAS_IEEE_EXCEPTIONS)
111 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
112 : #endif
113 :
114 9487 : CPASSERT(info == 0)
115 :
116 107914 : DO i = 1, ndim
117 2251278 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
118 107914 : evals(i) = CMPLX(diag(i), 0.0, dp)
119 : END DO
120 9487 : END SUBROUTINE arnoldi_tridiag_local_diag
121 :
122 : ! **************************************************************************************************
123 : !> \brief ...
124 : !> \param jobvl ...
125 : !> \param jobvr ...
126 : !> \param matrix ...
127 : !> \param ndim ...
128 : !> \param evals ...
129 : !> \param revec ...
130 : !> \param levec ...
131 : ! **************************************************************************************************
132 134489 : SUBROUTINE arnoldi_general_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
133 : CHARACTER(1) :: jobvl, jobvr
134 : REAL(dp), DIMENSION(:, :) :: matrix
135 : INTEGER :: ndim
136 : COMPLEX(dp), DIMENSION(:) :: evals
137 : COMPLEX(dp), DIMENSION(:, :) :: revec, levec
138 :
139 : INTEGER :: i, info, lwork
140 268978 : LOGICAL :: selects(ndim)
141 268978 : REAL(dp) :: norm, tmp_array(ndim, ndim), &
142 268978 : work(20*ndim)
143 268978 : REAL(dp), DIMENSION(ndim) :: eval1, eval2
144 134489 : REAL(dp), DIMENSION(ndim, ndim) :: evec_l, evec_r
145 :
146 : MARK_USED(jobvr) !the argument has to be here for the template to work
147 : MARK_USED(jobvl) !the argument has to be here for the template to work
148 :
149 1268987 : eval1 = REAL(0.0, dp); eval2 = REAL(0.0, dp)
150 6924021 : tmp_array(:, :) = matrix(:, :)
151 : ! ask lapack how much space it would like in the work vector, don't ask me why
152 134489 : lwork = -1
153 134489 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
154 :
155 134489 : lwork = MIN(20*ndim, INT(work(1)))
156 134489 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
157 134489 : CALL dtrevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info)
158 :
159 : ! compose the eigenvectors, lapacks way of storing them is a pain
160 : ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec
161 : ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized
162 134489 : i = 1
163 701736 : DO WHILE (i <= ndim)
164 701736 : IF (ABS(eval2(i)) < EPSILON(REAL(0.0, dp))) THEN
165 13011739 : evec_r(:, i) = evec_r(:, i)/NORM2(evec_r(:, i))
166 6789492 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
167 6789492 : levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, dp), dp)
168 567245 : i = i + 1
169 2 : ELSE IF (eval2(i) > EPSILON(REAL(0.0, dp))) THEN
170 38 : norm = SQRT(SUM(evec_r(:, i)**2.0_dp) + SUM(evec_r(:, i + 1)**2.0_dp))
171 20 : revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), dp)/norm
172 20 : revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), dp)/norm
173 20 : levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), dp)
174 20 : levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), dp)
175 2 : i = i + 2
176 : ELSE
177 0 : CPABORT('something went wrong while sorting the EV in arnoldi_geev')
178 : END IF
179 : END DO
180 :
181 : ! this is to keep the interface consistent with complex geev
182 701738 : DO i = 1, ndim
183 701738 : evals(i) = CMPLX(eval1(i), eval2(i), dp)
184 : END DO
185 :
186 134489 : END SUBROUTINE arnoldi_general_local_diag
187 :
188 : END MODULE arnoldi_geev
|
