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 (__parallel)
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 9009 : 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 9009 : fallback_applied = .FALSE.
142 :
143 9009 : IF (diag_lib == "ScaLAPACK") THEN
144 188 : diag_type = FM_DIAG_TYPE_SCALAPACK
145 8821 : ELSE IF (diag_lib == "ELPA") THEN
146 : #if defined (__ELPA)
147 : ! ELPA is requested and available
148 8821 : 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 9009 : IF (diag_type == FM_DIAG_TYPE_ELPA) THEN
168 8821 : CALL initialize_elpa_library()
169 8821 : CALL set_elpa_kernel(elpa_kernel)
170 8821 : CALL set_elpa_qr(elpa_qr, elpa_qr_unsafe)
171 8821 : CALL set_elpa_print(elpa_print)
172 : END IF
173 :
174 9009 : 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 9009 : eps_check_diag = eps_check_diag_input
181 :
182 9009 : END SUBROUTINE diag_init
183 :
184 : ! **************************************************************************************************
185 : !> \brief Finalize the diagonalization library
186 : ! **************************************************************************************************
187 8799 : SUBROUTINE diag_finalize()
188 : #if defined (__ELPA)
189 8799 : IF (diag_type == FM_DIAG_TYPE_ELPA) &
190 8611 : CALL finalize_elpa_library()
191 : #endif
192 8799 : 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 160625 : 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 160625 : CALL m_memory()
217 :
218 160625 : IF (PRESENT(info)) info = 0 ! Default for solvers that do not return an info.
219 :
220 160625 : IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN
221 4114 : CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
222 :
223 156511 : ELSE IF (diag_type == FM_DIAG_TYPE_ELPA) THEN
224 156511 : IF (matrix%matrix_struct%nrow_global < elpa_neigvec_min) THEN
225 : ! We don't trust ELPA with very small matrices.
226 97571 : CALL cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
227 : ELSE
228 58940 : 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 160625 : CALL check_diag(matrix, eigenvectors, nvec=SIZE(eigenvalues))
254 :
255 160625 : 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 264908 : 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(__parallel)
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 264908 : CALL timeset(routineN, handle)
282 :
283 264908 : IF (diag_type == FM_DIAG_TYPE_SCALAPACK) THEN
284 8228 : diag_type_name = "SYEVD"
285 256680 : ELSE IF (diag_type == FM_DIAG_TYPE_ELPA) THEN
286 256680 : 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 264908 : output_unit = default_output_unit
296 :
297 264908 : 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 264908 : 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 264908 : eps_abort = 10.0_dp*eps_warning
312 :
313 264908 : IF (check_eigenvectors) THEN
314 : #if defined(__parallel)
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 264908 : CALL timestop(handle)
396 :
397 264908 : 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 104243 : 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(__parallel)
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 104243 : CALL timeset(routineN, handle)
433 :
434 104243 : myinfo = 0
435 :
436 104243 : n = matrix%matrix_struct%nrow_global
437 312729 : ALLOCATE (eig(n))
438 :
439 : #if defined(__parallel)
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 104243 : CALL cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new)
445 :
446 : ! Call scalapack on CPUs that hold the new matrix
447 104243 : IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
448 52779 : IF (PRESENT(info)) THEN
449 1109 : CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig, myinfo)
450 : ELSE
451 51670 : CALL cp_fm_syevd_base(matrix_new, eigenvectors_new, eig)
452 : END IF
453 : END IF
454 : ! Redistribute results and clean up
455 104243 : 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 104243 : IF (PRESENT(info)) info = myinfo
509 :
510 104243 : nmo = SIZE(eigenvalues, 1)
511 104243 : IF (nmo > n) THEN
512 4744 : eigenvalues(1:n) = eig(1:n)
513 : ELSE
514 722756 : eigenvalues(1:nmo) = eig(1:nmo)
515 : END IF
516 :
517 104243 : DEALLOCATE (eig)
518 :
519 104243 : CALL check_diag(matrix, eigenvectors, n)
520 :
521 104243 : CALL timestop(handle)
522 :
523 208486 : END SUBROUTINE cp_fm_syevd
524 :
525 : ! **************************************************************************************************
526 : !> \brief ...
527 : !> \param matrix ...
528 : !> \param eigenvectors ...
529 : !> \param eigenvalues ...
530 : !> \param info ...
531 : ! **************************************************************************************************
532 52779 : 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(__parallel)
543 : TYPE(cp_blacs_env_type), POINTER :: context
544 : INTEGER :: liwork, lwork, &
545 : mypcol, myprow, n
546 : INTEGER, DIMENSION(9) :: descm, descv
547 52779 : INTEGER, DIMENSION(:), POINTER :: iwork
548 52779 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
549 52779 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: m, v
550 : #if defined (__HAS_IEEE_EXCEPTIONS)
551 : LOGICAL, DIMENSION(5) :: halt
552 : #endif
553 : #endif
554 :
555 52779 : CALL timeset(routineN, handle)
556 :
557 52779 : myinfo = 0
558 :
559 : #if defined(__parallel)
560 :
561 52779 : n = matrix%matrix_struct%nrow_global
562 52779 : m => matrix%local_data
563 52779 : context => matrix%matrix_struct%context
564 52779 : myprow = context%mepos(1)
565 52779 : mypcol = context%mepos(2)
566 527790 : descm(:) = matrix%matrix_struct%descriptor(:)
567 :
568 52779 : v => eigenvectors%local_data
569 527790 : descv(:) = eigenvectors%matrix_struct%descriptor(:)
570 :
571 52779 : liwork = 7*n + 8*context%num_pe(2) + 2
572 158337 : ALLOCATE (iwork(liwork))
573 :
574 : ! Work space query
575 52779 : lwork = -1
576 52779 : 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 52779 : work(1), lwork, iwork(1), liwork, myinfo)
580 :
581 52779 : 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 52779 : lwork = NINT(work(1) + 100000)
597 : ! lwork = NINT(work(1))
598 52779 : DEALLOCATE (work)
599 158337 : 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 52779 : 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 52779 : 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 52779 : IF (PRESENT(info)) info = myinfo
624 :
625 52779 : DEALLOCATE (work)
626 52779 : 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 52779 : CALL timestop(handle)
639 :
640 52779 : 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(__parallel)
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(__parallel)
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(__parallel)
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 1092 : 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 1092 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: eigenvectors
921 : TYPE(cp_blacs_env_type), POINTER :: context
922 :
923 : #if defined(__parallel)
924 : INTEGER :: icol_local, ipcol, iprow, irow_global, &
925 : irow_local
926 : #endif
927 :
928 1092 : CALL timeset(routineN, handle)
929 :
930 1092 : my_verbose = .FALSE.
931 1092 : IF (PRESENT(verbose)) my_verbose = verbose
932 :
933 1092 : context => matrix%matrix_struct%context
934 1092 : myprow = context%mepos(1)
935 1092 : mypcol = context%mepos(2)
936 1092 : nprow = context%num_pe(1)
937 1092 : npcol = context%num_pe(2)
938 1092 : n_dependent = 0
939 1092 : p = 0.5_dp*exponent
940 :
941 1092 : nrow_global = matrix%matrix_struct%nrow_global
942 1092 : ncol_global = matrix%matrix_struct%ncol_global
943 :
944 3276 : ALLOCATE (eigenvalues(ncol_global))
945 :
946 23755 : eigenvalues(:) = 0.0_dp
947 :
948 : ! Compute the eigenvectors and eigenvalues
949 :
950 1092 : CALL choose_eigv_solver(matrix, work, eigenvalues)
951 :
952 1092 : IF (PRESENT(eigvals)) THEN
953 330 : eigvals(1) = eigenvalues(1)
954 330 : eigvals(2) = eigenvalues(ncol_global)
955 : END IF
956 :
957 : #if defined(__parallel)
958 1092 : eigenvectors => work%local_data
959 :
960 : ! Build matrix**exponent with eigenvector quenching
961 :
962 23755 : DO icol_global = 1, ncol_global
963 :
964 23755 : IF (eigenvalues(icol_global) < threshold) THEN
965 :
966 82 : n_dependent = n_dependent + 1
967 :
968 82 : ipcol = work%matrix_struct%g2p_col(icol_global)
969 :
970 82 : IF (mypcol == ipcol) THEN
971 82 : icol_local = work%matrix_struct%g2l_col(icol_global)
972 4714 : DO irow_global = 1, nrow_global
973 4632 : iprow = work%matrix_struct%g2p_row(irow_global)
974 4714 : IF (myprow == iprow) THEN
975 2316 : irow_local = work%matrix_struct%g2l_row(irow_global)
976 2316 : eigenvectors(irow_local, icol_local) = 0.0_dp
977 : END IF
978 : END DO
979 : END IF
980 :
981 : ELSE
982 :
983 22581 : f = eigenvalues(icol_global)**p
984 :
985 22581 : ipcol = work%matrix_struct%g2p_col(icol_global)
986 :
987 22581 : IF (mypcol == ipcol) THEN
988 22581 : icol_local = work%matrix_struct%g2l_col(icol_global)
989 1554160 : DO irow_global = 1, nrow_global
990 1531579 : iprow = work%matrix_struct%g2p_row(irow_global)
991 1554160 : IF (myprow == iprow) THEN
992 809864 : irow_local = work%matrix_struct%g2l_row(irow_global)
993 : eigenvectors(irow_local, icol_local) = &
994 809864 : f*eigenvectors(irow_local, icol_local)
995 : END IF
996 : END DO
997 : END IF
998 :
999 : END IF
1000 :
1001 : END DO
1002 :
1003 : #else
1004 :
1005 : eigenvectors => work%local_data
1006 :
1007 : ! Build matrix**exponent with eigenvector quenching
1008 :
1009 : DO icol_global = 1, ncol_global
1010 :
1011 : IF (eigenvalues(icol_global) < threshold) THEN
1012 :
1013 : n_dependent = n_dependent + 1
1014 : eigenvectors(1:nrow_global, icol_global) = 0.0_dp
1015 :
1016 : ELSE
1017 :
1018 : f = eigenvalues(icol_global)**p
1019 : eigenvectors(1:nrow_global, icol_global) = &
1020 : f*eigenvectors(1:nrow_global, icol_global)
1021 :
1022 : END IF
1023 :
1024 : END DO
1025 :
1026 : #endif
1027 1092 : CALL cp_fm_syrk("U", "N", ncol_global, 1.0_dp, work, 1, 1, 0.0_dp, matrix)
1028 1092 : CALL cp_fm_upper_to_full(matrix, work)
1029 :
1030 : ! Print some warnings/notes
1031 1092 : IF (matrix%matrix_struct%para_env%is_source() .AND. my_verbose) THEN
1032 0 : condition_number = ABS(eigenvalues(ncol_global)/eigenvalues(1))
1033 : WRITE (UNIT=cp_logger_get_default_unit_nr(), FMT="(/,(T2,A,ES15.6))") &
1034 0 : "CP_FM_POWER: smallest eigenvalue:", eigenvalues(1), &
1035 0 : "CP_FM_POWER: largest eigenvalue: ", eigenvalues(ncol_global), &
1036 0 : "CP_FM_POWER: condition number: ", condition_number
1037 0 : IF (eigenvalues(1) <= 0.0_dp) THEN
1038 : WRITE (UNIT=cp_logger_get_default_unit_nr(), FMT="(/,T2,A)") &
1039 0 : "WARNING: matrix has a negative eigenvalue, tighten EPS_DEFAULT"
1040 : END IF
1041 0 : IF (condition_number > 1.0E12_dp) THEN
1042 : WRITE (UNIT=cp_logger_get_default_unit_nr(), FMT="(/,T2,A)") &
1043 0 : "WARNING: high condition number => possibly ill-conditioned matrix"
1044 : END IF
1045 : END IF
1046 :
1047 1092 : DEALLOCATE (eigenvalues)
1048 :
1049 1092 : CALL timestop(handle)
1050 :
1051 1092 : END SUBROUTINE cp_fm_power
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief ...
1055 : !> \param matrix ...
1056 : !> \param eigenvectors ...
1057 : !> \param eigval ...
1058 : !> \param thresh ...
1059 : !> \param start_sec_block ...
1060 : ! **************************************************************************************************
1061 18 : SUBROUTINE cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, &
1062 : start_sec_block)
1063 :
1064 : ! Calculates block diagonalization of a full symmetric matrix
1065 : ! It has its origin in cp_fm_syevx. This routine rotates only elements
1066 : ! which are larger than a threshold values "thresh".
1067 : ! start_sec_block is the start of the second block.
1068 : ! IT DOES ONLY ONE SWEEP!
1069 :
1070 : ! - Creation (07.10.2002, Martin Fengler)
1071 : ! - Cosmetics (05.04.06, MK)
1072 :
1073 : TYPE(cp_fm_type), INTENT(IN) :: eigenvectors, matrix
1074 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigval
1075 : INTEGER, INTENT(IN) :: start_sec_block
1076 : REAL(KIND=dp), INTENT(IN) :: thresh
1077 :
1078 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_block_jacobi'
1079 :
1080 : INTEGER :: handle
1081 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, ev
1082 :
1083 : REAL(KIND=dp) :: tan_theta, tau, c, s
1084 : INTEGER :: q, p, N
1085 18 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: c_ip
1086 :
1087 : #if defined(__parallel)
1088 : TYPE(cp_blacs_env_type), POINTER :: context
1089 :
1090 : INTEGER :: myprow, mypcol, nprow, npcol, block_dim_row, block_dim_col, &
1091 : info, ev_row_block_size, iam, mynumrows, mype, npe, &
1092 : q_loc
1093 18 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: a_loc, ev_loc
1094 : INTEGER, DIMENSION(9) :: desca, descz, desc_a_block, &
1095 : desc_ev_loc
1096 : TYPE(mp_comm_type):: allgrp
1097 : TYPE(cp_blacs_type) :: ictxt_loc
1098 : INTEGER, EXTERNAL :: numroc
1099 : #endif
1100 :
1101 : ! -------------------------------------------------------------------------
1102 :
1103 18 : CALL timeset(routineN, handle)
1104 :
1105 : #if defined(__parallel)
1106 18 : context => matrix%matrix_struct%context
1107 18 : allgrp = matrix%matrix_struct%para_env
1108 :
1109 18 : myprow = context%mepos(1)
1110 18 : mypcol = context%mepos(2)
1111 18 : nprow = context%num_pe(1)
1112 18 : npcol = context%num_pe(2)
1113 :
1114 18 : N = matrix%matrix_struct%nrow_global
1115 :
1116 18 : A => matrix%local_data
1117 180 : desca(:) = matrix%matrix_struct%descriptor(:)
1118 18 : EV => eigenvectors%local_data
1119 180 : descz(:) = eigenvectors%matrix_struct%descriptor(:)
1120 :
1121 : ! Copy the block to be rotated to the master process firstly and broadcast to all processes
1122 : ! start_sec_block defines where the second block starts!
1123 : ! Block will be processed together with the OO block
1124 :
1125 18 : block_dim_row = start_sec_block - 1
1126 18 : block_dim_col = N - block_dim_row
1127 72 : ALLOCATE (A_loc(block_dim_row, block_dim_col))
1128 :
1129 18 : mype = matrix%matrix_struct%para_env%mepos
1130 18 : npe = matrix%matrix_struct%para_env%num_pe
1131 : ! Get a new context
1132 18 : CALL ictxt_loc%gridinit(matrix%matrix_struct%para_env, 'Row-major', NPROW*NPCOL, 1)
1133 :
1134 : CALL descinit(desc_a_block, block_dim_row, block_dim_col, block_dim_row, &
1135 18 : block_dim_col, 0, 0, ictxt_loc%get_handle(), block_dim_row, info)
1136 :
1137 : CALL pdgemr2d(block_dim_row, block_dim_col, A, 1, start_sec_block, desca, &
1138 18 : A_loc, 1, 1, desc_a_block, context%get_handle())
1139 : ! Only the master (root) process received data yet
1140 18 : CALL allgrp%bcast(A_loc, 0)
1141 :
1142 : ! Since each process owns now the upper block, the eigenvectors can be re-sorted in such a way that
1143 : ! each process has a NN*1 grid, i.e. the process owns a bunch of rows which can be modified locally
1144 :
1145 : ! Initialize distribution of the eigenvectors
1146 18 : iam = mype
1147 18 : ev_row_block_size = n/(nprow*npcol)
1148 18 : mynumrows = NUMROC(N, ev_row_block_size, iam, 0, NPROW*NPCOL)
1149 :
1150 108 : ALLOCATE (EV_loc(mynumrows, N), c_ip(mynumrows))
1151 :
1152 : CALL descinit(desc_ev_loc, N, N, ev_row_block_size, N, 0, 0, ictxt_loc%get_handle(), &
1153 18 : mynumrows, info)
1154 :
1155 18 : CALL pdgemr2d(N, N, EV, 1, 1, descz, EV_loc, 1, 1, desc_ev_loc, context%get_handle())
1156 :
1157 : ! Start block diagonalization of matrix
1158 :
1159 18 : q_loc = 0
1160 :
1161 1170 : DO q = start_sec_block, N
1162 1152 : q_loc = q_loc + 1
1163 148626 : DO p = 1, (start_sec_block - 1)
1164 :
1165 148608 : IF (ABS(A_loc(p, q_loc)) > thresh) THEN
1166 :
1167 117566 : tau = (eigval(q) - eigval(p))/(2.0_dp*A_loc(p, q_loc))
1168 :
1169 117566 : tan_theta = SIGN(1.0_dp, tau)/(ABS(tau) + SQRT(1.0_dp + tau*tau))
1170 :
1171 : ! Cos(theta)
1172 117566 : c = 1.0_dp/SQRT(1.0_dp + tan_theta*tan_theta)
1173 117566 : s = tan_theta*c
1174 :
1175 : ! Calculate eigenvectors: Q*J
1176 : ! c_ip = c*EV_loc(:,p) - s*EV_loc(:,q)
1177 : ! c_iq = s*EV_loc(:,p) + c*EV_loc(:,q)
1178 : ! EV(:,p) = c_ip
1179 : ! EV(:,q) = c_iq
1180 117566 : CALL dcopy(mynumrows, EV_loc(1, p), 1, c_ip(1), 1)
1181 117566 : CALL dscal(mynumrows, c, EV_loc(1, p), 1)
1182 117566 : CALL daxpy(mynumrows, -s, EV_loc(1, q), 1, EV_loc(1, p), 1)
1183 117566 : CALL dscal(mynumrows, c, EV_loc(1, q), 1)
1184 117566 : CALL daxpy(mynumrows, s, c_ip(1), 1, EV_loc(1, q), 1)
1185 :
1186 : END IF
1187 :
1188 : END DO
1189 : END DO
1190 :
1191 : ! Copy eigenvectors back to the original distribution
1192 18 : CALL pdgemr2d(N, N, EV_loc, 1, 1, desc_ev_loc, EV, 1, 1, descz, context%get_handle())
1193 :
1194 : ! Release work storage
1195 18 : DEALLOCATE (A_loc, EV_loc, c_ip)
1196 :
1197 18 : CALL ictxt_loc%gridexit()
1198 :
1199 : #else
1200 :
1201 : N = matrix%matrix_struct%nrow_global
1202 :
1203 : ALLOCATE (c_ip(N)) ! Local eigenvalue vector
1204 :
1205 : A => matrix%local_data ! Contains the Matrix to be worked on
1206 : EV => eigenvectors%local_data ! Contains the eigenvectors up to blocksize, rest is garbage
1207 :
1208 : ! Start matrix diagonalization
1209 :
1210 : tan_theta = 0.0_dp
1211 : tau = 0.0_dp
1212 :
1213 : DO q = start_sec_block, N
1214 : DO p = 1, (start_sec_block - 1)
1215 :
1216 : IF (ABS(A(p, q)) > thresh) THEN
1217 :
1218 : tau = (eigval(q) - eigval(p))/(2.0_dp*A(p, q))
1219 :
1220 : tan_theta = SIGN(1.0_dp, tau)/(ABS(tau) + SQRT(1.0_dp + tau*tau))
1221 :
1222 : ! Cos(theta)
1223 : c = 1.0_dp/SQRT(1.0_dp + tan_theta*tan_theta)
1224 : s = tan_theta*c
1225 :
1226 : ! Calculate eigenvectors: Q*J
1227 : ! c_ip = c*EV(:,p) - s*EV(:,q)
1228 : ! c_iq = s*EV(:,p) + c*EV(:,q)
1229 : ! EV(:,p) = c_ip
1230 : ! EV(:,q) = c_iq
1231 : CALL dcopy(N, EV(1, p), 1, c_ip(1), 1)
1232 : CALL dscal(N, c, EV(1, p), 1)
1233 : CALL daxpy(N, -s, EV(1, q), 1, EV(1, p), 1)
1234 : CALL dscal(N, c, EV(1, q), 1)
1235 : CALL daxpy(N, s, c_ip(1), 1, EV(1, q), 1)
1236 :
1237 : END IF
1238 :
1239 : END DO
1240 : END DO
1241 :
1242 : ! Release work storage
1243 :
1244 : DEALLOCATE (c_ip)
1245 :
1246 : #endif
1247 :
1248 18 : CALL timestop(handle)
1249 :
1250 90 : END SUBROUTINE cp_fm_block_jacobi
1251 :
1252 : ! **************************************************************************************************
1253 : !> \brief General Eigenvalue Problem AX = BXE
1254 : !> Single option version: Cholesky decomposition of B
1255 : !> \param amatrix ...
1256 : !> \param bmatrix ...
1257 : !> \param eigenvectors ...
1258 : !> \param eigenvalues ...
1259 : !> \param work ...
1260 : ! **************************************************************************************************
1261 262 : SUBROUTINE cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
1262 :
1263 : TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1264 : REAL(KIND=dp), DIMENSION(:) :: eigenvalues
1265 : TYPE(cp_fm_type), INTENT(IN) :: work
1266 :
1267 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geeig'
1268 :
1269 : INTEGER :: handle, nao, nmo
1270 :
1271 262 : CALL timeset(routineN, handle)
1272 :
1273 262 : CALL cp_fm_get_info(amatrix, nrow_global=nao)
1274 262 : nmo = SIZE(eigenvalues)
1275 : ! Cholesky decompose S=U(T)U
1276 262 : CALL cp_fm_cholesky_decompose(bmatrix)
1277 : ! Invert to get U^(-1)
1278 262 : CALL cp_fm_triangular_invert(bmatrix)
1279 : ! Reduce to get U^(-T) * H * U^(-1)
1280 262 : CALL cp_fm_triangular_multiply(bmatrix, amatrix, side="R")
1281 262 : CALL cp_fm_triangular_multiply(bmatrix, amatrix, transpose_tr=.TRUE.)
1282 : ! Diagonalize
1283 : CALL choose_eigv_solver(matrix=amatrix, eigenvectors=work, &
1284 262 : eigenvalues=eigenvalues)
1285 : ! Restore vectors C = U^(-1) * C*
1286 262 : CALL cp_fm_triangular_multiply(bmatrix, work)
1287 262 : CALL cp_fm_to_fm(work, eigenvectors, nmo)
1288 :
1289 262 : CALL timestop(handle)
1290 :
1291 262 : END SUBROUTINE cp_fm_geeig
1292 :
1293 : ! **************************************************************************************************
1294 : !> \brief General Eigenvalue Problem AX = BXE
1295 : !> Use canonical diagonalization : U*s**(-1/2)
1296 : !> \param amatrix ...
1297 : !> \param bmatrix ...
1298 : !> \param eigenvectors ...
1299 : !> \param eigenvalues ...
1300 : !> \param work ...
1301 : !> \param epseig ...
1302 : ! **************************************************************************************************
1303 64 : SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
1304 :
1305 : TYPE(cp_fm_type), INTENT(IN) :: amatrix, bmatrix, eigenvectors
1306 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1307 : TYPE(cp_fm_type), INTENT(IN) :: work
1308 : REAL(KIND=dp), INTENT(IN) :: epseig
1309 :
1310 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geeig_canon'
1311 :
1312 : INTEGER :: handle, i, icol, irow, nao, nc, ncol, &
1313 : nmo, nx
1314 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: seigval
1315 :
1316 64 : CALL timeset(routineN, handle)
1317 :
1318 : ! Test sizees
1319 64 : CALL cp_fm_get_info(amatrix, nrow_global=nao)
1320 64 : nmo = SIZE(eigenvalues)
1321 192 : ALLOCATE (seigval(nao))
1322 :
1323 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
1324 64 : CALL cp_fm_scale(-1.0_dp, bmatrix)
1325 64 : CALL choose_eigv_solver(matrix=bmatrix, eigenvectors=work, eigenvalues=seigval)
1326 4784 : seigval(:) = -seigval(:)
1327 64 : nc = nao
1328 4492 : DO i = 1, nao
1329 4492 : IF (seigval(i) < epseig) THEN
1330 40 : nc = i - 1
1331 40 : EXIT
1332 : END IF
1333 : END DO
1334 64 : CPASSERT(nc /= 0)
1335 :
1336 64 : IF (nc /= nao) THEN
1337 40 : IF (nc < nmo) THEN
1338 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
1339 0 : ncol = nmo - nc
1340 0 : CALL cp_fm_to_fm(work, eigenvectors, ncol, nc + 1, nc + 1)
1341 : END IF
1342 : ! Set NULL space in eigenvector matrix of S to zero
1343 332 : DO icol = nc + 1, nao
1344 36172 : DO irow = 1, nao
1345 36132 : CALL cp_fm_set_element(work, irow, icol, 0.0_dp)
1346 : END DO
1347 : END DO
1348 : ! Set small eigenvalues to a dummy save value
1349 332 : seigval(nc + 1:nao) = 1.0_dp
1350 : END IF
1351 : ! Calculate U*s**(-1/2)
1352 4784 : seigval(:) = 1.0_dp/SQRT(seigval(:))
1353 64 : CALL cp_fm_column_scale(work, seigval)
1354 : ! Reduce to get U^(-T) * H * U^(-1)
1355 64 : CALL cp_fm_gemm("T", "N", nao, nao, nao, 1.0_dp, work, amatrix, 0.0_dp, bmatrix)
1356 64 : CALL cp_fm_gemm("N", "N", nao, nao, nao, 1.0_dp, bmatrix, work, 0.0_dp, amatrix)
1357 64 : IF (nc /= nao) THEN
1358 : ! set diagonal values to save large value
1359 332 : DO icol = nc + 1, nao
1360 332 : CALL cp_fm_set_element(amatrix, icol, icol, 10000.0_dp)
1361 : END DO
1362 : END IF
1363 : ! Diagonalize
1364 64 : CALL choose_eigv_solver(matrix=amatrix, eigenvectors=bmatrix, eigenvalues=eigenvalues)
1365 64 : nx = MIN(nc, nmo)
1366 : ! Restore vectors C = U^(-1) * C*
1367 64 : CALL cp_fm_gemm("N", "N", nao, nx, nc, 1.0_dp, work, bmatrix, 0.0_dp, eigenvectors)
1368 :
1369 64 : DEALLOCATE (seigval)
1370 :
1371 64 : CALL timestop(handle)
1372 :
1373 64 : END SUBROUTINE cp_fm_geeig_canon
1374 :
1375 : END MODULE cp_fm_diag
|