(git:d18deda)
Loading...
Searching...
No Matches
cp_cfm_diag.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
46CONTAINS
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 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 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 CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
74#if defined(__DLAF)
75 END IF
76#endif
77
78 CALL timestop(handle)
79
80 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 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 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
99 COMPLEX(KIND=dp), DIMENSION(:, :), &
100 POINTER :: m
101 INTEGER :: handle, info, liwork, &
102 lrwork, lwork, n
103 INTEGER, DIMENSION(:), POINTER :: iwork
104 REAL(kind=dp), DIMENSION(:), POINTER :: rwork
105#if defined(__parallel)
106 INTEGER, DIMENSION(9) :: descm, descv
107 COMPLEX(KIND=dp), DIMENSION(:, :), &
108 POINTER :: v
109#if defined (__HAS_IEEE_EXCEPTIONS)
110 LOGICAL, DIMENSION(5) :: halt
111#endif
112#endif
113
114 CALL timeset(routinen, handle)
115
116 n = matrix%matrix_struct%nrow_global
117 m => matrix%local_data
118 ALLOCATE (iwork(1), rwork(1), work(1))
119 ! work space query
120 lwork = -1
121 lrwork = -1
122 liwork = -1
123
124#if defined(__parallel)
125 v => eigenvectors%local_data
126 descm(:) = matrix%matrix_struct%descriptor(:)
127 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 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 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 lrwork = ceiling(rwork(1)) + 1000000
135 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 DEALLOCATE (iwork, rwork, work)
145 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 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 DEALLOCATE (iwork, rwork, work)
168 IF (info /= 0) &
169 cpabort("Diagonalisation of a complex matrix failed")
170
171 CALL timestop(handle)
172
173 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 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 CALL timeset(routinen, handle)
198
199 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
200 ALLOCATE (evals(nao))
201 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 CALL cp_cfm_cholesky_decompose(bmatrix)
211 ! Invert to get U^(-1)
212 CALL cp_cfm_triangular_invert(bmatrix)
213 ! Reduce to get U^(-T) * H * U^(-1)
214 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
215 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
216 ! Diagonalize
217 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
218 ! Restore vectors C = U^(-1) * C*
219 CALL cp_cfm_triangular_multiply(bmatrix, work)
220#if defined(__DLAF)
221 END IF
222#endif
223
224 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
225 eigenvalues(1:nmo) = evals(1:nmo)
226
227 DEALLOCATE (evals)
228
229 CALL timestop(handle)
230
231 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 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 CALL timeset(routinen, handle)
260
261 ! Test sizes
262 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
263 nmo = SIZE(eigenvalues)
264 ALLOCATE (evals(nao), cevals(nao))
265
266 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
267 CALL cp_cfm_scale(-cone, bmatrix)
268 CALL cp_cfm_heevd(bmatrix, work, evals)
269 evals(:) = -evals(:)
270 nc = nao
271 DO i = 1, nao
272 IF (evals(i) < epseig) THEN
273 nc = i - 1
274 EXIT
275 END IF
276 END DO
277 cpassert(nc /= 0)
278
279 IF (nc /= nao) THEN
280 IF (nc < nmo) THEN
281 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
282 ncol = nmo - nc
283 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 DO icol = nc + 1, nao
287 DO irow = 1, nao
288 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 evals(nc + 1:nao) = 1.0_dp
293 END IF
294 ! Calculate U*s**(-1/2)
295 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
296 CALL cp_cfm_column_scale(work, cevals)
297 ! Reduce to get U^(-C) * H * U^(-1)
298 CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
299 CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
300 IF (nc /= nao) THEN
301 ! set diagonal values to save large value
302 DO icol = nc + 1, nao
303 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 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
308 eigenvalues(1:nmo) = evals(1:nmo)
309 nx = min(nc, nmo)
310 ! Restore vectors C = U^(-1) * C*
311 CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
312
313 DEALLOCATE (evals)
314
315 CALL timestop(handle)
316
317 END SUBROUTINE cp_cfm_geeig_canon
318
319END MODULE cp_cfm_diag
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:55
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:563
cp_cfm_basic_linalg::cp_cfm_triangular_multiply
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_cfm_basic_linalg.F:979
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:633
cp_cfm_basic_linalg::cp_cfm_triangular_invert
subroutine, public cp_cfm_triangular_invert(matrix_a, uplo, info_out)
Inverts a triangular matrix.
Definition cp_cfm_basic_linalg.F:1061
cp_cfm_cholesky
various cholesky decomposition related routines
Definition cp_cfm_cholesky.F:13
cp_cfm_cholesky::cp_cfm_cholesky_decompose
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Definition cp_cfm_cholesky.F:51
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_geeig
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_cfm_diag.F:187
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:58
cp_cfm_diag::cp_cfm_geeig_canon
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Definition cp_cfm_diag.F:244
cp_cfm_dlaf_api
Definition cp_cfm_dlaf_api.F:8
cp_cfm_dlaf_api::cp_cfm_diag_dlaf
subroutine, public cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
DLA-Future eigensolver for complex Hermitian matrices.
Definition cp_cfm_dlaf_api.F:78
cp_cfm_dlaf_api::cp_cfm_diag_gen_dlaf
subroutine, public cp_cfm_diag_gen_dlaf(amatrix, bmatrix, eigenvectors, eigenvalues)
DLA-Future generalized eigensolver for complex Hermitian matrices.
Definition cp_cfm_dlaf_api.F:116
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition cp_cfm_types.F:275
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:603
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::fm_diag_type_dlaf
integer, parameter, public fm_diag_type_dlaf
Definition cp_fm_diag.F:104
cp_fm_diag::diag_type
integer, save, public diag_type
Definition cp_fm_diag.F:90
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68