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