Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : USE kinds, ONLY: dp
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 : CONTAINS
29 :
30 : ! **************************************************************************************************
31 : !> \brief ...
32 : !> \param jobvr ...
33 : !> \param matrix ...
34 : !> \param ndim ...
35 : !> \param evals ...
36 : !> \param revec ...
37 : ! **************************************************************************************************
38 4503 : SUBROUTINE arnoldi_symm_local_diag(jobvr, matrix, ndim, evals, revec)
39 : CHARACTER(1) :: jobvr
40 : REAL(dp), DIMENSION(:, :) :: matrix
41 : INTEGER :: ndim
42 : COMPLEX(dp), DIMENSION(:) :: evals
43 : COMPLEX(dp), DIMENSION(:, :) :: revec
44 :
45 9006 : INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim)
46 9006 : REAL(dp) :: tmp_array(ndim, ndim), &
47 9006 : work(1 + 6*ndim + 2*ndim**2)
48 9006 : REAL(dp), DIMENSION(ndim) :: eval
49 :
50 4503 : lwork = 1 + 6*ndim + 2*ndim**2
51 4503 : liwork = 3 + 5*ndim
52 :
53 1392311 : tmp_array(:, :) = matrix(:, :)
54 4503 : CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
55 :
56 79276 : DO i = 1, ndim
57 1387808 : revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, dp), dp)
58 79276 : evals(i) = CMPLX(eval(i), 0.0, dp)
59 : END DO
60 :
61 4503 : END SUBROUTINE arnoldi_symm_local_diag
62 :
63 : ! **************************************************************************************************
64 : !> \brief ...
65 : !> \param jobvl ...
66 : !> \param jobvr ...
67 : !> \param matrix ...
68 : !> \param ndim ...
69 : !> \param evals ...
70 : !> \param revec ...
71 : !> \param levec ...
72 : ! **************************************************************************************************
73 8701 : SUBROUTINE arnoldi_tridiag_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
74 : CHARACTER(1) :: jobvl, jobvr
75 : REAL(dp), DIMENSION(:, :) :: matrix
76 : INTEGER :: ndim
77 : COMPLEX(dp), DIMENSION(:) :: evals
78 : COMPLEX(dp), DIMENSION(:, :) :: revec, levec
79 :
80 : INTEGER :: i, info
81 17402 : REAL(dp) :: work(20*ndim)
82 17402 : REAL(dp), DIMENSION(ndim) :: diag, offdiag
83 17402 : REAL(dp), DIMENSION(ndim, ndim) :: evec_r
84 :
85 : MARK_USED(jobvl) !the argument has to be here for the template to work
86 :
87 8701 : levec(1, 1) = CMPLX(0.0, 0.0, dp)
88 8701 : info = 0
89 8701 : diag(ndim) = matrix(ndim, ndim)
90 96203 : DO i = 1, ndim - 1
91 87502 : diag(i) = matrix(i, i)
92 96203 : offdiag(i) = matrix(i + 1, i)
93 :
94 : END DO
95 :
96 8701 : CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
97 :
98 104904 : DO i = 1, ndim
99 2231494 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
100 104904 : evals(i) = CMPLX(diag(i), 0.0, dp)
101 : END DO
102 8701 : END SUBROUTINE arnoldi_tridiag_local_diag
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param jobvl ...
107 : !> \param jobvr ...
108 : !> \param matrix ...
109 : !> \param ndim ...
110 : !> \param evals ...
111 : !> \param revec ...
112 : !> \param levec ...
113 : ! **************************************************************************************************
114 123025 : SUBROUTINE arnoldi_general_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
115 : CHARACTER(1) :: jobvl, jobvr
116 : REAL(dp), DIMENSION(:, :) :: matrix
117 : INTEGER :: ndim
118 : COMPLEX(dp), DIMENSION(:) :: evals
119 : COMPLEX(dp), DIMENSION(:, :) :: revec, levec
120 :
121 : INTEGER :: i, info, lwork
122 246050 : LOGICAL :: selects(ndim)
123 246050 : REAL(dp) :: norm, tmp_array(ndim, ndim), &
124 246050 : work(20*ndim)
125 246050 : REAL(dp), DIMENSION(ndim) :: eval1, eval2
126 123025 : REAL(dp), DIMENSION(ndim, ndim) :: evec_l, evec_r
127 :
128 : MARK_USED(jobvr) !the argument has to be here for the template to work
129 : MARK_USED(jobvl) !the argument has to be here for the template to work
130 :
131 1192555 : eval1 = REAL(0.0, dp); eval2 = REAL(0.0, dp)
132 6688773 : tmp_array(:, :) = matrix(:, :)
133 : ! ask lapack how much space it would like in the work vector, don't ask me why
134 123025 : lwork = -1
135 123025 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
136 :
137 123025 : lwork = MIN(20*ndim, INT(work(1)))
138 123025 : CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
139 123025 : CALL dtrevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info)
140 :
141 : ! compose the eigenvectors, lapacks way of storing them is a pain
142 : ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec
143 : ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized
144 123025 : i = 1
145 657790 : DO WHILE (i .LE. ndim)
146 657790 : IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, dp))) THEN
147 12596731 : evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i)))
148 6565748 : revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
149 6565748 : levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, dp), dp)
150 534765 : i = i + 1
151 0 : ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, dp))) THEN
152 0 : norm = SQRT(SUM(evec_r(:, i)**2.0_dp) + SUM(evec_r(:, i + 1)**2.0_dp))
153 0 : revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), dp)/norm
154 0 : revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), dp)/norm
155 0 : levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), dp)
156 0 : levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), dp)
157 0 : i = i + 2
158 : ELSE
159 0 : CPABORT('something went wrong while sorting the EV in arnoldi_geev')
160 : END IF
161 : END DO
162 :
163 : ! this is to keep the interface consistent with complex geev
164 657790 : DO i = 1, ndim
165 657790 : evals(i) = CMPLX(eval1(i), eval2(i), dp)
166 : END DO
167 :
168 123025 : END SUBROUTINE arnoldi_general_local_diag
169 :
170 : END MODULE arnoldi_geev
|
