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 used for collecting diagonalization schemes available for cp_cfm_type
10 : !> \note
11 : !> first version : only one routine right now
12 : !> \author Joost VandeVondele (2003-09)
13 : ! **************************************************************************************************
14 : MODULE cp_cfm_diag
15 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose, &
16 : cp_cfm_gemm, &
17 : cp_cfm_column_scale, &
18 : cp_cfm_scale, &
19 : cp_cfm_triangular_invert, &
20 : cp_cfm_triangular_multiply
21 : USE cp_cfm_types, ONLY: cp_cfm_get_info, &
22 : cp_cfm_set_element, &
23 : cp_cfm_to_cfm, &
24 : cp_cfm_type
25 : USE kinds, ONLY: dp
26 : #if defined (__HAS_IEEE_EXCEPTIONS)
27 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
28 : ieee_set_halting_mode, &
29 : ieee_all
30 : #endif
31 :
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 : PRIVATE
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
37 :
38 : PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief Perform a diagonalisation of a complex matrix
44 : !> \param matrix ...
45 : !> \param eigenvectors ...
46 : !> \param eigenvalues ...
47 : !> \par History
48 : !> - (De)Allocation checks updated (15.02.2011,MK)
49 : !> \author Joost VandeVondele
50 : ! **************************************************************************************************
51 25518 : SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
52 :
53 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
54 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
55 :
56 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
57 :
58 25518 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
59 : COMPLEX(KIND=dp), DIMENSION(:, :), &
60 25518 : POINTER :: m
61 : INTEGER :: handle, info, liwork, &
62 : lrwork, lwork, n
63 25518 : INTEGER, DIMENSION(:), POINTER :: iwork
64 25518 : REAL(KIND=dp), DIMENSION(:), POINTER :: rwork
65 : #if defined(__SCALAPACK)
66 : INTEGER, DIMENSION(9) :: descm, descv
67 : COMPLEX(KIND=dp), DIMENSION(:, :), &
68 25518 : POINTER :: v
69 : #if defined (__HAS_IEEE_EXCEPTIONS)
70 : LOGICAL, DIMENSION(5) :: halt
71 : #endif
72 : #endif
73 :
74 25518 : CALL timeset(routineN, handle)
75 :
76 25518 : n = matrix%matrix_struct%nrow_global
77 25518 : m => matrix%local_data
78 25518 : ALLOCATE (iwork(1), rwork(1), work(1))
79 : ! work space query
80 25518 : lwork = -1
81 25518 : lrwork = -1
82 25518 : liwork = -1
83 :
84 : #if defined(__SCALAPACK)
85 25518 : v => eigenvectors%local_data
86 255180 : descm(:) = matrix%matrix_struct%descriptor(:)
87 255180 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
88 : CALL PZHEEVD('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
89 25518 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
90 : ! The work space query for lwork does not return always sufficiently large values.
91 : ! Let's add some margin to avoid crashes.
92 25518 : lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
93 : ! needed to correct for a bug in scalapack, unclear how much the right number is
94 25518 : lrwork = CEILING(REAL(rwork(1), KIND=dp)) + 1000000
95 25518 : liwork = iwork(1)
96 : #else
97 : CALL ZHEEVD('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
98 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
99 : lwork = CEILING(REAL(work(1), KIND=dp))
100 : lrwork = CEILING(REAL(rwork(1), KIND=dp))
101 : liwork = iwork(1)
102 : #endif
103 :
104 25518 : DEALLOCATE (iwork, rwork, work)
105 178626 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
106 :
107 : #if defined(__SCALAPACK)
108 : ! Scalapack takes advantage of IEEE754 exceptions for speedup.
109 : ! Therefore, we disable floating point traps temporarily.
110 : #if defined (__HAS_IEEE_EXCEPTIONS)
111 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
112 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
113 : #endif
114 :
115 : CALL PZHEEVD('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
116 25518 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
117 :
118 : #if defined (__HAS_IEEE_EXCEPTIONS)
119 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
120 : #endif
121 : #else
122 : CALL ZHEEVD('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
123 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
124 : eigenvectors%local_data = matrix%local_data
125 : #endif
126 :
127 25518 : DEALLOCATE (iwork, rwork, work)
128 25518 : IF (info /= 0) &
129 0 : CPABORT("Diagonalisation of a complex matrix failed")
130 :
131 25518 : CALL timestop(handle)
132 :
133 25518 : END SUBROUTINE cp_cfm_heevd
134 :
135 : ! **************************************************************************************************
136 : !> \brief General Eigenvalue Problem AX = BXE
137 : !> Single option version: Cholesky decomposition of B
138 : !> \param amatrix ...
139 : !> \param bmatrix ...
140 : !> \param eigenvectors ...
141 : !> \param eigenvalues ...
142 : !> \param work ...
143 : ! **************************************************************************************************
144 15128 : SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
145 :
146 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
147 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
148 : TYPE(cp_cfm_type), INTENT(IN) :: work
149 :
150 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig'
151 :
152 : INTEGER :: handle, nao, nmo
153 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
154 :
155 15128 : CALL timeset(routineN, handle)
156 :
157 15128 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
158 45384 : ALLOCATE (evals(nao))
159 : ! Cholesky decompose S=U(T)U
160 15128 : CALL cp_cfm_cholesky_decompose(bmatrix)
161 : ! Invert to get U^(-1)
162 15128 : CALL cp_cfm_triangular_invert(bmatrix)
163 : ! Reduce to get U^(-T) * H * U^(-1)
164 15128 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
165 15128 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
166 : ! Diagonalize
167 15128 : CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
168 : ! Restore vectors C = U^(-1) * C*
169 15128 : CALL cp_cfm_triangular_multiply(bmatrix, work)
170 15128 : nmo = SIZE(eigenvalues)
171 15128 : CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
172 236198 : eigenvalues(1:nmo) = evals(1:nmo)
173 :
174 15128 : DEALLOCATE (evals)
175 :
176 15128 : CALL timestop(handle)
177 :
178 15128 : END SUBROUTINE cp_cfm_geeig
179 :
180 : ! **************************************************************************************************
181 : !> \brief General Eigenvalue Problem AX = BXE
182 : !> Use canonical orthogonalization
183 : !> \param amatrix ...
184 : !> \param bmatrix ...
185 : !> \param eigenvectors ...
186 : !> \param eigenvalues ...
187 : !> \param work ...
188 : !> \param epseig ...
189 : ! **************************************************************************************************
190 274 : SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
191 :
192 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
193 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
194 : TYPE(cp_cfm_type), INTENT(IN) :: work
195 : REAL(KIND=dp), INTENT(IN) :: epseig
196 :
197 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
198 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
199 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
200 :
201 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
202 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
203 : nmo, nx
204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
205 :
206 274 : CALL timeset(routineN, handle)
207 :
208 : ! Test sizes
209 274 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
210 274 : nmo = SIZE(eigenvalues)
211 1370 : ALLOCATE (evals(nao), cevals(nao))
212 :
213 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
214 274 : CALL cp_cfm_scale(-cone, bmatrix)
215 274 : CALL cp_cfm_heevd(bmatrix, work, evals)
216 17602 : evals(:) = -evals(:)
217 274 : nc = nao
218 17038 : DO i = 1, nao
219 17038 : IF (evals(i) < epseig) THEN
220 72 : nc = i - 1
221 72 : EXIT
222 : END IF
223 : END DO
224 274 : CPASSERT(nc /= 0)
225 :
226 274 : IF (nc /= nao) THEN
227 72 : IF (nc < nmo) THEN
228 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
229 0 : ncol = nmo - nc
230 0 : CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
231 : END IF
232 : ! Set NULL space in eigenvector matrix of S to zero
233 636 : DO icol = nc + 1, nao
234 80028 : DO irow = 1, nao
235 79956 : CALL cp_cfm_set_element(work, irow, icol, czero)
236 : END DO
237 : END DO
238 : ! Set small eigenvalues to a dummy save value
239 636 : evals(nc + 1:nao) = 1.0_dp
240 : END IF
241 : ! calculate U*s**(-1/2)
242 17602 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
243 274 : CALL cp_cfm_column_scale(work, cevals)
244 : ! Reduce to get U^(-C) * H * U^(-1)
245 274 : CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
246 274 : CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
247 274 : IF (nc /= nao) THEN
248 : ! set diagonal values to save large value
249 636 : DO icol = nc + 1, nao
250 636 : CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp))
251 : END DO
252 : END IF
253 : ! Diagonalize
254 274 : CALL cp_cfm_heevd(amatrix, bmatrix, evals)
255 4778 : eigenvalues(1:nmo) = evals(1:nmo)
256 274 : nx = MIN(nc, nmo)
257 : ! Restore vectors C = U^(-1) * C*
258 274 : CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
259 :
260 274 : DEALLOCATE (evals)
261 :
262 274 : CALL timestop(handle)
263 :
264 548 : END SUBROUTINE cp_cfm_geeig_canon
265 :
266 : END MODULE cp_cfm_diag
|
