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