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