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