(git:374b731)
Loading...
Searching...
No Matches
cp_fm_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 some of the diagonalization schemes available for
10!> cp_fm_type. cp_fm_power also moved here as it is very related
11!> \note
12!> first version : most routines imported
13!> \par History
14!> - unused Jacobi routines removed, cosmetics (05.04.06,MK)
15!> \author Joost VandeVondele (2003-08)
16! **************************************************************************************************
17MODULE cp_fm_diag
18
19 USE cp_blacs_types, ONLY: cp_blacs_type
20 USE cp_blacs_env, ONLY: cp_blacs_env_type
21 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale, &
22 cp_fm_gemm, &
23 cp_fm_scale, &
24 cp_fm_syrk, &
25 cp_fm_triangular_invert, &
26 cp_fm_triangular_multiply, &
27 cp_fm_upper_to_full
28 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
29 USE cp_fm_diag_utils, ONLY: cp_fm_redistribute_end, &
30 cp_fm_redistribute_start
31 USE cp_fm_elpa, ONLY: cp_fm_diag_elpa, &
32 finalize_elpa_library, &
33 initialize_elpa_library, &
34 set_elpa_kernel, &
35 set_elpa_print, &
36 set_elpa_qr
37 USE cp_fm_cusolver_api, ONLY: cp_fm_diag_cusolver
38#if defined(__DLAF)
39 USE cp_fm_dlaf_api, ONLY: cp_fm_diag_dlaf
40#endif
41 USE cp_fm_types, ONLY: cp_fm_get_info, &
42 cp_fm_set_element, &
43 cp_fm_to_fm, &
44 cp_fm_type
45 USE cp_log_handling, ONLY: cp_get_default_logger, &
46 cp_logger_get_default_unit_nr, &
47 cp_logger_get_unit_nr, &
48 cp_logger_type
49 USE kinds, ONLY: default_string_length, &
50 dp
51 USE machine, ONLY: default_output_unit, &
52 m_memory
53#if defined (__SCALAPACK)
54 USE message_passing, ONLY: mp_comm_type
55#endif
56#if defined (__HAS_IEEE_EXCEPTIONS)
57 USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
58 ieee_set_halting_mode, &
59 ieee_all
60#endif
61#include "../base/base_uses.f90"
62
63 IMPLICIT NONE
64
65 PRIVATE
66
67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_diag'
68
69 REAL(kind=dp), PARAMETER, PUBLIC :: eps_check_diag_default = 5.0e-14_dp
70
71 ! The following saved variables are diagonalization global
72 ! Stores the default library for diagonalization
73 INTEGER, SAVE :: diag_type = 0
74 ! Minimum number of eigenvectors for the use of the ELPA eigensolver.
75 ! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
76 INTEGER, SAVE :: elpa_neigvec_min = 0
77#if defined(__DLAF)
78 ! Minimum number of eigenvectors for the use of the DLAF eigensolver.
79 ! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
80 INTEGER, SAVE :: dlaf_neigvec_min = 0
81#endif
82 ! Threshold value for the orthonormality check of the eigenvectors obtained
83 ! after a diagonalization. A negative value disables the check.
84 REAL(kind=dp), SAVE :: eps_check_diag = -1.0_dp
85
86 ! Constants for the diag_type above
87 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_scalapack = 101, &
88 fm_diag_type_elpa = 102, &
89 fm_diag_type_cusolver = 103, &
90 fm_diag_type_dlaf = 104
91#if defined(__CUSOLVERMP)
92 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_cusolver
93#elif defined(__ELPA)
94 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_elpa
95#elif defined(__DLAF)
96 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_dlaf
97#else
98 INTEGER, PARAMETER, PUBLIC :: fm_diag_type_default = fm_diag_type_scalapack
99#endif
100
101 ! Public subroutines
102 PUBLIC :: choose_eigv_solver, &
103 cp_fm_block_jacobi, &
104 cp_fm_power, &
105 cp_fm_syevd, &
106 cp_fm_syevx, &
107 cp_fm_geeig, &
108 cp_fm_geeig_canon, &
109 diag_init, &
110 diag_finalize
111
112CONTAINS
113
114! **************************************************************************************************
115!> \brief Setup the diagonalization library to be used
116!> \param diag_lib diag_library flag from GLOBAL section in input
117!> \param fallback_applied .TRUE. if support for the requested library was not compiled-in and fallback
118!> to ScaLAPACK was applied, .FALSE. otherwise.
119!> \param elpa_kernel integer that determines which ELPA kernel to use for diagonalization
120!> \param elpa_neigvec_min_input ...
121!> \param elpa_qr logical that determines if ELPA should try to use QR to accelerate the
122!> diagonalization procedure of suitably sized matrices
123!> \param elpa_print logical that determines if information about the ELPA diagonalization should
124!> be printed
125!> \param elpa_qr_unsafe logical that enables potentially unsafe ELPA options
126!> \param dlaf_neigvec_min_input ...
127!> \param eps_check_diag_input ...
128!> \par History
129!> - Add support for DLA-Future (05.09.2023, RMeli)
130!> \author MI 11.2013
131! **************************************************************************************************
132 SUBROUTINE diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_neigvec_min_input, elpa_qr, &
133 elpa_print, elpa_qr_unsafe, dlaf_neigvec_min_input, eps_check_diag_input)
134 CHARACTER(LEN=*), INTENT(IN) :: diag_lib
135 LOGICAL, INTENT(OUT) :: fallback_applied
136 INTEGER, INTENT(IN) :: elpa_kernel, elpa_neigvec_min_input
137 LOGICAL, INTENT(IN) :: elpa_qr, elpa_print, elpa_qr_unsafe
138 INTEGER, INTENT(IN) :: dlaf_neigvec_min_input
139 REAL(kind=dp), INTENT(IN) :: eps_check_diag_input
140
141 fallback_applied = .false.
142
143 IF (diag_lib == "ScaLAPACK") THEN
144 diag_type = fm_diag_type_scalapack
145 ELSE IF (diag_lib == "ELPA") THEN
146#if defined (__ELPA)
147 ! ELPA is requested and available
148 diag_type = fm_diag_type_elpa
149#else
150 ! ELPA library requested but not linked, switch back to SL
151 diag_type = fm_diag_type_scalapack
152 fallback_applied = .true.
153#endif
154 ELSE IF (diag_lib == "cuSOLVER") THEN
155 diag_type = fm_diag_type_cusolver
156 ELSE IF (diag_lib == "DLAF") THEN
157#if defined (__DLAF)
158 diag_type = fm_diag_type_dlaf
159#else
160 cpabort("ERROR in diag_init: CP2K was not compiled with DLA-Future support")
161#endif
162 ELSE
163 cpabort("ERROR in diag_init: Initialization of unknown diagonalization library requested")
164 END IF
165
166 ! Initialization of requested diagonalization library:
167 IF (diag_type == fm_diag_type_elpa) THEN
168 CALL initialize_elpa_library()
169 CALL set_elpa_kernel(elpa_kernel)
170 CALL set_elpa_qr(elpa_qr, elpa_qr_unsafe)
171 CALL set_elpa_print(elpa_print)
172 END IF
173
174 elpa_neigvec_min = elpa_neigvec_min_input
175#if defined(__DLAF)
176 dlaf_neigvec_min = dlaf_neigvec_min_input
177#else
178 mark_used(dlaf_neigvec_min_input)
179#endif
180 eps_check_diag = eps_check_diag_input
181
182 END SUBROUTINE diag_init
183
184! **************************************************************************************************
185!> \brief Finalize the diagonalization library
186! **************************************************************************************************
187 SUBROUTINE diag_finalize()
188#if defined (__ELPA)
189 IF (diag_type == fm_diag_type_elpa) &
190 CALL finalize_elpa_library()
191#endif
192 END SUBROUTINE diag_finalize
193
194! **************************************************************************************************
195!> \brief Choose the Eigensolver depending on which library is available
196!> ELPA seems to be unstable for small systems
197!> \param matrix ...
198!> \param eigenvectors ...
199!> \param eigenvalues ...
200!> \param info ...
201!> \par info If present returns error code and prevents program stops.
202!> Works currently only for cp_fm_syevd with ScaLAPACK.
203!> Other solvers will end the program regardless of PRESENT(info).
204!> \par History
205!> - Do not use ELPA for small matrices and use instead ScaLAPACK as fallback (10.05.2021, MK)
206! **************************************************************************************************
207 SUBROUTINE choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
208
209 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
210 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
211 INTEGER, INTENT(OUT), OPTIONAL :: info
212
213 CHARACTER(LEN=*), PARAMETER :: routinen = 'choose_eigv_solver'
214
215 ! Sample peak memory
216 CALL m_memory()
217
218 IF (PRESENT(info)) info = 0 ! Default for solvers that do not return an info.
219
220 IF (diag_type == fm_diag_type_scalapack) THEN
221 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
222
223 ELSE IF (diag_type == fm_diag_type_elpa) THEN
224 IF (matrix%matrix_struct%nrow_global < elpa_neigvec_min) THEN
225 ! We don't trust ELPA with very small matrices.
226 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
227 ELSE
228 CALL cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
229 END IF
230
231 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
232 IF (matrix%matrix_struct%nrow_global < 64) THEN
233 ! We don't trust cuSolver with very small matrices.
234 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
235 ELSE
236 CALL cp_fm_diag_cusolver(matrix, eigenvectors, eigenvalues)
237 END IF
238
239#if defined(__DLAF)
240 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
241 IF (matrix%matrix_struct%nrow_global < dlaf_neigvec_min) THEN
242 ! Use ScaLAPACK for small matrices
243 CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
244 ELSE
245 CALL cp_fm_diag_dlaf(matrix, eigenvectors, eigenvalues)
246 END IF
247#endif
248
249 ELSE
250 cpabort("ERROR in "//routinen//": Invalid diagonalization type requested")
251 END IF
252
253 CALL check_diag(matrix, eigenvectors, nvec=SIZE(eigenvalues))
254
255 END SUBROUTINE choose_eigv_solver
256
257! **************************************************************************************************
258!> \brief Check result of diagonalization, i.e. the orthonormality of the eigenvectors
259!> \param matrix Work matrix
260!> \param eigenvectors Eigenvectors to be checked
261!> \param nvec ...
262! **************************************************************************************************
263 SUBROUTINE check_diag(matrix, eigenvectors, nvec)
264
265 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
266 INTEGER, INTENT(IN) :: nvec
267
268 CHARACTER(LEN=*), PARAMETER :: routinen = 'check_diag'
269
270 CHARACTER(LEN=default_string_length) :: diag_type_name
271 REAL(kind=dp) :: eps_abort, eps_warning
272 INTEGER :: handle, i, j, ncol, nrow, output_unit
273 LOGICAL :: check_eigenvectors
274#if defined(__SCALAPACK)
275 TYPE(cp_blacs_env_type), POINTER :: context
276 INTEGER :: il, jl, ipcol, iprow, &
277 mypcol, myprow, npcol, nprow
278 INTEGER, DIMENSION(9) :: desca
279#endif
280
281 CALL timeset(routinen, handle)
282
283 IF (diag_type == fm_diag_type_scalapack) THEN
284 diag_type_name = "SYEVD"
285 ELSE IF (diag_type == fm_diag_type_elpa) THEN
286 diag_type_name = "ELPA"
287 ELSE IF (diag_type == fm_diag_type_cusolver) THEN
288 diag_type_name = "CUSOLVER"
289 ELSE IF (diag_type == fm_diag_type_dlaf) THEN
290 diag_type_name = "DLAF"
291 ELSE
292 cpabort("Unknown diag_type")
293 END IF
294
295 output_unit = default_output_unit
296
297 eps_warning = eps_check_diag_default
298#if defined(__CHECK_DIAG)
299 check_eigenvectors = .true.
300 IF (eps_check_diag >= 0.0_dp) THEN
301 eps_warning = eps_check_diag
302 END IF
303#else
304 IF (eps_check_diag >= 0.0_dp) THEN
305 check_eigenvectors = .true.
306 eps_warning = eps_check_diag
307 ELSE
308 check_eigenvectors = .false.
309 END IF
310#endif
311 eps_abort = 10.0_dp*eps_warning
312
313 IF (check_eigenvectors) THEN
314#if defined(__SCALAPACK)
315 nrow = eigenvectors%matrix_struct%nrow_global
316 ncol = min(eigenvectors%matrix_struct%ncol_global, nvec)
317 CALL cp_fm_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, eigenvectors, eigenvectors, 0.0_dp, matrix)
318 context => matrix%matrix_struct%context
319 myprow = context%mepos(1)
320 mypcol = context%mepos(2)
321 nprow = context%num_pe(1)
322 npcol = context%num_pe(2)
323 desca(:) = matrix%matrix_struct%descriptor(:)
324 DO i = 1, ncol
325 DO j = 1, ncol
326 CALL infog2l(i, j, desca, nprow, npcol, myprow, mypcol, il, jl, iprow, ipcol)
327 IF ((iprow == myprow) .AND. (ipcol == mypcol)) THEN
328 IF (i == j) THEN
329 IF (abs(matrix%local_data(il, jl) - 1.0_dp) > eps_warning) THEN
330 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,ES10.3)") &
331 "The eigenvectors returned by "//trim(diag_type_name)//" are not orthonormal", &
332 "Matrix element (", i, ", ", j, ") = ", matrix%local_data(il, jl), &
333 "The deviation from the expected value 1 is", abs(matrix%local_data(il, jl) - 1.0_dp)
334 IF (abs(matrix%local_data(il, jl) - 1.0_dp) > eps_abort) THEN
335 cpabort("ERROR in "//routinen//": Check of matrix diagonalization failed")
336 ELSE
337 cpwarn("Check of matrix diagonalization failed in routine "//routinen)
338 END IF
339 END IF
340 ELSE
341 IF (abs(matrix%local_data(il, jl)) > eps_warning) THEN
342 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,ES10.3)") &
343 "The eigenvectors returned by "//trim(diag_type_name)//" are not orthonormal", &
344 "Matrix element (", i, ", ", j, ") = ", matrix%local_data(il, jl), &
345 "The deviation from the expected value 0 is", abs(matrix%local_data(il, jl))
346 IF (abs(matrix%local_data(il, jl)) > eps_abort) THEN
347 cpabort("ERROR in "//routinen//": Check of matrix diagonalization failed")
348 ELSE
349 cpwarn("Check of matrix diagonalization failed in routine "//routinen)
350 END IF
351 END IF
352 END IF
353 END IF
354 END DO
355 END DO
356#else
357 nrow = SIZE(eigenvectors%local_data, 1)
358 ncol = min(SIZE(eigenvectors%local_data, 2), nvec)
359 CALL dgemm("T", "N", ncol, ncol, nrow, 1.0_dp, &
360 eigenvectors%local_data(1, 1), nrow, &
361 eigenvectors%local_data(1, 1), nrow, &
362 0.0_dp, matrix%local_data(1, 1), nrow)
363 DO i = 1, ncol
364 DO j = 1, ncol
365 IF (i == j) THEN
366 IF (abs(matrix%local_data(i, j) - 1.0_dp) > eps_warning) THEN
367 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,ES10.3)") &
368 "The eigenvectors returned by "//trim(diag_type_name)//" are not orthonormal", &
369 "Matrix element (", i, ", ", j, ") = ", matrix%local_data(i, j), &
370 "The deviation from the expected value 1 is", abs(matrix%local_data(i, j) - 1.0_dp)
371 IF (abs(matrix%local_data(i, j) - 1.0_dp) > eps_abort) THEN
372 cpabort("ERROR in "//routinen//": Check of matrix diagonalization failed")
373 ELSE
374 cpwarn("Check of matrix diagonalization failed in routine "//routinen)
375 END IF
376 END IF
377 ELSE
378 IF (abs(matrix%local_data(i, j)) > eps_warning) THEN
379 WRITE (unit=output_unit, fmt="(/,T2,A,/,T2,A,I0,A,I0,A,F0.15,/,T2,A,ES10.3)") &
380 "The eigenvectors returned by "//trim(diag_type_name)//" are not orthonormal", &
381 "Matrix element (", i, ", ", j, ") = ", matrix%local_data(i, j), &
382 "The deviation from the expected value 0 is", abs(matrix%local_data(i, j))
383 IF (abs(matrix%local_data(i, j)) > eps_abort) THEN
384 cpabort("ERROR in "//routinen//": Check of matrix diagonalization failed")
385 ELSE
386 cpwarn("Check of matrix diagonalization failed in routine "//routinen)
387 END IF
388 END IF
389 END IF
390 END DO
391 END DO
392#endif
393 END IF
394
395 CALL timestop(handle)
396
397 END SUBROUTINE check_diag
398
399! **************************************************************************************************
400!> \brief Computes all eigenvalues and vectors of a real symmetric matrix
401!> significantly faster than syevx, scales also much better.
402!> Needs workspace to allocate all the eigenvectors
403!> \param matrix ...
404!> \param eigenvectors ...
405!> \param eigenvalues ...
406!> \param info ...
407!> \par matrix is supposed to be in upper triangular form, and overwritten by this routine
408!> \par info If present returns error code and prevents program stops.
409!> Works currently only for scalapack.
410!> Other solvers will end the program regardless of PRESENT(info).
411! **************************************************************************************************
412 SUBROUTINE cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
413
414 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
415 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
416 INTEGER, INTENT(OUT), OPTIONAL :: info
417
418 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_syevd'
419
420 INTEGER :: handle, myinfo, n, nmo
421 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig
422#if defined(__SCALAPACK)
423 TYPE(cp_fm_type) :: eigenvectors_new, matrix_new
424#else
425 CHARACTER(LEN=2*default_string_length) :: message
426 INTEGER :: liwork, lwork, nl
427 INTEGER, DIMENSION(:), POINTER :: iwork
428 REAL(kind=dp), DIMENSION(:, :), POINTER :: m
429 REAL(kind=dp), DIMENSION(:), POINTER :: work
430#endif
431
432 CALL timeset(routinen, handle)
433
434 myinfo = 0
435
436 n = matrix%matrix_struct%nrow_global
437 ALLOCATE (eig(n))
438
439#if defined(__SCALAPACK)
440 ! Determine if the input matrix needs to be redistributed before diagonalization.
441 ! Heuristics are used to determine the optimal number of CPUs for diagonalization.
442 ! The redistributed matrix is stored in matrix_new, which is just a pointer
443 ! to the original matrix if no redistribution is required
444 CALL cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new)
445
446 ! Call scalapack on CPUs that hold the new matrix
447 IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
448 IF (PRESENT(info)) THEN
449 CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig, myinfo)
450 ELSE
451 CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig)
452 END IF
453 END IF
454 ! Redistribute results and clean up
455 CALL cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
456#else
457 ! Retrieve the optimal work array sizes first
458 lwork = -1
459 liwork = -1
460 m => matrix%local_data
461 eig(:) = 0.0_dp
462
463 ALLOCATE (work(1))
464 work(:) = 0.0_dp
465 ALLOCATE (iwork(1))
466 iwork(:) = 0
467 nl = SIZE(m, 1)
468
469 CALL dsyevd('V', 'U', n, m(1, 1), nl, eig(1), work(1), lwork, iwork(1), liwork, myinfo)
470
471 IF (myinfo /= 0) THEN
472 WRITE (message, "(A,I0,A)") "ERROR in DSYEVD: Work space query failed (INFO = ", myinfo, ")"
473 IF (PRESENT(info)) THEN
474 cpwarn(trim(message))
475 ELSE
476 cpabort(trim(message))
477 END IF
478 END IF
479
480 ! Reallocate work arrays and perform diagonalisation
481 lwork = int(work(1))
482 DEALLOCATE (work)
483 ALLOCATE (work(lwork))
484 work(:) = 0.0_dp
485
486 liwork = iwork(1)
487 DEALLOCATE (iwork)
488 ALLOCATE (iwork(liwork))
489 iwork(:) = 0
490
491 CALL dsyevd('V', 'U', n, m(1, 1), nl, eig(1), work(1), lwork, iwork(1), liwork, myinfo)
492
493 IF (myinfo /= 0) THEN
494 WRITE (message, "(A,I0,A)") "ERROR in DSYEVD: Matrix diagonalization failed (INFO = ", myinfo, ")"
495 IF (PRESENT(info)) THEN
496 cpwarn(trim(message))
497 ELSE
498 cpabort(trim(message))
499 END IF
500 END IF
501
502 CALL cp_fm_to_fm(matrix, eigenvectors)
503
504 DEALLOCATE (iwork)
505 DEALLOCATE (work)
506#endif
507
508 IF (PRESENT(info)) info = myinfo
509
510 nmo = SIZE(eigenvalues, 1)
511 IF (nmo > n) THEN
512 eigenvalues(1:n) = eig(1:n)
513 ELSE
514 eigenvalues(1:nmo) = eig(1:nmo)
515 END IF
516
517 DEALLOCATE (eig)
518
519 CALL check_diag(matrix, eigenvectors, n)
520
521 CALL timestop(handle)
522
523 END SUBROUTINE cp_fm_syevd
524
525! **************************************************************************************************
526!> \brief ...
527!> \param matrix ...
528!> \param eigenvectors ...
529!> \param eigenvalues ...
530!> \param info ...
531! **************************************************************************************************
532 SUBROUTINE cp_fm_syevd_base(matrix, eigenvectors, eigenvalues, info)
533
534 TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
535 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
536 INTEGER, INTENT(OUT), OPTIONAL :: info
537
538 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_syevd_base'
539
540 CHARACTER(LEN=2*default_string_length) :: message
541 INTEGER :: handle, myinfo
542#if defined(__SCALAPACK)
543 TYPE(cp_blacs_env_type), POINTER :: context
544 INTEGER :: liwork, lwork, &
545 mypcol, myprow, n
546 INTEGER, DIMENSION(9) :: descm, descv
547 INTEGER, DIMENSION(:), POINTER :: iwork
548 REAL(kind=dp), DIMENSION(:), POINTER :: work
549 REAL(kind=dp), DIMENSION(:, :), POINTER :: m, v
550#if defined (__HAS_IEEE_EXCEPTIONS)
551 LOGICAL, DIMENSION(5) :: halt
552#endif
553#endif
554
555 CALL timeset(routinen, handle)
556
557 myinfo = 0
558
559#if defined(__SCALAPACK)
560
561 n = matrix%matrix_struct%nrow_global
562 m => matrix%local_data
563 context => matrix%matrix_struct%context
564 myprow = context%mepos(1)
565 mypcol = context%mepos(2)
566 descm(:) = matrix%matrix_struct%descriptor(:)
567
568 v => eigenvectors%local_data
569 descv(:) = eigenvectors%matrix_struct%descriptor(:)
570
571 liwork = 7*n + 8*context%num_pe(2) + 2
572 ALLOCATE (iwork(liwork))
573
574 ! Work space query
575 lwork = -1
576 ALLOCATE (work(1))
577
578 CALL pdsyevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
579 work(1), lwork, iwork(1), liwork, myinfo)
580
581 IF (matrix%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
582 WRITE (message, "(A,I0,A)") "ERROR in PDSYEVD: Work space query failed (INFO = ", myinfo, ")"
583 IF (PRESENT(info)) THEN
584 cpwarn(trim(message))
585 ELSE
586 cpabort(trim(message))
587 END IF
588 END IF
589
590 ! look here for a PDORMTR warning :-)
591 ! this routine seems to need more workspace than reported
592 ! (? lapack bug ...)
593 ! arbitrary additional memory ... we give 100000 more words
594 ! (it seems to depend on the block size used)
595
596 lwork = nint(work(1) + 100000)
597 ! lwork = NINT(work(1))
598 DEALLOCATE (work)
599 ALLOCATE (work(lwork))
600
601 ! Scalapack takes advantage of IEEE754 exceptions for speedup.
602 ! Therefore, we disable floating point traps temporarily.
603#if defined (__HAS_IEEE_EXCEPTIONS)
604 CALL ieee_get_halting_mode(ieee_all, halt)
605 CALL ieee_set_halting_mode(ieee_all, .false.)
606#endif
607
608 CALL pdsyevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
609 work(1), lwork, iwork(1), liwork, myinfo)
610
611#if defined (__HAS_IEEE_EXCEPTIONS)
612 CALL ieee_set_halting_mode(ieee_all, halt)
613#endif
614 IF (matrix%matrix_struct%para_env%is_source() .AND. (myinfo /= 0)) THEN
615 WRITE (message, "(A,I0,A)") "ERROR in PDSYEVD: Matrix diagonalization failed (INFO = ", myinfo, ")"
616 IF (PRESENT(info)) THEN
617 cpwarn(trim(message))
618 ELSE
619 cpabort(trim(message))
620 END IF
621 END IF
622
623 IF (PRESENT(info)) info = myinfo
624
625 DEALLOCATE (work)
626 DEALLOCATE (iwork)
627#else
628 mark_used(matrix)
629 mark_used(eigenvectors)
630 mark_used(eigenvalues)
631 myinfo = -1
632 IF (PRESENT(info)) info = myinfo
633 message = "ERROR in "//trim(routinen)// &
634 ": Matrix diagonalization using PDSYEVD requested without ScaLAPACK support"
635 cpabort(trim(message))
636#endif
637
638 CALL timestop(handle)
639
640 END SUBROUTINE cp_fm_syevd_base
641
642! **************************************************************************************************
643!> \brief compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack.
644!> If eigenvectors are required this routine will replicate a full matrix on each CPU...
645!> if more than a handful of vectors are needed, use cp_fm_syevd instead
646!> \param matrix ...
647!> \param eigenvectors ...
648!> \param eigenvalues ...
649!> \param neig ...
650!> \param work_syevx ...
651!> \par matrix is supposed to be in upper triangular form, and overwritten by this routine
652!> neig is the number of vectors needed (default all)
653!> work_syevx evec calculation only, is the fraction of the working buffer allowed (1.0 use full buffer)
654!> reducing this saves time, but might cause the routine to fail
655! **************************************************************************************************
656 SUBROUTINE cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
657
658 ! Diagonalise the symmetric n by n matrix using the LAPACK library.
659
660 TYPE(cp_fm_type), INTENT(IN) :: matrix
661 TYPE(cp_fm_type), OPTIONAL, INTENT(IN) :: eigenvectors
662 REAL(kind=dp), OPTIONAL, INTENT(IN) :: work_syevx
663 INTEGER, INTENT(IN), OPTIONAL :: neig
664 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
665
666 CHARACTER(LEN=*), PARAMETER :: routinen = "cp_fm_syevx"
667
668#if defined(__SCALAPACK)
669 REAL(kind=dp), PARAMETER :: orfac = -1.0_dp
670#endif
671 REAL(kind=dp), PARAMETER :: vl = 0.0_dp, &
672 vu = 0.0_dp
673
674 TYPE(cp_blacs_env_type), POINTER :: context
675 TYPE(cp_logger_type), POINTER :: logger
676 CHARACTER(LEN=1) :: job_type
677 REAL(kind=dp) :: abstol, work_syevx_local
678 INTEGER :: handle, info, &
679 liwork, lwork, m, mypcol, &
680 myprow, n, nb, npcol, &
681 nprow, output_unit, &
682 neig_local
683 LOGICAL :: ionode, needs_evecs
684 INTEGER, DIMENSION(:), ALLOCATABLE :: ifail, iwork
685 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: w, work
686 REAL(kind=dp), DIMENSION(:, :), POINTER :: a, z
687
688 REAL(kind=dp), EXTERNAL :: dlamch
689
690#if defined(__SCALAPACK)
691 INTEGER :: nn, np0, npe, nq0, nz
692 INTEGER, DIMENSION(9) :: desca, descz
693 INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr
694 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: gap
695 INTEGER, EXTERNAL :: iceil, numroc
696#if defined (__HAS_IEEE_EXCEPTIONS)
697 LOGICAL, DIMENSION(5) :: halt
698#endif
699#else
700 INTEGER :: nla, nlz
701 INTEGER, EXTERNAL :: ilaenv
702#endif
703
704 ! by default all
705 n = matrix%matrix_struct%nrow_global
706 neig_local = n
707 IF (PRESENT(neig)) neig_local = neig
708 IF (neig_local == 0) RETURN
709
710 CALL timeset(routinen, handle)
711
712 needs_evecs = PRESENT(eigenvectors)
713
714 logger => cp_get_default_logger()
715 ionode = logger%para_env%is_source()
716 n = matrix%matrix_struct%nrow_global
717
718 ! by default allocate all needed space
719 work_syevx_local = 1.0_dp
720 IF (PRESENT(work_syevx)) work_syevx_local = work_syevx
721
722 ! set scalapack job type
723 IF (needs_evecs) THEN
724 job_type = "V"
725 ELSE
726 job_type = "N"
727 END IF
728
729 ! target the most accurate calculation of the eigenvalues
730 abstol = 2.0_dp*dlamch("S")
731
732 context => matrix%matrix_struct%context
733 myprow = context%mepos(1)
734 mypcol = context%mepos(2)
735 nprow = context%num_pe(1)
736 npcol = context%num_pe(2)
737
738 ALLOCATE (w(n))
739 eigenvalues(:) = 0.0_dp
740#if defined(__SCALAPACK)
741
742 IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block) THEN
743 cpabort("ERROR in "//routinen//": Invalid blocksize (no square blocks) found")
744 END IF
745
746 a => matrix%local_data
747 desca(:) = matrix%matrix_struct%descriptor(:)
748
749 IF (needs_evecs) THEN
750 z => eigenvectors%local_data
751 descz(:) = eigenvectors%matrix_struct%descriptor(:)
752 ELSE
753 ! z will not be referenced
754 z => matrix%local_data
755 descz = desca
756 END IF
757
758 ! Get the optimal work storage size
759
760 npe = nprow*npcol
761 nb = matrix%matrix_struct%nrow_block
762 nn = max(n, nb, 2)
763 np0 = numroc(nn, nb, 0, 0, nprow)
764 nq0 = max(numroc(nn, nb, 0, 0, npcol), nb)
765
766 IF (needs_evecs) THEN
767 lwork = 5*n + max(5*nn, np0*nq0) + iceil(neig_local, npe)*nn + 2*nb*nb + &
768 int(work_syevx_local*real((neig_local - 1)*n, dp)) !!!! allocates a full matrix on every CPU !!!!!
769 ELSE
770 lwork = 5*n + max(5*nn, nb*(np0 + 1))
771 END IF
772 liwork = 6*max(n, npe + 1, 4)
773
774 ALLOCATE (gap(npe))
775 gap = 0.0_dp
776 ALLOCATE (iclustr(2*npe))
777 iclustr = 0
778 ALLOCATE (ifail(n))
779 ifail = 0
780 ALLOCATE (iwork(liwork))
781 ALLOCATE (work(lwork))
782
783 ! Scalapack takes advantage of IEEE754 exceptions for speedup.
784 ! Therefore, we disable floating point traps temporarily.
785#if defined (__HAS_IEEE_EXCEPTIONS)
786 CALL ieee_get_halting_mode(ieee_all, halt)
787 CALL ieee_set_halting_mode(ieee_all, .false.)
788#endif
789
790 CALL pdsyevx(job_type, "I", "U", n, a(1, 1), 1, 1, desca, vl, vu, 1, neig_local, abstol, m, nz, w(1), orfac, &
791 z(1, 1), 1, 1, descz, work(1), lwork, iwork(1), liwork, ifail(1), iclustr(1), gap, info)
792
793#if defined (__HAS_IEEE_EXCEPTIONS)
794 CALL ieee_set_halting_mode(ieee_all, halt)
795#endif
796
797 ! Error handling
798
799 IF (info /= 0) THEN
800 IF (ionode) THEN
801 output_unit = cp_logger_get_unit_nr(logger, local=.false.)
802 WRITE (unit=output_unit, fmt="(/,(T3,A,T12,1X,I10))") &
803 "info = ", info, &
804 "lwork = ", lwork, &
805 "liwork = ", liwork, &
806 "nz = ", nz
807 IF (info > 0) THEN
808 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
809 "ifail = ", ifail
810 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
811 "iclustr = ", iclustr
812 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,E10.3)))") &
813 "gap = ", gap
814 END IF
815 END IF
816 cpabort("ERROR in PDSYEVX (ScaLAPACK)")
817 END IF
818
819 ! Release work storage
820
821 DEALLOCATE (gap)
822 DEALLOCATE (iclustr)
823 DEALLOCATE (ifail)
824 DEALLOCATE (iwork)
825 DEALLOCATE (work)
826
827#else
828
829 a => matrix%local_data
830 IF (needs_evecs) THEN
831 z => eigenvectors%local_data
832 ELSE
833 ! z will not be referenced
834 z => matrix%local_data
835 END IF
836
837 ! Get the optimal work storage size
838
839 nb = max(ilaenv(1, "DSYTRD", "U", n, -1, -1, -1), &
840 ilaenv(1, "DORMTR", "U", n, -1, -1, -1))
841
842 lwork = max((nb + 3)*n, 8*n) + n ! sun bug fix
843 liwork = 5*n
844
845 ALLOCATE (ifail(n))
846 ifail = 0
847 ALLOCATE (iwork(liwork))
848 ALLOCATE (work(lwork))
849 info = 0
850 nla = SIZE(a, 1)
851 nlz = SIZE(z, 1)
852
853 CALL dsyevx(job_type, "I", "U", n, a(1, 1), nla, vl, vu, 1, neig_local, &
854 abstol, m, w, z(1, 1), nlz, work(1), lwork, iwork(1), ifail(1), info)
855
856 ! Error handling
857
858 IF (info /= 0) THEN
859 output_unit = cp_logger_get_unit_nr(logger, local=.false.)
860 WRITE (unit=output_unit, fmt="(/,(T3,A,T12,1X,I10))") &
861 "info = ", info
862 IF (info > 0) THEN
863 WRITE (unit=output_unit, fmt="(/,T3,A,(T12,6(1X,I10)))") &
864 "ifail = ", ifail
865 END IF
866 cpabort("Error in DSYEVX (ScaLAPACK)")
867 END IF
868
869 ! Release work storage
870
871 DEALLOCATE (ifail)
872 DEALLOCATE (iwork)
873 DEALLOCATE (work)
874
875#endif
876 eigenvalues(1:neig_local) = w(1:neig_local)
877 DEALLOCATE (w)
878
879 IF (needs_evecs) CALL check_diag(matrix, eigenvectors, neig_local)
880
881 CALL timestop(handle)
882
883 END SUBROUTINE cp_fm_syevx
884
885! **************************************************************************************************
886!> \brief ...
887!> \param matrix ...
888!> \param work ...
889!> \param exponent ...
890!> \param threshold ...
891!> \param n_dependent ...
892!> \param verbose ...
893!> \param eigvals ...
894! **************************************************************************************************
895 SUBROUTINE cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
896
897 ! Raise the real symmetric n by n matrix to the power given by
898 ! the exponent. All eigenvectors with a corresponding eigenvalue lower
899 ! than threshold are quenched. result in matrix
900
901 ! - Creation (29.03.1999, Matthias Krack)
902 ! - Parallelised using BLACS and ScaLAPACK (06.06.2001,MK)
903
904 TYPE(cp_fm_type), INTENT(IN) :: matrix, work
905 REAL(kind=dp), INTENT(IN) :: exponent, threshold
906 INTEGER, INTENT(OUT) :: n_dependent
907 LOGICAL, INTENT(IN), OPTIONAL :: verbose
908 REAL(kind=dp), DIMENSION(2), INTENT(OUT), &
909 OPTIONAL :: eigvals
910
911 CHARACTER(LEN=*), PARAMETER :: routinen = 'cp_fm_power'
912
913 INTEGER :: handle, icol_global, &
914 mypcol, myprow, &
915 ncol_global, npcol, nprow, &
916 nrow_global
917 LOGICAL :: my_verbose
918 REAL(kind=dp) :: condition_number, f, p
919 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues
920 REAL(kind=dp), DIMENSION(:, :), POINTER :: eigenvectors
921 TYPE(cp_blacs_env_type), POINTER :: context
922
923#if defined(__SCALAPACK)
924 INTEGER :: icol_local, ipcol, iprow, irow_global, &
925 irow_local, ncol_block, nrow_block
926 INTEGER, EXTERNAL :: indxg2l, indxg2p
927#endif
928
929 CALL timeset(routinen, handle)
930
931 my_verbose = .false.
932 IF (PRESENT(verbose)) my_verbose = verbose
933
934 context => matrix%matrix_struct%context
935 myprow = context%mepos(1)
936 mypcol = context%mepos(2)
937 nprow = context%num_pe(1)
938 npcol = context%num_pe(2)
939 n_dependent = 0
940 p = 0.5_dp*exponent
941
942 nrow_global = matrix%matrix_struct%nrow_global
943 ncol_global = matrix%matrix_struct%ncol_global
944
945 ALLOCATE (eigenvalues(ncol_global))
946
947 eigenvalues(:) = 0.0_dp
948
949 ! Compute the eigenvectors and eigenvalues
950
951 CALL choose_eigv_solver(matrix, work, eigenvalues)
952
953 IF (PRESENT(eigvals)) THEN
954 eigvals(1) = eigenvalues(1)
955 eigvals(2) = eigenvalues(ncol_global)
956 END IF
957
958#if defined(__SCALAPACK)
959 nrow_block = work%matrix_struct%nrow_block
960 ncol_block = work%matrix_struct%ncol_block
961
962 eigenvectors => work%local_data
963
964 ! Build matrix**exponent with eigenvector quenching
965
966 DO icol_global = 1, ncol_global
967
968 IF (eigenvalues(icol_global) < threshold) THEN
969
970 n_dependent = n_dependent + 1
971
972 ipcol = indxg2p(icol_global, ncol_block, mypcol, &
973 work%matrix_struct%first_p_pos(2), npcol)
974
975 IF (mypcol == ipcol) THEN
976 icol_local = indxg2l(icol_global, ncol_block, mypcol, &
977 work%matrix_struct%first_p_pos(2), npcol)
978 DO irow_global = 1, nrow_global
979 iprow = indxg2p(irow_global, nrow_block, myprow, &
980 work%matrix_struct%first_p_pos(1), nprow)
981 IF (myprow == iprow) THEN
982 irow_local = indxg2l(irow_global, nrow_block, myprow, &
983 work%matrix_struct%first_p_pos(1), nprow)
984 eigenvectors(irow_local, icol_local) = 0.0_dp
985 END IF
986 END DO
987 END IF
988
989 ELSE
990
991 f = eigenvalues(icol_global)**p
992
993 ipcol = indxg2p(icol_global, ncol_block, mypcol, &
994 work%matrix_struct%first_p_pos(2), npcol)
995
996 IF (mypcol == ipcol) THEN
997 icol_local = indxg2l(icol_global, ncol_block, mypcol, &
998 work%matrix_struct%first_p_pos(2), npcol)
999 DO irow_global = 1, nrow_global
1000 iprow = indxg2p(irow_global, nrow_block, myprow, &
1001 work%matrix_struct%first_p_pos(1), nprow)
1002 IF (myprow == iprow) THEN
1003 irow_local = indxg2l(irow_global, nrow_block, myprow, &
1004 work%matrix_struct%first_p_pos(1), nprow)
1005 eigenvectors(irow_local, icol_local) = &
1006 f*eigenvectors(irow_local, icol_local)
1007 END IF
1008 END DO
1009 END IF
1010
1011 END IF
1012
1013 END DO
1014
1015#else
1016
1017 eigenvectors => work%local_data
1018
1019 ! Build matrix**exponent with eigenvector quenching
1020
1021 DO icol_global = 1, ncol_global
1022
1023 IF (eigenvalues(icol_global) < threshold) THEN
1024
1025 n_dependent = n_dependent + 1
1026 eigenvectors(1:nrow_global, icol_global) = 0.0_dp
1027
1028 ELSE
1029
1030 f = eigenvalues(icol_global)**p
1031 eigenvectors(1:nrow_global, icol_global) = &
1032 f*eigenvectors(1:nrow_global, icol_global)
1033
1034 END IF
1035
1036 END DO
1037
1038#endif
1039 CALL cp_fm_syrk("U", "N", ncol_global, 1.0_dp, work, 1, 1, 0.0_dp, matrix)
1040 CALL cp_fm_upper_to_full(matrix, work)
1041
1042 ! Print some warnings/notes
1043 IF (matrix%matrix_struct%para_env%is_source() .AND. my_verbose) THEN
1044 condition_number = abs(eigenvalues(ncol_global)/eigenvalues(1))
1045 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,(T2,A,ES15.6))") &
1046 "CP_FM_POWER: smallest eigenvalue:", eigenvalues(1), &
1047 "CP_FM_POWER: largest eigenvalue: ", eigenvalues(ncol_global), &
1048 "CP_FM_POWER: condition number: ", condition_number
1049 IF (eigenvalues(1) <= 0.0_dp) THEN
1050 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
1051 "WARNING: matrix has a negative eigenvalue, tighten EPS_DEFAULT"
1052 END IF
1053 IF (condition_number > 1.0e12_dp) THEN
1054 WRITE (unit=cp_logger_get_default_unit_nr(), fmt="(/,T2,A)") &
1055 "WARNING: high condition number => possibly ill-conditioned matrix"
1056 END IF
1057 END IF
1058
1059 DEALLOCATE (eigenvalues)
1060
1061 CALL timestop(handle)
1062
1063 END SUBROUTINE cp_fm_power
1064
1065! **************************************************************************************************
1066!> \brief ...
1067!> \param matrix ...
1068!> \param eigenvectors ...
1069!> \param eigval ...
1070!> \param thresh ...
1071!> \param start_sec_block ...
1072! **************************************************************************************************
1073 SUBROUTINE cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, &
1074 start_sec_block)
1075
1076 ! Calculates block diagonalization of a full symmetric matrix
1077 ! It has its origin in cp_fm_syevx. This routine rotates only elements
1078 ! which are larger than a threshold values "thresh".
1079 ! start_sec_block is the start of the second block.
1080 ! IT DOES ONLY ONE SWEEP!
1081
1082 ! - Creation (07.10.2002, Martin Fengler)
1083 ! - Cosmetics (05.04.06, MK)
1084
1085 TYPE(cp_fm_type), INTENT(IN) :: eigenvectors, matrix
1086 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigval
1087 INTEGER, INTENT(IN) :: start_sec_block
1088 REAL(kind=dp), INTENT(IN) :: thresh
1089
1090 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_block_jacobi'
1091
1092 INTEGER :: handle
1093 REAL(kind=dp), DIMENSION(:, :), POINTER :: a, ev
1094
1095 REAL(kind=dp) :: tan_theta, tau, c, s
1096 INTEGER :: q, p, n
1097 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: c_ip
1098
1099#if defined(__SCALAPACK)
1100 TYPE(cp_blacs_env_type), POINTER :: context
1101
1102 INTEGER :: myprow, mypcol, nprow, npcol, block_dim_row, block_dim_col, &
1103 info, ev_row_block_size, iam, mynumrows, mype, npe, &
1104 q_loc
1105 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: a_loc, ev_loc
1106 INTEGER, DIMENSION(9) :: desca, descz, desc_a_block, &
1107 desc_ev_loc
1108 TYPE(mp_comm_type):: allgrp
1109 TYPE(cp_blacs_type) :: ictxt_loc
1110 INTEGER, EXTERNAL :: numroc, indxl2g, indxg2l
1111#endif
1112
1113 ! -------------------------------------------------------------------------
1114
1115 CALL timeset(routinen, handle)
1116
1117#if defined(__SCALAPACK)
1118 context => matrix%matrix_struct%context
1119 allgrp = matrix%matrix_struct%para_env
1120
1121 myprow = context%mepos(1)
1122 mypcol = context%mepos(2)
1123 nprow = context%num_pe(1)
1124 npcol = context%num_pe(2)
1125
1126 n = matrix%matrix_struct%nrow_global
1127
1128 a => matrix%local_data
1129 desca(:) = matrix%matrix_struct%descriptor(:)
1130 ev => eigenvectors%local_data
1131 descz(:) = eigenvectors%matrix_struct%descriptor(:)
1132
1133 ! Copy the block to be rotated to the master process firstly and broadcast to all processes
1134 ! start_sec_block defines where the second block starts!
1135 ! Block will be processed together with the OO block
1136
1137 block_dim_row = start_sec_block - 1
1138 block_dim_col = n - block_dim_row
1139 ALLOCATE (a_loc(block_dim_row, block_dim_col))
1140
1141 mype = matrix%matrix_struct%para_env%mepos
1142 npe = matrix%matrix_struct%para_env%num_pe
1143 ! Get a new context
1144 CALL ictxt_loc%gridinit(matrix%matrix_struct%para_env, 'Row-major', nprow*npcol, 1)
1145
1146 CALL descinit(desc_a_block, block_dim_row, block_dim_col, block_dim_row, &
1147 block_dim_col, 0, 0, ictxt_loc%get_handle(), block_dim_row, info)
1148
1149 CALL pdgemr2d(block_dim_row, block_dim_col, a, 1, start_sec_block, desca, &
1150 a_loc, 1, 1, desc_a_block, context%get_handle())
1151 ! Only the master (root) process received data yet
1152 CALL allgrp%bcast(a_loc, 0)
1153
1154 ! Since each process owns now the upper block, the eigenvectors can be re-sorted in such a way that
1155 ! each process has a NN*1 grid, i.e. the process owns a bunch of rows which can be modified locally
1156
1157 ! Initialize distribution of the eigenvectors
1158 iam = mype
1159 ev_row_block_size = n/(nprow*npcol)
1160 mynumrows = numroc(n, ev_row_block_size, iam, 0, nprow*npcol)
1161
1162 ALLOCATE (ev_loc(mynumrows, n), c_ip(mynumrows))
1163
1164 CALL descinit(desc_ev_loc, n, n, ev_row_block_size, n, 0, 0, ictxt_loc%get_handle(), &
1165 mynumrows, info)
1166
1167 CALL pdgemr2d(n, n, ev, 1, 1, descz, ev_loc, 1, 1, desc_ev_loc, context%get_handle())
1168
1169 ! Start block diagonalization of matrix
1170
1171 q_loc = 0
1172
1173 DO q = start_sec_block, n
1174 q_loc = q_loc + 1
1175 DO p = 1, (start_sec_block - 1)
1176
1177 IF (abs(a_loc(p, q_loc)) > thresh) THEN
1178
1179 tau = (eigval(q) - eigval(p))/(2.0_dp*a_loc(p, q_loc))
1180
1181 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1182
1183 ! Cos(theta)
1184 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
1185 s = tan_theta*c
1186
1187 ! Calculate eigenvectors: Q*J
1188 ! c_ip = c*EV_loc(:,p) - s*EV_loc(:,q)
1189 ! c_iq = s*EV_loc(:,p) + c*EV_loc(:,q)
1190 ! EV(:,p) = c_ip
1191 ! EV(:,q) = c_iq
1192 CALL dcopy(mynumrows, ev_loc(1, p), 1, c_ip(1), 1)
1193 CALL dscal(mynumrows, c, ev_loc(1, p), 1)
1194 CALL daxpy(mynumrows, -s, ev_loc(1, q), 1, ev_loc(1, p), 1)
1195 CALL dscal(mynumrows, c, ev_loc(1, q), 1)
1196 CALL daxpy(mynumrows, s, c_ip(1), 1, ev_loc(1, q), 1)
1197
1198 END IF
1199
1200 END DO
1201 END DO
1202
1203 ! Copy eigenvectors back to the original distribution
1204 CALL pdgemr2d(n, n, ev_loc, 1, 1, desc_ev_loc, ev, 1, 1, descz, context%get_handle())
1205
1206 ! Release work storage
1207 DEALLOCATE (a_loc, ev_loc, c_ip)
1208
1209 CALL ictxt_loc%gridexit()
1210
1211#else
1212
1213 n = matrix%matrix_struct%nrow_global
1214
1215 ALLOCATE (c_ip(n)) ! Local eigenvalue vector
1216
1217 a => matrix%local_data ! Contains the Matrix to be worked on
1218 ev => eigenvectors%local_data ! Contains the eigenvectors up to blocksize, rest is garbage
1219
1220 ! Start matrix diagonalization
1221
1222 tan_theta = 0.0_dp
1223 tau = 0.0_dp
1224
1225 DO q = start_sec_block, n
1226 DO p = 1, (start_sec_block - 1)
1227
1228 IF (abs(a(p, q)) > thresh) THEN
1229
1230 tau = (eigval(q) - eigval(p))/(2.0_dp*a(p, q))
1231
1232 tan_theta = sign(1.0_dp, tau)/(abs(tau) + sqrt(1.0_dp + tau*tau))
1233
1234 ! Cos(theta)
1235 c = 1.0_dp/sqrt(1.0_dp + tan_theta*tan_theta)
1236 s = tan_theta*c
1237
1238 ! Calculate eigenvectors: Q*J
1239 ! c_ip = c*EV(:,p) - s*EV(:,q)
1240 ! c_iq = s*EV(:,p) + c*EV(:,q)
1241 ! EV(:,p) = c_ip
1242 ! EV(:,q) = c_iq
1243 CALL dcopy(n, ev(1, p), 1, c_ip(1), 1)
1244 CALL dscal(n, c, ev(1, p), 1)
1245 CALL daxpy(n, -s, ev(1, q), 1, ev(1, p), 1)
1246 CALL dscal(n, c, ev(1, q), 1)
1247 CALL daxpy(n, s, c_ip(1), 1, ev(1, q), 1)
1248
1249 END IF
1250
1251 END DO
1252 END DO
1253
1254 ! Release work storage
1255
1256 DEALLOCATE (c_ip)
1257
1258#endif
1259
1260 CALL timestop(handle)
1261
1262 END SUBROUTINE cp_fm_block_jacobi
1263
1264! **************************************************************************************************
1265!> \brief General Eigenvalue Problem AX = BXE
1266!> Single option version: Cholesky decomposition of B
1267!> \param amatrix ...
1268!> \param bmatrix ...
1269!> \param eigenvectors ...
1270!> \param eigenvalues ...
1271!> \param work ...
1272! **************************************************************************************************
1273 SUBROUTINE cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
1274
1275 TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1276 REAL(kind=dp), DIMENSION(:) :: eigenvalues
1277 TYPE(cp_fm_type), INTENT(IN) :: work
1278
1279 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_geeig'
1280
1281 INTEGER :: handle, nao, nmo
1282
1283 CALL timeset(routinen, handle)
1284
1285 CALL cp_fm_get_info(amatrix, nrow_global=nao)
1286 nmo = SIZE(eigenvalues)
1287 ! Cholesky decompose S=U(T)U
1288 CALL cp_fm_cholesky_decompose(bmatrix)
1289 ! Invert to get U^(-1)
1290 CALL cp_fm_triangular_invert(bmatrix)
1291 ! Reduce to get U^(-T) * H * U^(-1)
1292 CALL cp_fm_triangular_multiply(bmatrix, amatrix, side="R")
1293 CALL cp_fm_triangular_multiply(bmatrix, amatrix, transpose_tr=.true.)
1294 ! Diagonalize
1295 CALL choose_eigv_solver(matrix=amatrix, eigenvectors=work, &
1296 eigenvalues=eigenvalues)
1297 ! Restore vectors C = U^(-1) * C*
1298 CALL cp_fm_triangular_multiply(bmatrix, work)
1299 CALL cp_fm_to_fm(work, eigenvectors, nmo)
1300
1301 CALL timestop(handle)
1302
1303 END SUBROUTINE cp_fm_geeig
1304
1305! **************************************************************************************************
1306!> \brief General Eigenvalue Problem AX = BXE
1307!> Use canonical diagonalization : U*s**(-1/2)
1308!> \param amatrix ...
1309!> \param bmatrix ...
1310!> \param eigenvectors ...
1311!> \param eigenvalues ...
1312!> \param work ...
1313!> \param epseig ...
1314! **************************************************************************************************
1315 SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
1316
1317 TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1318 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1319 TYPE(cp_fm_type), INTENT(IN) :: work
1320 REAL(kind=dp), INTENT(IN) :: epseig
1321
1322 CHARACTER(len=*), PARAMETER :: routinen = 'cp_fm_geeig_canon'
1323
1324 INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
1325 nmo, nx
1326 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: seigval
1327
1328 CALL timeset(routinen, handle)
1329
1330 ! Test sizees
1331 CALL cp_fm_get_info(amatrix, nrow_global=nao)
1332 nmo = SIZE(eigenvalues)
1333 ALLOCATE (seigval(nao))
1334
1335 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
1336 CALL cp_fm_scale(-1.0_dp, bmatrix)
1337 CALL choose_eigv_solver(matrix=bmatrix, eigenvectors=work, eigenvalues=seigval)
1338 seigval(:) = -seigval(:)
1339 nc = nao
1340 DO i = 1, nao
1341 IF (seigval(i) < epseig) THEN
1342 nc = i - 1
1343 EXIT
1344 END IF
1345 END DO
1346 cpassert(nc /= 0)
1347
1348 IF (nc /= nao) THEN
1349 IF (nc < nmo) THEN
1350 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
1351 ncol = nmo - nc
1352 CALL cp_fm_to_fm(work, eigenvectors, ncol, nc + 1, nc + 1)
1353 END IF
1354 ! Set NULL space in eigenvector matrix of S to zero
1355 DO icol = nc + 1, nao
1356 DO irow = 1, nao
1357 CALL cp_fm_set_element(work, irow, icol, 0.0_dp)
1358 END DO
1359 END DO
1360 ! Set small eigenvalues to a dummy save value
1361 seigval(nc + 1:nao) = 1.0_dp
1362 END IF
1363 ! Calculate U*s**(-1/2)
1364 seigval(:) = 1.0_dp/sqrt(seigval(:))
1365 CALL cp_fm_column_scale(work, seigval)
1366 ! Reduce to get U^(-T) * H * U^(-1)
1367 CALL cp_fm_gemm("T", "N", nao, nao, nao, 1.0_dp, work, amatrix, 0.0_dp, bmatrix)
1368 CALL cp_fm_gemm("N", "N", nao, nao, nao, 1.0_dp, bmatrix, work, 0.0_dp, amatrix)
1369 IF (nc /= nao) THEN
1370 ! set diagonal values to save large value
1371 DO icol = nc + 1, nao
1372 CALL cp_fm_set_element(amatrix, icol, icol, 10000.0_dp)
1373 END DO
1374 END IF
1375 ! Diagonalize
1376 CALL choose_eigv_solver(matrix=amatrix, eigenvectors=bmatrix, eigenvalues=eigenvalues)
1377 nx = min(nc, nmo)
1378 ! Restore vectors C = U^(-1) * C*
1379 CALL cp_fm_gemm("N", "N", nao, nx, nc, 1.0_dp, work, bmatrix, 0.0_dp, eigenvectors)
1380
1381 DEALLOCATE (seigval)
1382
1383 CALL timestop(handle)
1384
1385 END SUBROUTINE cp_fm_geeig_canon
1386
1387END MODULE cp_fm_diag
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:195
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_types
wrappers for the actual blacs calls. all functionality needed in the code should actually be provide ...
Definition cp_blacs_types.F:18
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_gemm
subroutine, public cp_fm_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)
computes matrix_c = beta * matrix_c + alpha * ( matrix_a ** transa ) * ( matrix_b ** transb )
Definition cp_fm_basic_linalg.F:440
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_syrk
subroutine, public cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
performs a rank-k update of a symmetric matrix_c matrix_c = beta * matrix_c + alpha * matrix_a * tran...
Definition cp_fm_basic_linalg.F:675
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1759
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1680
cp_fm_basic_linalg::cp_fm_triangular_invert
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
Definition cp_fm_basic_linalg.F:2193
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_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_fm_basic_linalg.F:1592
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:50
cp_fm_cusolver_api
Wrapper for cuSOLVERMp.
Definition cp_fm_cusolver_api.F:12
cp_fm_cusolver_api::cp_fm_diag_cusolver
subroutine, public cp_fm_diag_cusolver(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the cuSOLVERMp library.
Definition cp_fm_cusolver_api.F:36
cp_fm_diag_utils
Auxiliary tools to redistribute cp_fm_type matrices before and after diagonalization....
Definition cp_fm_diag_utils.F:16
cp_fm_diag_utils::cp_fm_redistribute_end
subroutine, public cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
Redistributes eigenvectors and eigenvalues back to the original communicator group.
Definition cp_fm_diag_utils.F:456
cp_fm_diag_utils::cp_fm_redistribute_start
subroutine, public cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, caller_is_elpa, redist_info)
Determines the optimal number of CPUs for matrix diagonalization and redistributes the input matrices...
Definition cp_fm_diag_utils.F:264
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::cp_fm_block_jacobi
subroutine, public cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, start_sec_block)
...
Definition cp_fm_diag.F:1075
cp_fm_diag::eps_check_diag_default
real(kind=dp), parameter, public eps_check_diag_default
Definition cp_fm_diag.F:69
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:896
cp_fm_diag::fm_diag_type_cusolver
integer, parameter, public fm_diag_type_cusolver
Definition cp_fm_diag.F:87
cp_fm_diag::fm_diag_type_dlaf
integer, parameter, public fm_diag_type_dlaf
Definition cp_fm_diag.F:87
cp_fm_diag::fm_diag_type_scalapack
integer, parameter, public fm_diag_type_scalapack
Definition cp_fm_diag.F:87
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition cp_fm_diag.F:413
cp_fm_diag::diag_finalize
subroutine, public diag_finalize()
Finalize the diagonalization library.
Definition cp_fm_diag.F:188
cp_fm_diag::diag_init
subroutine, public diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_neigvec_min_input, elpa_qr, elpa_print, elpa_qr_unsafe, dlaf_neigvec_min_input, eps_check_diag_input)
Setup the diagonalization library to be used.
Definition cp_fm_diag.F:134
cp_fm_diag::fm_diag_type_default
integer, parameter, public fm_diag_type_default
Definition cp_fm_diag.F:98
cp_fm_diag::cp_fm_geeig
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
Definition cp_fm_diag.F:1274
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:208
cp_fm_diag::fm_diag_type_elpa
integer, parameter, public fm_diag_type_elpa
Definition cp_fm_diag.F:87
cp_fm_diag::cp_fm_geeig_canon
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2)
Definition cp_fm_diag.F:1316
cp_fm_diag::cp_fm_syevx
subroutine, public cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)
compute eigenvalues and optionally eigenvectors of a real symmetric matrix using scalapack....
Definition cp_fm_diag.F:657
cp_fm_dlaf_api
Definition cp_fm_dlaf_api.F:8
cp_fm_dlaf_api::cp_fm_diag_dlaf
subroutine, public cp_fm_diag_dlaf(matrix, eigenvectors, eigenvalues)
...
Definition cp_fm_dlaf_api.F:145
cp_fm_elpa
Wrapper for ELPA.
Definition cp_fm_elpa.F:12
cp_fm_elpa::set_elpa_print
subroutine, public set_elpa_print(flag)
Sets a flag that determines if additional information about the ELPA diagonalization should be printe...
Definition cp_fm_elpa.F:266
cp_fm_elpa::set_elpa_kernel
subroutine, public set_elpa_kernel(requested_kernel)
Sets the active ELPA kernel.
Definition cp_fm_elpa.F:201
cp_fm_elpa::cp_fm_diag_elpa
subroutine, public cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the ELPA library.
Definition cp_fm_elpa.F:278
cp_fm_elpa::set_elpa_qr
subroutine, public set_elpa_qr(use_qr, use_qr_unsafe)
Sets flags that determines if ELPA should try to use QR during diagonalization If use_qr = ....
Definition cp_fm_elpa.F:254
cp_fm_elpa::initialize_elpa_library
subroutine, public initialize_elpa_library()
Initialize the ELPA library.
Definition cp_fm_elpa.F:177
cp_fm_elpa::finalize_elpa_library
subroutine, public finalize_elpa_library()
Finalize the ELPA library.
Definition cp_fm_elpa.F:189
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:700
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_logger_get_unit_nr
integer function, public cp_logger_get_unit_nr(logger, local)
returns the unit nr for the requested kind of log.
Definition cp_log_handling.F:498
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::default_output_unit
integer, parameter, public default_output_unit
Definition machine.F:45
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:332
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_blacs_types::cp_blacs_type
Definition cp_blacs_types.F:34
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:189