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