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