(git:b195825)
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 ! **************************************************************************************************
17 MODULE 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 
112 CONTAINS
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 
1387 END 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_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:374
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:1816
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:609
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:1693
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:1614
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:2127
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:1526
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