(git:374b731)
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-2024 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_basic_linalg, ONLY: cp_cfm_cholesky_decompose, &
16 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 USE kinds, ONLY: dp
26#if defined (__HAS_IEEE_EXCEPTIONS)
27 USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
28 ieee_set_halting_mode, &
29 ieee_all
30#endif
31
32#include "../base/base_uses.f90"
33
34 IMPLICIT NONE
35 PRIVATE
36 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
37
38 PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
39
40CONTAINS
41
42! **************************************************************************************************
43!> \brief Perform a diagonalisation of a complex matrix
44!> \param matrix ...
45!> \param eigenvectors ...
46!> \param eigenvalues ...
47!> \par History
48!> - (De)Allocation checks updated (15.02.2011,MK)
49!> \author Joost VandeVondele
50! **************************************************************************************************
51 SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
52
53 TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
54 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
55
56 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_heevd'
57
58 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
59 COMPLEX(KIND=dp), DIMENSION(:, :), &
60 POINTER :: m
61 INTEGER :: handle, info, liwork, &
62 lrwork, lwork, n
63 INTEGER, DIMENSION(:), POINTER :: iwork
64 REAL(kind=dp), DIMENSION(:), POINTER :: rwork
65#if defined(__SCALAPACK)
66 INTEGER, DIMENSION(9) :: descm, descv
67 COMPLEX(KIND=dp), DIMENSION(:, :), &
68 POINTER :: v
69#if defined (__HAS_IEEE_EXCEPTIONS)
70 LOGICAL, DIMENSION(5) :: halt
71#endif
72#endif
73
74 CALL timeset(routinen, handle)
75
76 n = matrix%matrix_struct%nrow_global
77 m => matrix%local_data
78 ALLOCATE (iwork(1), rwork(1), work(1))
79 ! work space query
80 lwork = -1
81 lrwork = -1
82 liwork = -1
83
84#if defined(__SCALAPACK)
85 v => eigenvectors%local_data
86 descm(:) = matrix%matrix_struct%descriptor(:)
87 descv(:) = eigenvectors%matrix_struct%descriptor(:)
88 CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
89 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
90 ! The work space query for lwork does not return always sufficiently large values.
91 ! Let's add some margin to avoid crashes.
92 lwork = ceiling(real(work(1), kind=dp)) + 1000
93 ! needed to correct for a bug in scalapack, unclear how much the right number is
94 lrwork = ceiling(real(rwork(1), kind=dp)) + 1000000
95 liwork = iwork(1)
96#else
97 CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
98 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
99 lwork = ceiling(real(work(1), kind=dp))
100 lrwork = ceiling(real(rwork(1), kind=dp))
101 liwork = iwork(1)
102#endif
103
104 DEALLOCATE (iwork, rwork, work)
105 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
106
107#if defined(__SCALAPACK)
108! Scalapack takes advantage of IEEE754 exceptions for speedup.
109! Therefore, we disable floating point traps temporarily.
110#if defined (__HAS_IEEE_EXCEPTIONS)
111 CALL ieee_get_halting_mode(ieee_all, halt)
112 CALL ieee_set_halting_mode(ieee_all, .false.)
113#endif
114
115 CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
116 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
117
118#if defined (__HAS_IEEE_EXCEPTIONS)
119 CALL ieee_set_halting_mode(ieee_all, halt)
120#endif
121#else
122 CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
123 work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
124 eigenvectors%local_data = matrix%local_data
125#endif
126
127 DEALLOCATE (iwork, rwork, work)
128 IF (info /= 0) &
129 cpabort("Diagonalisation of a complex matrix failed")
130
131 CALL timestop(handle)
132
133 END SUBROUTINE cp_cfm_heevd
134
135! **************************************************************************************************
136!> \brief General Eigenvalue Problem AX = BXE
137!> Single option version: Cholesky decomposition of B
138!> \param amatrix ...
139!> \param bmatrix ...
140!> \param eigenvectors ...
141!> \param eigenvalues ...
142!> \param work ...
143! **************************************************************************************************
144 SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
145
146 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
147 REAL(kind=dp), DIMENSION(:) :: eigenvalues
148 TYPE(cp_cfm_type), INTENT(IN) :: work
149
150 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig'
151
152 INTEGER :: handle, nao, nmo
153 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
154
155 CALL timeset(routinen, handle)
156
157 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
158 ALLOCATE (evals(nao))
159 ! Cholesky decompose S=U(T)U
160 CALL cp_cfm_cholesky_decompose(bmatrix)
161 ! Invert to get U^(-1)
162 CALL cp_cfm_triangular_invert(bmatrix)
163 ! Reduce to get U^(-T) * H * U^(-1)
164 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
165 CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
166 ! Diagonalize
167 CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
168 ! Restore vectors C = U^(-1) * C*
169 CALL cp_cfm_triangular_multiply(bmatrix, work)
170 nmo = SIZE(eigenvalues)
171 CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
172 eigenvalues(1:nmo) = evals(1:nmo)
173
174 DEALLOCATE (evals)
175
176 CALL timestop(handle)
177
178 END SUBROUTINE cp_cfm_geeig
179
180! **************************************************************************************************
181!> \brief General Eigenvalue Problem AX = BXE
182!> Use canonical orthogonalization
183!> \param amatrix ...
184!> \param bmatrix ...
185!> \param eigenvectors ...
186!> \param eigenvalues ...
187!> \param work ...
188!> \param epseig ...
189! **************************************************************************************************
190 SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
191
192 TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
193 REAL(kind=dp), DIMENSION(:) :: eigenvalues
194 TYPE(cp_cfm_type), INTENT(IN) :: work
195 REAL(kind=dp), INTENT(IN) :: epseig
196
197 CHARACTER(len=*), PARAMETER :: routinen = 'cp_cfm_geeig_canon'
198 COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
199 czero = cmplx(0.0_dp, 0.0_dp, kind=dp)
200
201 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
202 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
203 nmo, nx
204 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
205
206 CALL timeset(routinen, handle)
207
208 ! Test sizes
209 CALL cp_cfm_get_info(amatrix, nrow_global=nao)
210 nmo = SIZE(eigenvalues)
211 ALLOCATE (evals(nao), cevals(nao))
212
213 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
214 CALL cp_cfm_scale(-cone, bmatrix)
215 CALL cp_cfm_heevd(bmatrix, work, evals)
216 evals(:) = -evals(:)
217 nc = nao
218 DO i = 1, nao
219 IF (evals(i) < epseig) THEN
220 nc = i - 1
221 EXIT
222 END IF
223 END DO
224 cpassert(nc /= 0)
225
226 IF (nc /= nao) THEN
227 IF (nc < nmo) THEN
228 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
229 ncol = nmo - nc
230 CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
231 END IF
232 ! Set NULL space in eigenvector matrix of S to zero
233 DO icol = nc + 1, nao
234 DO irow = 1, nao
235 CALL cp_cfm_set_element(work, irow, icol, czero)
236 END DO
237 END DO
238 ! Set small eigenvalues to a dummy save value
239 evals(nc + 1:nao) = 1.0_dp
240 END IF
241 ! calculate U*s**(-1/2)
242 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
243 CALL cp_cfm_column_scale(work, cevals)
244 ! Reduce to get U^(-C) * H * U^(-1)
245 CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
246 CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
247 IF (nc /= nao) THEN
248 ! set diagonal values to save large value
249 DO icol = nc + 1, nao
250 CALL cp_cfm_set_element(amatrix, icol, icol, cmplx(10000.0_dp, 0.0_dp, kind=dp))
251 END DO
252 END IF
253 ! Diagonalize
254 CALL cp_cfm_heevd(amatrix, bmatrix, evals)
255 eigenvalues(1:nmo) = evals(1:nmo)
256 nx = min(nc, nmo)
257 ! Restore vectors C = U^(-1) * C*
258 CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
259
260 DEALLOCATE (evals)
261
262 CALL timestop(handle)
263
264 END SUBROUTINE cp_cfm_geeig_canon
265
266END MODULE cp_cfm_diag
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:61
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:567
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:1107
cp_cfm_basic_linalg::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_basic_linalg.F:926
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:649
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:1189
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:145
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:52
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:191
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:279
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:607
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