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