(git:e68d6fd)
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-2026 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#endif
118#if defined (__HAS_IEEE_EXCEPTIONS)
119 LOGICAL, DIMENSION(5) :: halt
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! (Sca-)LAPACK 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#if defined(__parallel)
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#else
165 CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
166 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
167 eigenvectors%local_data = matrix%local_data
168#endif
169#if defined (__HAS_IEEE_EXCEPTIONS)
170 CALL ieee_set_halting_mode(ieee_all, halt)
171#endif
172
173 DEALLOCATE (iwork, rwork, work)
174 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
175
176 CALL timestop(handle)
177
178 END SUBROUTINE cp_cfm_heevd_base
179
180! **************************************************************************************************
181!> \brief General Eigenvalue Problem AX = BXE
182!> Single option version: Cholesky decomposition of B
183!> \param amatrix ...
184!> \param bmatrix ...
185!> \param eigenvectors ...
186!> \param eigenvalues ...
187!> \param work ...
188!> \par History
189!> 12.2024 Added DLA-Future support [Rocco Meli]
190! **************************************************************************************************
191 SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
192
193 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
194 REAL(kind=dp), DIMENSION(:) :: eigenvalues
195 TYPE(cp_cfm_type), INTENT(IN) :: work
196
197 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig'
198
199 INTEGER :: handle, nao, nmo
200 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
201
202 CALL timeset(routinen, handle)
203
204 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
205 ALLOCATE (evals(nao))
206 nmo = SIZE(eigenvalues)
207
208#if defined(__DLAF)
209 IF (diag_type == fm_diag_type_dlaf .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
210 ! Initialize DLA-Future on-demand; if already initialized, does nothing
211 CALL cp_dlaf_initialize()
212
213 ! Create DLAF grid from BLACS context; if already present, does nothing
214 CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
215 CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
216 CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
217
218 ! Use DLA-Future generalized eigenvalue solver for large matrices
219 CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
220 ELSE
221#endif
222 ! Cholesky decompose S=U(T)U
223 CALL cp_cfm_cholesky_decompose(bmatrix)
224 ! Invert to get U^(-1)
225 CALL cp_cfm_triangular_invert(bmatrix)
226 ! Reduce to get U^(-T) * H * U^(-1)
227 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
228 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
229 ! Diagonalize
230 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
231 ! Restore vectors C = U^(-1) * C*
232 CALL cp_cfm_triangular_multiply(bmatrix, work)
233#if defined(__DLAF)
234 END IF
235#endif
236
237 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
238 eigenvalues(1:nmo) = evals(1:nmo)
239
240 DEALLOCATE (evals)
241
242 CALL timestop(handle)
243
244 END SUBROUTINE cp_cfm_geeig
245
246! **************************************************************************************************
247!> \brief General Eigenvalue Problem AX = BXE
248!> Use canonical orthogonalization
249!> \param amatrix ...
250!> \param bmatrix ...
251!> \param eigenvectors ...
252!> \param eigenvalues ...
253!> \param work ...
254!> \param epseig ...
255! **************************************************************************************************
256 SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
257
258 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
259 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
260 TYPE(cp_cfm_type), INTENT(IN) :: work
261 REAL(kind=dp), INTENT(IN) :: epseig
262
263 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig_canon'
264
265 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
266 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
267 nmo, nx
268 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
269
270 CALL timeset(routinen, handle)
271
272 ! Test sizes
273 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
274 nmo = SIZE(eigenvalues)
275 ALLOCATE (evals(nao), cevals(nao))
276
277 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
278 CALL cp_cfm_scale(-z_one, bmatrix)
279 CALL cp_cfm_heevd(bmatrix, work, evals)
280 evals(:) = -evals(:)
281 nc = nao
282 DO i = 1, nao
283 IF (evals(i) < epseig) THEN
284 nc = i - 1
285 EXIT
286 END IF
287 END DO
288 cpassert(nc /= 0)
289
290 IF (nc /= nao) THEN
291 IF (nc < nmo) THEN
292 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
293 ncol = nmo - nc
294 CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
295 END IF
296 ! Set NULL space in eigenvector matrix of S to zero
297 DO icol = nc + 1, nao
298 DO irow = 1, nao
299 CALL cp_cfm_set_element(work, irow, icol, z_zero)
300 END DO
301 END DO
302 ! Set small eigenvalues to a dummy save value
303 evals(nc + 1:nao) = 1.0_dp
304 END IF
305 ! Calculate U*s**(-1/2)
306 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
307 CALL cp_cfm_column_scale(work, cevals)
308 ! Reduce to get U^(-C) * H * U^(-1)
309 CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
310 CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
311 IF (nc /= nao) THEN
312 ! set diagonal values to save large value
313 DO icol = nc + 1, nao
314 CALL cp_cfm_set_element(amatrix, icol, icol, cmplx(10000.0_dp, 0.0_dp, kind=dp))
315 END DO
316 END IF
317 ! Diagonalize
318 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
319 eigenvalues(1:nmo) = evals(1:nmo)
320 nx = min(nc, nmo)
321 ! Restore vectors C = U^(-1) * C*
322 CALL cp_cfm_gemm("N", "N", nao, nx, nc, z_one, work, bmatrix, z_zero, eigenvectors)
323
324 DEALLOCATE (evals)
325
326 CALL timestop(handle)
327
328 END SUBROUTINE cp_cfm_geeig_canon
329
330END 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:56
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:192
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:257
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:302
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:630
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:69