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