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 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_cholesky, ONLY: cp_cfm_cholesky_decompose
16 : USE cp_cfm_basic_linalg, ONLY: 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 : #if defined(__DLAF)
26 : USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
27 : cp_cfm_diag_dlaf
28 : USE cp_fm_diag, ONLY: diag_type, dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
29 : #endif
30 : USE kinds, ONLY: dp
31 : #if defined (__HAS_IEEE_EXCEPTIONS)
32 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
33 : ieee_set_halting_mode, &
34 : IEEE_ALL
35 : #endif
36 :
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
43 :
44 : PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Perform a diagonalisation of a complex matrix
50 : !> \param matrix ...
51 : !> \param eigenvectors ...
52 : !> \param eigenvalues ...
53 : !> \par History
54 : !> 12.2024 Added DLA-Future support [Rocco Meli]
55 : !> \author Joost VandeVondele
56 : ! **************************************************************************************************
57 25662 : SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
58 :
59 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
60 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
61 :
62 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
63 :
64 : INTEGER :: handle
65 :
66 25662 : CALL timeset(routineN, handle)
67 :
68 : #if defined(__DLAF)
69 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
70 : CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
71 : ELSE
72 : #endif
73 25662 : CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
74 : #if defined(__DLAF)
75 : END IF
76 : #endif
77 :
78 25662 : CALL timestop(handle)
79 :
80 25662 : END SUBROUTINE cp_cfm_heevd
81 :
82 : ! **************************************************************************************************
83 : !> \brief Perform a diagonalisation of a complex matrix
84 : !> \param matrix ...
85 : !> \param eigenvectors ...
86 : !> \param eigenvalues ...
87 : !> \par History
88 : !> - (De)Allocation checks updated (15.02.2011,MK)
89 : !> \author Joost VandeVondele
90 : ! **************************************************************************************************
91 25662 : SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
92 :
93 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
94 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
95 :
96 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
97 :
98 25662 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
99 : COMPLEX(KIND=dp), DIMENSION(:, :), &
100 25662 : POINTER :: m
101 : INTEGER :: handle, info, liwork, &
102 : lrwork, lwork, n
103 25662 : INTEGER, DIMENSION(:), POINTER :: iwork
104 25662 : REAL(KIND=dp), DIMENSION(:), POINTER :: rwork
105 : #if defined(__parallel)
106 : INTEGER, DIMENSION(9) :: descm, descv
107 : COMPLEX(KIND=dp), DIMENSION(:, :), &
108 25662 : POINTER :: v
109 : #if defined (__HAS_IEEE_EXCEPTIONS)
110 : LOGICAL, DIMENSION(5) :: halt
111 : #endif
112 : #endif
113 :
114 25662 : CALL timeset(routineN, handle)
115 :
116 25662 : n = matrix%matrix_struct%nrow_global
117 25662 : m => matrix%local_data
118 25662 : ALLOCATE (iwork(1), rwork(1), work(1))
119 : ! work space query
120 25662 : lwork = -1
121 25662 : lrwork = -1
122 25662 : liwork = -1
123 :
124 : #if defined(__parallel)
125 25662 : v => eigenvectors%local_data
126 256620 : descm(:) = matrix%matrix_struct%descriptor(:)
127 256620 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
128 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
129 25662 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
130 : ! The work space query for lwork does not return always sufficiently large values.
131 : ! Let's add some margin to avoid crashes.
132 25662 : lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
133 : ! needed to correct for a bug in scalapack, unclear how much the right number is
134 25662 : lrwork = CEILING(rwork(1)) + 1000000
135 25662 : liwork = iwork(1)
136 : #else
137 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
138 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
139 : lwork = CEILING(REAL(work(1), KIND=dp))
140 : lrwork = CEILING(rwork(1))
141 : liwork = iwork(1)
142 : #endif
143 :
144 25662 : DEALLOCATE (iwork, rwork, work)
145 179634 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
146 :
147 : #if defined(__parallel)
148 : ! Scalapack takes advantage of IEEE754 exceptions for speedup.
149 : ! Therefore, we disable floating point traps temporarily.
150 : #if defined (__HAS_IEEE_EXCEPTIONS)
151 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
152 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
153 : #endif
154 :
155 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
156 25662 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
157 :
158 : #if defined (__HAS_IEEE_EXCEPTIONS)
159 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
160 : #endif
161 : #else
162 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
163 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
164 : eigenvectors%local_data = matrix%local_data
165 : #endif
166 :
167 25662 : DEALLOCATE (iwork, rwork, work)
168 25662 : IF (info /= 0) &
169 0 : CPABORT("Diagonalisation of a complex matrix failed")
170 :
171 25662 : CALL timestop(handle)
172 :
173 25662 : END SUBROUTINE cp_cfm_heevd_base
174 :
175 : ! **************************************************************************************************
176 : !> \brief General Eigenvalue Problem AX = BXE
177 : !> Single option version: Cholesky decomposition of B
178 : !> \param amatrix ...
179 : !> \param bmatrix ...
180 : !> \param eigenvectors ...
181 : !> \param eigenvalues ...
182 : !> \param work ...
183 : !> \par History
184 : !> 12.2024 Added DLA-Future support [Rocco Meli]
185 : ! **************************************************************************************************
186 15568 : SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
187 :
188 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
189 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
190 : TYPE(cp_cfm_type), INTENT(IN) :: work
191 :
192 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig'
193 :
194 : INTEGER :: handle, nao, nmo
195 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
196 :
197 15568 : CALL timeset(routineN, handle)
198 :
199 15568 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
200 46704 : ALLOCATE (evals(nao))
201 15568 : nmo = SIZE(eigenvalues)
202 :
203 : #if defined(__DLAF)
204 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
205 : ! Use DLA-Future generalized eigenvalue solver for large matrices
206 : CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
207 : ELSE
208 : #endif
209 : ! Cholesky decompose S=U(T)U
210 15568 : CALL cp_cfm_cholesky_decompose(bmatrix)
211 : ! Invert to get U^(-1)
212 15568 : CALL cp_cfm_triangular_invert(bmatrix)
213 : ! Reduce to get U^(-T) * H * U^(-1)
214 15568 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
215 15568 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
216 : ! Diagonalize
217 15568 : CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
218 : ! Restore vectors C = U^(-1) * C*
219 15568 : CALL cp_cfm_triangular_multiply(bmatrix, work)
220 : #if defined(__DLAF)
221 : END IF
222 : #endif
223 :
224 15568 : CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
225 305694 : eigenvalues(1:nmo) = evals(1:nmo)
226 :
227 15568 : DEALLOCATE (evals)
228 :
229 15568 : CALL timestop(handle)
230 :
231 15568 : END SUBROUTINE cp_cfm_geeig
232 :
233 : ! **************************************************************************************************
234 : !> \brief General Eigenvalue Problem AX = BXE
235 : !> Use canonical orthogonalization
236 : !> \param amatrix ...
237 : !> \param bmatrix ...
238 : !> \param eigenvectors ...
239 : !> \param eigenvalues ...
240 : !> \param work ...
241 : !> \param epseig ...
242 : ! **************************************************************************************************
243 494 : SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
244 :
245 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
246 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
247 : TYPE(cp_cfm_type), INTENT(IN) :: work
248 : REAL(KIND=dp), INTENT(IN) :: epseig
249 :
250 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
251 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
252 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
253 :
254 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
255 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
256 : nmo, nx
257 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
258 :
259 494 : CALL timeset(routineN, handle)
260 :
261 : ! Test sizes
262 494 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
263 494 : nmo = SIZE(eigenvalues)
264 2470 : ALLOCATE (evals(nao), cevals(nao))
265 :
266 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
267 494 : CALL cp_cfm_scale(-cone, bmatrix)
268 494 : CALL cp_cfm_heevd(bmatrix, work, evals)
269 19962 : evals(:) = -evals(:)
270 494 : nc = nao
271 19398 : DO i = 1, nao
272 19398 : IF (evals(i) < epseig) THEN
273 72 : nc = i - 1
274 72 : EXIT
275 : END IF
276 : END DO
277 494 : CPASSERT(nc /= 0)
278 :
279 494 : IF (nc /= nao) THEN
280 72 : IF (nc < nmo) THEN
281 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
282 0 : ncol = nmo - nc
283 0 : CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
284 : END IF
285 : ! Set NULL space in eigenvector matrix of S to zero
286 636 : DO icol = nc + 1, nao
287 80028 : DO irow = 1, nao
288 79956 : CALL cp_cfm_set_element(work, irow, icol, czero)
289 : END DO
290 : END DO
291 : ! Set small eigenvalues to a dummy save value
292 636 : evals(nc + 1:nao) = 1.0_dp
293 : END IF
294 : ! Calculate U*s**(-1/2)
295 19962 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
296 494 : CALL cp_cfm_column_scale(work, cevals)
297 : ! Reduce to get U^(-C) * H * U^(-1)
298 494 : CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
299 494 : CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
300 494 : IF (nc /= nao) THEN
301 : ! set diagonal values to save large value
302 636 : DO icol = nc + 1, nao
303 636 : CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp))
304 : END DO
305 : END IF
306 : ! Diagonalize
307 494 : CALL cp_cfm_heevd(amatrix, bmatrix, evals)
308 7138 : eigenvalues(1:nmo) = evals(1:nmo)
309 494 : nx = MIN(nc, nmo)
310 : ! Restore vectors C = U^(-1) * C*
311 494 : CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
312 :
313 494 : DEALLOCATE (evals)
314 :
315 494 : CALL timestop(handle)
316 :
317 988 : END SUBROUTINE cp_cfm_geeig_canon
318 :
319 : END MODULE cp_cfm_diag
|
