Line data Source code
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 0 : 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 0 : iwork(3 + 5*ndim)
61 0 : COMPLEX(real_8) :: work(2*ndim + ndim**2), &
62 0 : tmp_array(ndim, ndim)
63 0 : REAL(real_8) :: rwork(1 + 5*ndim + 2*ndim**2)
64 :
65 0 : tmp_array(:, :) = matrix(:, :)
66 0 : lwork = 2*ndim + ndim**2
67 0 : lrwork = 1 + 5*ndim + 2*ndim**2
68 0 : liwork = 3 + 5*ndim
69 :
70 0 : CALL zheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info)
71 :
72 0 : DO i = 1, ndim
73 0 : revec(:, i) = tmp_array(:, i)
74 : END DO
75 :
76 0 : 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 0 : 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 0 : iwork(3 + 5*ndim)
95 0 : COMPLEX(real_4) :: work(2*ndim + ndim**2), &
96 0 : tmp_array(ndim, ndim)
97 0 : REAL(real_4) :: rwork(1 + 5*ndim + 2*ndim**2)
98 :
99 0 : tmp_array(:, :) = matrix(:, :)
100 0 : lwork = 2*ndim + ndim**2
101 0 : lrwork = 1 + 5*ndim + 2*ndim**2
102 0 : liwork = 3 + 5*ndim
103 :
104 0 : CALL cheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info)
105 :
106 0 : DO i = 1, ndim
107 0 : revec(:, i) = tmp_array(:, i)
108 : END DO
109 :
110 0 : 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 4605 : 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 9210 : INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
128 9210 : REAL(real_8) :: tmp_array(ndim, ndim), &
129 9210 : work(1 + 6*ndim + 2*ndim**2)
130 9210 : REAL(real_8), DIMENSION(ndim) :: eval
131 :
132 4605 : lwork = 1 + 6*ndim + 2*ndim**2
133 4605 : liwork = 3 + 5*ndim
134 :
135 1429005 : tmp_array(:, :) = matrix(:, :)
136 4605 : CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
137 :
138 81228 : DO i = 1, ndim
139 1424400 : revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, real_8), real_8)
140 81228 : evals(i) = CMPLX(eval(i), 0.0, real_8)
141 : END DO
142 :
143 4605 : 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 0 : 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 0 : INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
161 0 : REAL(real_4) :: tmp_array(ndim, ndim), &
162 0 : work(1 + 6*ndim + 2*ndim**2)
163 0 : REAL(real_4), DIMENSION(ndim) :: eval
164 :
165 : MARK_USED(jobvr) !the argument has to be here for the template to work
166 0 : lwork = 1 + 6*ndim + 2*ndim**2
167 0 : liwork = 3 + 5*ndim
168 :
169 0 : tmp_array(:, :) = matrix(:, :)
170 0 : CALL ssyevd("V", "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
171 :
172 0 : DO i = 1, ndim
173 0 : revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, real_4), real_4)
174 0 : evals(i) = CMPLX(eval(i), 0.0, real_4)
175 : END DO
176 :
177 0 : 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 0 : 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 0 : REAL(real_4) :: work(20*ndim)
198 0 : REAL(real_4), DIMENSION(ndim) :: diag, offdiag
199 0 : 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 0 : levec(1, 1) = CMPLX(0.0, 0.0, real_4)
204 0 : info = 0
205 0 : diag(ndim) = matrix(ndim, ndim)
206 0 : DO i = 1, ndim - 1
207 0 : diag(i) = matrix(i, i)
208 0 : offdiag(i) = matrix(i + 1, i)
209 : END DO
210 :
211 0 : CALL sstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
212 :
213 0 : DO i = 1, ndim
214 0 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_4), real_4)
215 0 : evals(i) = CMPLX(diag(i), 0.0, real_4)
216 : END DO
217 0 : 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 8650 : 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 17300 : REAL(real_8) :: work(20*ndim)
238 17300 : REAL(real_8), DIMENSION(ndim) :: diag, offdiag
239 17300 : 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 8650 : levec(1, 1) = CMPLX(0.0, 0.0, real_8)
244 8650 : info = 0
245 8650 : diag(ndim) = matrix(ndim, ndim)
246 95373 : DO i = 1, ndim - 1
247 86723 : diag(i) = matrix(i, i)
248 95373 : offdiag(i) = matrix(i + 1, i)
249 :
250 : END DO
251 :
252 8650 : CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
253 :
254 104023 : DO i = 1, ndim
255 2214314 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_8), real_8)
256 104023 : evals(i) = CMPLX(diag(i), 0.0, real_8)
257 : END DO
258 8650 : 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 0 : 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 0 : LOGICAL :: selects(ndim)
278 0 : REAL(real_4) :: norm, tmp_array(ndim, ndim), &
279 0 : work(20*ndim)
280 0 : REAL(real_4), DIMENSION(ndim) :: eval1, eval2
281 0 : 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 0 : eval1 = REAL(0.0, real_4); eval2 = REAL(0.0, real_4)
287 0 : tmp_array(:, :) = matrix(:, :)
288 : ! ask lapack how much space it would like in the work vector, don't ask me why
289 0 : lwork = -1
290 0 : CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
291 :
292 0 : lwork = MIN(20*ndim, INT(work(1)))
293 0 : CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
294 0 : 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 0 : i = 1
300 0 : DO WHILE (i .LE. ndim)
301 0 : IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, real_4))) THEN
302 0 : evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i)))
303 0 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_4), real_4)
304 0 : levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, real_4), real_4)
305 0 : i = i + 1
306 0 : ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, real_4))) THEN
307 0 : norm = SQRT(SUM(evec_r(:, i)**2.0_real_4) + SUM(evec_r(:, i + 1)**2.0_real_4))
308 0 : revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), real_4)/norm
309 0 : revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), real_4)/norm
310 0 : levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), real_4)
311 0 : levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), real_4)
312 0 : i = i + 2
313 : ELSE
314 0 : 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 0 : DO i = 1, ndim
320 0 : evals(i) = CMPLX(eval1(i), eval2(i), real_4)
321 : END DO
322 :
323 0 : 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 116659 : 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 233318 : LOGICAL :: selects(ndim)
344 233318 : REAL(real_8) :: norm, tmp_array(ndim, ndim), &
345 233318 : work(20*ndim)
346 233318 : REAL(real_8), DIMENSION(ndim) :: eval1, eval2
347 116659 : 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 1112429 : eval1 = REAL(0.0, real_8); eval2 = REAL(0.0, real_8)
353 6098613 : tmp_array(:, :) = matrix(:, :)
354 : ! ask lapack how much space it would like in the work vector, don't ask me why
355 116659 : lwork = -1
356 116659 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
357 :
358 116659 : lwork = MIN(20*ndim, INT(work(1)))
359 116659 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
360 116659 : 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 116659 : i = 1
366 614544 : DO WHILE (i .LE. ndim)
367 614544 : IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, real_8))) THEN
368 11466023 : evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i)))
369 5981954 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_8), real_8)
370 5981954 : levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, real_8), real_8)
371 497885 : i = i + 1
372 0 : ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, real_8))) THEN
373 0 : norm = SQRT(SUM(evec_r(:, i)**2.0_real_8) + SUM(evec_r(:, i + 1)**2.0_real_8))
374 0 : revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), real_8)/norm
375 0 : revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), real_8)/norm
376 0 : levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), real_8)
377 0 : levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), real_8)
378 0 : i = i + 2
379 : ELSE
380 0 : 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 614544 : DO i = 1, ndim
386 614544 : evals(i) = CMPLX(eval1(i), eval2(i), real_8)
387 : END DO
388 :
389 116659 : 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 0 : 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 0 : COMPLEX(real_8) :: work(20*ndim), tmp_array(ndim, ndim)
410 0 : REAL(real_8) :: work2(2*ndim)
411 :
412 0 : 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 0 : lwork = -1
415 0 : CALL ZGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
416 0 : lwork = MIN(20*ndim, INT(work(1)))
417 :
418 0 : tmp_array(:, :) = matrix(:, :)
419 0 : CALL ZGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
420 :
421 0 : 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 0 : 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 0 : COMPLEX(real_4) :: work(20*ndim), tmp_array(ndim, ndim)
442 0 : REAL(real_4) :: work2(2*ndim)
443 :
444 0 : 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 0 : lwork = -1
447 0 : CALL CGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
448 0 : lwork = MIN(20*ndim, INT(work(1)))
449 :
450 0 : tmp_array(:, :) = matrix(:, :)
451 0 : CALL CGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info)
452 :
453 0 : END SUBROUTINE arnoldi_cgeev
454 :
455 : END MODULE arnoldi_geev
|
