Line data Source code
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 : ! **************************************************************************************************
14 : MODULE cp_cfm_diag
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
17 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, &
18 : cp_cfm_column_scale, &
19 : cp_cfm_scale, &
20 : cp_cfm_triangular_invert, &
21 : cp_cfm_triangular_multiply
22 : USE cp_cfm_types, ONLY: cp_cfm_create, &
23 : cp_cfm_get_info, &
24 : cp_cfm_release, &
25 : cp_cfm_set_element, &
26 : cp_cfm_to_cfm, &
27 : cp_cfm_type
28 : USE cp_fm_diag, ONLY: diag_check_requested, &
29 : diag_check_warning_threshold, &
30 : diag_type, &
31 : direct_generalized_diagonalization, &
32 : cusolver_n_min, &
33 : FM_DIAG_TYPE_CUSOLVER, &
34 : FM_DIAG_TYPE_SCALAPACK
35 : USE cp_fm_cusolver_api, ONLY: cp_cfm_general_cusolver
36 : #if defined(__DLAF)
37 : USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
38 : cp_cfm_diag_dlaf
39 : USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
40 : USE cp_fm_diag, ONLY: dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
41 : #endif
42 : USE cp_log_handling, ONLY: cp_to_string
43 : USE kinds, ONLY: default_string_length, &
44 : dp
45 : USE machine, ONLY: default_output_unit
46 : USE mathconstants, ONLY: z_one, &
47 : z_zero
48 : #if defined (__HAS_IEEE_EXCEPTIONS)
49 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
50 : ieee_set_halting_mode, &
51 : IEEE_ALL
52 : #endif
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 :
58 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
59 :
60 : PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief Perform a diagonalisation of a complex matrix
66 : !> \param matrix ...
67 : !> \param eigenvectors ...
68 : !> \param eigenvalues ...
69 : !> \par History
70 : !> 12.2024 Added DLA-Future support [Rocco Meli]
71 : !> \author Joost VandeVondele
72 : ! **************************************************************************************************
73 47424 : SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
74 :
75 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
76 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd'
79 :
80 : INTEGER :: handle
81 :
82 47424 : CALL timeset(routineN, handle)
83 :
84 : #if defined(__DLAF)
85 : IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
86 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
87 : CALL cp_dlaf_initialize()
88 :
89 : ! Create DLAF grid from BLACS context; if already present, does nothing
90 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
91 :
92 : CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
93 : ELSE
94 : #endif
95 47424 : CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
96 : #if defined(__DLAF)
97 : END IF
98 : #endif
99 :
100 47424 : CALL timestop(handle)
101 :
102 47424 : END SUBROUTINE cp_cfm_heevd
103 :
104 : ! **************************************************************************************************
105 : !> \brief Perform a diagonalisation of a complex matrix
106 : !> \param matrix ...
107 : !> \param eigenvectors ...
108 : !> \param eigenvalues ...
109 : !> \par History
110 : !> - (De)Allocation checks updated (15.02.2011,MK)
111 : !> \author Joost VandeVondele
112 : ! **************************************************************************************************
113 47424 : SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
114 :
115 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
116 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
117 :
118 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
119 :
120 47424 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: work
121 : COMPLEX(KIND=dp), DIMENSION(:, :), &
122 47424 : POINTER :: m
123 : INTEGER :: handle, info, liwork, &
124 : lrwork, lwork, n
125 47424 : INTEGER, DIMENSION(:), POINTER :: iwork
126 47424 : REAL(KIND=dp), DIMENSION(:), POINTER :: rwork
127 : #if defined(__parallel)
128 : INTEGER, DIMENSION(9) :: descm, descv
129 : COMPLEX(KIND=dp), DIMENSION(:, :), &
130 47424 : POINTER :: v
131 : #endif
132 : #if defined (__HAS_IEEE_EXCEPTIONS)
133 : LOGICAL, DIMENSION(5) :: halt
134 : #endif
135 :
136 47424 : CALL timeset(routineN, handle)
137 :
138 47424 : n = matrix%matrix_struct%nrow_global
139 47424 : m => matrix%local_data
140 47424 : ALLOCATE (iwork(1), rwork(1), work(1))
141 : ! work space query
142 47424 : lwork = -1
143 47424 : lrwork = -1
144 47424 : liwork = -1
145 :
146 : #if defined(__parallel)
147 47424 : v => eigenvectors%local_data
148 474240 : descm(:) = matrix%matrix_struct%descriptor(:)
149 474240 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
150 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
151 47424 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
152 : ! The work space query for lwork does not return always sufficiently large values.
153 : ! Let's add some margin to avoid crashes.
154 47424 : lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
155 : ! needed to correct for a bug in scalapack, unclear how much the right number is
156 47424 : lrwork = CEILING(rwork(1)) + 1000000
157 47424 : liwork = iwork(1)
158 : #else
159 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
160 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
161 : lwork = CEILING(REAL(work(1), KIND=dp))
162 : lrwork = CEILING(rwork(1))
163 : liwork = iwork(1)
164 : #endif
165 :
166 47424 : DEALLOCATE (iwork, rwork, work)
167 331968 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
168 :
169 : ! (Sca-)LAPACK takes advantage of IEEE754 exceptions for speedup.
170 : ! Therefore, we disable floating point traps temporarily.
171 : #if defined (__HAS_IEEE_EXCEPTIONS)
172 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
173 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
174 : #endif
175 : #if defined(__parallel)
176 : CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
177 47424 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
178 : #else
179 : CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
180 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
181 : eigenvectors%local_data = matrix%local_data
182 : #endif
183 : #if defined (__HAS_IEEE_EXCEPTIONS)
184 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
185 : #endif
186 :
187 47424 : DEALLOCATE (iwork, rwork, work)
188 47424 : IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
189 :
190 47424 : CALL timestop(handle)
191 :
192 47424 : END SUBROUTINE cp_cfm_heevd_base
193 :
194 : ! **************************************************************************************************
195 : !> \brief Check C^H*S*C = I for a generalized complex eigenvalue problem.
196 : !> \param overlap original overlap matrix S; used as work matrix and overwritten
197 : !> \param eigenvectors eigenvectors C to be checked
198 : !> \param scratch work matrix
199 : !> \param nvec ...
200 : ! **************************************************************************************************
201 16 : SUBROUTINE check_generalized_diag(overlap, eigenvectors, scratch, nvec)
202 :
203 : TYPE(cp_cfm_type), INTENT(IN) :: eigenvectors
204 : TYPE(cp_cfm_type), INTENT(INOUT) :: overlap, scratch
205 : INTEGER, INTENT(IN) :: nvec
206 :
207 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_generalized_diag'
208 :
209 : CHARACTER(LEN=default_string_length) :: diag_type_name
210 : COMPLEX(KIND=dp) :: gold, test
211 : INTEGER :: handle, i, j, ncol, nrow, output_unit
212 : REAL(KIND=dp) :: eps, eps_abort, eps_warning
213 : #if defined(__parallel)
214 : TYPE(cp_blacs_env_type), POINTER :: context
215 : INTEGER :: il, jl, ipcol, iprow, &
216 : mypcol, myprow, npcol, nprow
217 : INTEGER, DIMENSION(9) :: desca
218 : #endif
219 :
220 16 : CALL timeset(routineN, handle)
221 :
222 16 : IF (.NOT. diag_check_requested()) THEN
223 0 : CALL timestop(handle)
224 0 : RETURN
225 : END IF
226 :
227 16 : output_unit = default_output_unit
228 16 : eps_warning = diag_check_warning_threshold()
229 16 : eps_abort = 10.0_dp*eps_warning
230 :
231 16 : nrow = eigenvectors%matrix_struct%nrow_global
232 16 : ncol = MIN(eigenvectors%matrix_struct%ncol_global, nvec)
233 :
234 16 : CALL cp_cfm_gemm("N", "N", nrow, ncol, nrow, z_one, overlap, eigenvectors, z_zero, scratch)
235 16 : CALL cp_cfm_gemm("C", "N", ncol, ncol, nrow, z_one, eigenvectors, scratch, z_zero, overlap)
236 :
237 16 : gold = z_zero
238 16 : test = z_zero
239 16 : eps = 0.0_dp
240 :
241 : #if defined(__parallel)
242 16 : context => overlap%matrix_struct%context
243 16 : myprow = context%mepos(1)
244 16 : mypcol = context%mepos(2)
245 16 : nprow = context%num_pe(1)
246 16 : npcol = context%num_pe(2)
247 160 : desca(:) = overlap%matrix_struct%descriptor(:)
248 160 : outer: DO j = 1, ncol
249 1456 : DO i = 1, ncol
250 1296 : CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
251 1440 : IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
252 648 : gold = MERGE(z_zero, z_one, i /= j)
253 648 : test = overlap%local_data(il, jl)
254 648 : eps = ABS(test - gold)
255 648 : IF (eps > eps_warning) EXIT outer
256 : END IF
257 : END DO
258 : END DO outer
259 : #else
260 : outer: DO j = 1, ncol
261 : DO i = 1, ncol
262 : gold = MERGE(z_zero, z_one, i /= j)
263 : test = overlap%local_data(i, j)
264 : eps = ABS(test - gold)
265 : IF (eps > eps_warning) EXIT outer
266 : END DO
267 : END DO outer
268 : #endif
269 :
270 16 : IF (eps > eps_warning) THEN
271 0 : IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN
272 0 : diag_type_name = "HEGVX"
273 0 : ELSE IF (diag_type == FM_DIAG_TYPE_CUSOLVER) THEN
274 0 : diag_type_name = "CUSOLVER"
275 : #if defined(__DLAF)
276 : ELSE IF (diag_type == FM_DIAG_TYPE_DLAF) THEN
277 : diag_type_name = "DLAF"
278 : #endif
279 : ELSE
280 0 : diag_type_name = "generalized eigensolver"
281 : END IF
282 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/,T2,A,I0,A,I0,A,ES10.3,/,T2,A,F0.0,A,ES10.3)") &
283 0 : "The generalized eigenvectors returned by "//TRIM(diag_type_name)//" are not S-orthonormal", &
284 0 : "Absolute deviation of matrix element (", i, ", ", j, ") is ", eps, &
285 0 : "The deviation from the expected value ", REAL(gold, KIND=dp), " is", eps
286 0 : IF (eps > eps_abort) THEN
287 0 : CPABORT("ERROR in "//routineN//": Check of generalized matrix diagonalization failed")
288 : ELSE
289 0 : CPWARN("Check of generalized matrix diagonalization failed in routine "//routineN)
290 : END IF
291 : END IF
292 :
293 16 : CALL timestop(handle)
294 :
295 : END SUBROUTINE check_generalized_diag
296 :
297 : ! **************************************************************************************************
298 : !> \brief General Eigenvalue Problem AX = BXE
299 : !> Single option version: Cholesky decomposition of B
300 : !> \param amatrix ...
301 : !> \param bmatrix ...
302 : !> \param eigenvectors ...
303 : !> \param eigenvalues ...
304 : !> \param work ...
305 : !> \par History
306 : !> 12.2024 Added DLA-Future support [Rocco Meli]
307 : ! **************************************************************************************************
308 33882 : SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
309 :
310 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
311 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
312 : TYPE(cp_cfm_type), INTENT(IN) :: work
313 :
314 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig'
315 :
316 : INTEGER :: handle, nao, nmo
317 : LOGICAL :: check_eigenvectors
318 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
319 : TYPE(cp_cfm_type) :: overlap_check, scratch_check
320 :
321 33882 : CALL timeset(routineN, handle)
322 :
323 33882 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
324 101646 : ALLOCATE (evals(nao))
325 33882 : nmo = SIZE(eigenvalues)
326 33882 : check_eigenvectors = diag_check_requested()
327 :
328 33882 : IF (diag_type == FM_DIAG_TYPE_CUSOLVER .AND. direct_generalized_diagonalization .AND. &
329 : nao >= cusolver_n_min) THEN
330 : ! Use cuSolverMP generalized eigenvalue solver without a CP2K-side
331 : ! Cholesky reduction.
332 0 : IF (check_eigenvectors) THEN
333 0 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
334 0 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
335 0 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
336 : END IF
337 0 : CALL cp_cfm_general_cusolver(amatrix, bmatrix, work, evals)
338 0 : IF (check_eigenvectors) THEN
339 0 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
340 0 : CALL cp_cfm_release(scratch_check)
341 0 : CALL cp_cfm_release(overlap_check)
342 : END IF
343 : #if defined(__DLAF)
344 : ELSE IF (diag_type == FM_DIAG_TYPE_DLAF .AND. direct_generalized_diagonalization .AND. &
345 : nao >= dlaf_neigvec_min) THEN
346 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
347 : CALL cp_dlaf_initialize()
348 :
349 : ! Create DLAF grid from BLACS context; if already present, does nothing
350 : CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
351 : CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
352 : CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
353 :
354 : ! Use DLA-Future generalized eigenvalue solver for large matrices
355 : IF (check_eigenvectors) THEN
356 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
357 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
358 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
359 : END IF
360 : CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
361 : IF (check_eigenvectors) THEN
362 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
363 : CALL cp_cfm_release(scratch_check)
364 : CALL cp_cfm_release(overlap_check)
365 : END IF
366 : #endif
367 : #if defined(__parallel)
368 33882 : ELSE IF (diag_type == FM_DIAG_TYPE_SCALAPACK .AND. direct_generalized_diagonalization) THEN
369 : ! Use ScaLAPACK generalized eigenvalue solver without a CP2K-side
370 : ! Cholesky reduction.
371 16 : IF (check_eigenvectors) THEN
372 16 : CALL cp_cfm_create(overlap_check, bmatrix%matrix_struct)
373 16 : CALL cp_cfm_create(scratch_check, bmatrix%matrix_struct)
374 16 : CALL cp_cfm_to_cfm(bmatrix, overlap_check)
375 : END IF
376 16 : CALL cp_cfm_geeig_scalapack(amatrix, bmatrix, work, evals)
377 16 : IF (check_eigenvectors) THEN
378 16 : CALL check_generalized_diag(overlap_check, work, scratch_check, nmo)
379 16 : CALL cp_cfm_release(scratch_check)
380 16 : CALL cp_cfm_release(overlap_check)
381 : END IF
382 : #endif
383 : ELSE
384 : ! Cholesky decompose S=U(T)U
385 33866 : CALL cp_cfm_cholesky_decompose(bmatrix)
386 : ! Invert to get U^(-1)
387 33866 : CALL cp_cfm_triangular_invert(bmatrix)
388 : ! Reduce to get U^(-T) * H * U^(-1)
389 33866 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
390 33866 : CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
391 : ! Diagonalize
392 33866 : CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
393 : ! Restore vectors C = U^(-1) * C*
394 33866 : CALL cp_cfm_triangular_multiply(bmatrix, work)
395 : END IF
396 :
397 33882 : CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
398 738612 : eigenvalues(1:nmo) = evals(1:nmo)
399 :
400 33882 : DEALLOCATE (evals)
401 :
402 33882 : CALL timestop(handle)
403 :
404 33882 : END SUBROUTINE cp_cfm_geeig
405 :
406 : ! **************************************************************************************************
407 : !> \brief General Eigenvalue Problem AX = BXE using ScaLAPACK PZHEGVX.
408 : !> \param amatrix ...
409 : !> \param bmatrix ...
410 : !> \param eigenvectors ...
411 : !> \param eigenvalues ...
412 : ! **************************************************************************************************
413 16 : SUBROUTINE cp_cfm_geeig_scalapack(amatrix, bmatrix, eigenvectors, eigenvalues)
414 :
415 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
416 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
417 :
418 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_scalapack'
419 :
420 : #if defined(__parallel)
421 : REAL(KIND=dp), PARAMETER :: orfac = -1.0_dp, &
422 : vl = 0.0_dp, &
423 : vu = 0.0_dp
424 :
425 16 : COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
426 16 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
427 : INTEGER :: handle, info, liwork, lwork, lrwork, &
428 : m, n, nb, neig, npcol, nprow, nz
429 : INTEGER, DIMENSION(9) :: desca, descb, descz
430 16 : INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr, ifail, iwork
431 : REAL(KIND=dp) :: abstol
432 16 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: gap, rwork, w
433 :
434 : INTEGER :: mq0, nn, np0, npe
435 : INTEGER, EXTERNAL :: iceil, numroc
436 : REAL(KIND=dp), EXTERNAL :: dlamch
437 : #if defined (__HAS_IEEE_EXCEPTIONS)
438 : LOGICAL, DIMENSION(5) :: halt
439 : #endif
440 : #else
441 : INTEGER :: handle
442 : #endif
443 :
444 16 : CALL timeset(routineN, handle)
445 :
446 : #if defined(__parallel)
447 16 : n = amatrix%matrix_struct%nrow_global
448 16 : neig = MIN(SIZE(eigenvalues), n)
449 :
450 16 : IF (neig == 0) THEN
451 0 : CALL timestop(handle)
452 0 : RETURN
453 : END IF
454 :
455 16 : IF (amatrix%matrix_struct%nrow_block /= amatrix%matrix_struct%ncol_block) THEN
456 0 : CPABORT("ERROR in "//routineN//": Invalid blocksize (no square blocks) found")
457 : END IF
458 :
459 16 : a => amatrix%local_data
460 16 : b => bmatrix%local_data
461 16 : z => eigenvectors%local_data
462 160 : desca(:) = amatrix%matrix_struct%descriptor(:)
463 160 : descb(:) = bmatrix%matrix_struct%descriptor(:)
464 160 : descz(:) = eigenvectors%matrix_struct%descriptor(:)
465 :
466 16 : nprow = amatrix%matrix_struct%context%num_pe(1)
467 16 : npcol = amatrix%matrix_struct%context%num_pe(2)
468 16 : npe = nprow*npcol
469 16 : nb = amatrix%matrix_struct%nrow_block
470 16 : nn = MAX(n, nb, 2)
471 16 : np0 = numroc(nn, nb, 0, 0, nprow)
472 16 : mq0 = MAX(numroc(nn, nb, 0, 0, npcol), nb)
473 :
474 16 : lwork = n + (np0 + mq0 + nb)*nb
475 16 : lrwork = 4*n + MAX(5*nn, np0*mq0) + iceil(neig, npe)*nn + MAX(0, neig - 1)*n
476 16 : liwork = 6*MAX(n, npe + 1, 4)
477 :
478 48 : ALLOCATE (gap(npe))
479 48 : gap = 0.0_dp
480 48 : ALLOCATE (iclustr(2*npe))
481 80 : iclustr = 0
482 48 : ALLOCATE (ifail(n))
483 640 : ifail = 0
484 48 : ALLOCATE (iwork(liwork))
485 48 : ALLOCATE (rwork(lrwork))
486 48 : ALLOCATE (w(n))
487 48 : ALLOCATE (work(lwork))
488 :
489 16 : abstol = 2.0_dp*dlamch("S")
490 :
491 : #if defined (__HAS_IEEE_EXCEPTIONS)
492 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
493 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
494 : #endif
495 : CALL pzhegvx(1, "V", "I", "U", n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, &
496 : vl, vu, 1, neig, abstol, m, nz, w(1), orfac, z(1, 1), 1, 1, descz, &
497 : work(1), lwork, rwork(1), lrwork, iwork(1), liwork, ifail(1), &
498 16 : iclustr(1), gap(1), info)
499 : #if defined (__HAS_IEEE_EXCEPTIONS)
500 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
501 : #endif
502 :
503 16 : IF (info /= 0 .OR. m < neig .OR. nz < neig) THEN
504 0 : CPABORT("ERROR in PZHEGVX (ScaLAPACK), info="//TRIM(cp_to_string(info)))
505 : END IF
506 :
507 640 : eigenvalues(:) = 0.0_dp
508 640 : eigenvalues(1:neig) = w(1:neig)
509 :
510 16 : DEALLOCATE (gap, iclustr, ifail, iwork, rwork, w, work)
511 : #else
512 : MARK_USED(amatrix)
513 : MARK_USED(bmatrix)
514 : MARK_USED(eigenvectors)
515 : MARK_USED(eigenvalues)
516 : CPABORT("ERROR in "//routineN//": PZHEGVX requested without ScaLAPACK support")
517 : #endif
518 :
519 16 : CALL timestop(handle)
520 :
521 16 : END SUBROUTINE cp_cfm_geeig_scalapack
522 :
523 : ! **************************************************************************************************
524 : !> \brief General Eigenvalue Problem AX = BXE
525 : !> Use canonical orthogonalization
526 : !> \param amatrix ...
527 : !> \param bmatrix ...
528 : !> \param eigenvectors ...
529 : !> \param eigenvalues ...
530 : !> \param work ...
531 : !> \param epseig ...
532 : ! **************************************************************************************************
533 2232 : SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
534 :
535 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
536 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
537 : TYPE(cp_cfm_type), INTENT(IN) :: work
538 : REAL(KIND=dp), INTENT(IN) :: epseig
539 :
540 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
541 :
542 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
543 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
544 : nmo, nx
545 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
546 :
547 2232 : CALL timeset(routineN, handle)
548 :
549 : ! Test sizes
550 2232 : CALL cp_cfm_get_info(amatrix, nrow_global=nao)
551 2232 : nmo = SIZE(eigenvalues)
552 11160 : ALLOCATE (evals(nao), cevals(nao))
553 :
554 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
555 2232 : CALL cp_cfm_scale(-z_one, bmatrix)
556 2232 : CALL cp_cfm_heevd(bmatrix, work, evals)
557 62970 : evals(:) = -evals(:)
558 2232 : nc = nao
559 62406 : DO i = 1, nao
560 62406 : IF (evals(i) < epseig) THEN
561 72 : nc = i - 1
562 72 : EXIT
563 : END IF
564 : END DO
565 2232 : CPASSERT(nc /= 0)
566 :
567 2232 : IF (nc /= nao) THEN
568 72 : IF (nc < nmo) THEN
569 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
570 0 : ncol = nmo - nc
571 0 : CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
572 : END IF
573 : ! Set NULL space in eigenvector matrix of S to zero
574 636 : DO icol = nc + 1, nao
575 80028 : DO irow = 1, nao
576 79956 : CALL cp_cfm_set_element(work, irow, icol, z_zero)
577 : END DO
578 : END DO
579 : ! Set small eigenvalues to a dummy save value
580 636 : evals(nc + 1:nao) = 1.0_dp
581 : END IF
582 : ! Calculate U*s**(-1/2)
583 62970 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
584 2232 : CALL cp_cfm_column_scale(work, cevals)
585 : ! Reduce to get U^(-C) * H * U^(-1)
586 2232 : CALL cp_cfm_gemm("C", "N", nao, nao, nao, z_one, work, amatrix, z_zero, bmatrix)
587 2232 : CALL cp_cfm_gemm("N", "N", nao, nao, nao, z_one, bmatrix, work, z_zero, amatrix)
588 2232 : IF (nc /= nao) THEN
589 : ! set diagonal values to save large value
590 636 : DO icol = nc + 1, nao
591 636 : CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp))
592 : END DO
593 : END IF
594 : ! Diagonalize
595 2232 : CALL cp_cfm_heevd(amatrix, bmatrix, evals)
596 29248 : eigenvalues(1:nmo) = evals(1:nmo)
597 2232 : nx = MIN(nc, nmo)
598 : ! Restore vectors C = U^(-1) * C*
599 2232 : CALL cp_cfm_gemm("N", "N", nao, nx, nc, z_one, work, bmatrix, z_zero, eigenvectors)
600 :
601 2232 : DEALLOCATE (evals)
602 :
603 2232 : CALL timestop(handle)
604 :
605 4464 : END SUBROUTINE cp_cfm_geeig_canon
606 :
607 : END MODULE cp_cfm_diag
|
