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