Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Wrapper for ELPA (complex matrices, i.e. cp_cfm_type)
10 : ! **************************************************************************************************
11 : MODULE cp_cfm_elpa
12 : USE cp_blacs_env, ONLY: cp_blacs_env_create, &
13 : cp_blacs_env_release, &
14 : cp_blacs_env_type
15 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, &
16 : cp_cfm_uplo_to_full
17 : USE cp_cfm_types, ONLY: cp_cfm_create, &
18 : cp_cfm_release, &
19 : cp_cfm_to_cfm, &
20 : cp_cfm_type
21 : USE cp_fm_diag_utils, ONLY: cp_cfm_redistribute_start, &
22 : cp_cfm_redistribute_end, &
23 : cp_fm_redistribute_info
24 : USE cp_fm_elpa, ONLY: elpa_one_stage, &
25 : elpa_print
26 : USE cp_fm_struct, ONLY: cp_fm_struct_create, &
27 : cp_fm_struct_get, &
28 : cp_fm_struct_release, &
29 : cp_fm_struct_type
30 : USE cp_log_handling, ONLY: cp_get_default_logger, &
31 : cp_logger_get_default_io_unit, &
32 : cp_logger_type, &
33 : cp_to_string
34 : USE kinds, ONLY: default_string_length, dp
35 : USE machine, ONLY: m_cpuid_static, &
36 : MACHINE_CPU_GENERIC, &
37 : MACHINE_X86_SSE4, &
38 : MACHINE_X86_AVX, &
39 : MACHINE_X86_AVX2, &
40 : MACHINE_X86_AVX512
41 : USE message_passing, ONLY: mp_comm_self, &
42 : mp_comm_type, &
43 : mp_para_env_type
44 : USE OMP_LIB, ONLY: omp_get_max_threads
45 : USE parallel_rng_types, ONLY: rng_stream_type, &
46 : UNIFORM
47 : #if defined(__HAS_IEEE_EXCEPTIONS)
48 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
49 : ieee_set_halting_mode, &
50 : IEEE_ALL
51 : #endif
52 : #include "../base/base_uses.f90"
53 :
54 : #if defined(__ELPA)
55 : USE elpa_constants, ONLY: ELPA_SOLVER_1STAGE, ELPA_SOLVER_2STAGE, ELPA_OK, &
56 : ELPA_2STAGE_COMPLEX_INVALID, &
57 : ELPA_2STAGE_COMPLEX_DEFAULT, &
58 : ELPA_2STAGE_COMPLEX_GENERIC, &
59 : ELPA_2STAGE_COMPLEX_GENERIC_SIMPLE, &
60 : ELPA_2STAGE_COMPLEX_BGP, &
61 : ELPA_2STAGE_COMPLEX_BGQ, &
62 : ELPA_2STAGE_COMPLEX_SSE_BLOCK1, &
63 : ELPA_2STAGE_COMPLEX_SSE_BLOCK2, &
64 : ELPA_2STAGE_COMPLEX_AVX_BLOCK1, &
65 : ELPA_2STAGE_COMPLEX_AVX_BLOCK2, &
66 : ELPA_2STAGE_COMPLEX_AVX2_BLOCK1, &
67 : ELPA_2STAGE_COMPLEX_AVX2_BLOCK2, &
68 : ELPA_2STAGE_COMPLEX_AVX512_BLOCK1, &
69 : ELPA_2STAGE_COMPLEX_AVX512_BLOCK2, &
70 : ELPA_2STAGE_COMPLEX_NVIDIA_GPU, &
71 : ELPA_2STAGE_COMPLEX_AMD_GPU, &
72 : ELPA_2STAGE_COMPLEX_INTEL_GPU_SYCL
73 :
74 : USE elpa, ONLY: elpa_t, &
75 : elpa_allocate, elpa_deallocate
76 : #endif
77 :
78 : IMPLICIT NONE
79 :
80 : PRIVATE
81 :
82 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_elpa'
83 :
84 : #if defined(__ELPA)
85 : INTEGER, DIMENSION(16), PARAMETER :: elpa_c_kernel_ids = [ &
86 : ELPA_2STAGE_COMPLEX_INVALID, & ! auto
87 : ELPA_2STAGE_COMPLEX_GENERIC, &
88 : ELPA_2STAGE_COMPLEX_GENERIC_SIMPLE, &
89 : ELPA_2STAGE_COMPLEX_BGP, &
90 : ELPA_2STAGE_COMPLEX_BGQ, &
91 : ELPA_2STAGE_COMPLEX_SSE_BLOCK1, &
92 : ELPA_2STAGE_COMPLEX_SSE_BLOCK2, &
93 : ELPA_2STAGE_COMPLEX_AVX_BLOCK1, &
94 : ELPA_2STAGE_COMPLEX_AVX_BLOCK2, &
95 : ELPA_2STAGE_COMPLEX_AVX2_BLOCK1, &
96 : ELPA_2STAGE_COMPLEX_AVX2_BLOCK2, &
97 : ELPA_2STAGE_COMPLEX_AVX512_BLOCK1, &
98 : ELPA_2STAGE_COMPLEX_AVX512_BLOCK2, &
99 : ELPA_2STAGE_COMPLEX_NVIDIA_GPU, &
100 : ELPA_2STAGE_COMPLEX_AMD_GPU, &
101 : ELPA_2STAGE_COMPLEX_INTEL_GPU_SYCL]
102 :
103 : CHARACTER(len=14), DIMENSION(SIZE(elpa_c_kernel_ids)), PARAMETER :: &
104 : elpa_c_kernel_names = [CHARACTER(len=14) :: &
105 : "AUTO", &
106 : "GENERIC", &
107 : "GENERIC_SIMPLE", &
108 : "BGP", &
109 : "BGQ", &
110 : "SSE_BLOCK1", &
111 : "SSE_BLOCK2", &
112 : "AVX_BLOCK1", &
113 : "AVX_BLOCK2", &
114 : "AVX2_BLOCK1", &
115 : "AVX2_BLOCK2", &
116 : "AVX512_BLOCK1", &
117 : "AVX512_BLOCK2", &
118 : "NVIDIA_GPU", &
119 : "AMD_GPU", &
120 : "INTEL_GPU"]
121 :
122 : CHARACTER(len=44), DIMENSION(SIZE(elpa_c_kernel_ids)), PARAMETER :: &
123 : elpa_c_kernel_descriptions = [CHARACTER(len=44) :: &
124 : "Automatically selected kernel", &
125 : "Generic kernel", &
126 : "Simplified generic kernel", &
127 : "Kernel optimized for IBM BGP", &
128 : "Kernel optimized for IBM BGQ", &
129 : "Kernel optimized for x86_64/SSE (block=1)", &
130 : "Kernel optimized for x86_64/SSE (block=2)", &
131 : "Kernel optimized for Intel AVX (block=1)", &
132 : "Kernel optimized for Intel AVX (block=2)", &
133 : "Kernel optimized for Intel AVX2 (block=1)", &
134 : "Kernel optimized for Intel AVX2 (block=2)", &
135 : "Kernel optimized for Intel AVX-512 (block=1)", &
136 : "Kernel optimized for Intel AVX-512 (block=2)", &
137 : "Kernel targeting Nvidia GPUs", &
138 : "Kernel targeting AMD GPUs", &
139 : "Kernel targeting Intel GPUs"]
140 : #else
141 : INTEGER, DIMENSION(1), PARAMETER :: elpa_c_kernel_ids = [-1]
142 : CHARACTER(len=14), DIMENSION(1), PARAMETER :: elpa_c_kernel_names = ["AUTO"]
143 : CHARACTER(len=44), DIMENSION(1), PARAMETER :: elpa_c_kernel_descriptions = ["Automatically selected kernel"]
144 : #endif
145 :
146 : #if defined(__ELPA)
147 : INTEGER, SAVE :: elpa_c_kernel = elpa_c_kernel_ids(1) ! auto
148 : #endif
149 :
150 : ! Runtime correctness check state
151 : LOGICAL, SAVE, PRIVATE :: elpa_c_correctness_checked = .FALSE.
152 : LOGICAL, SAVE, PRIVATE :: elpa_c_broken = .FALSE.
153 :
154 : PUBLIC :: cp_cfm_diag_elpa, &
155 : set_elpa_c_kernel, &
156 : check_elpa_c_kernel_correctness, &
157 : is_elpa_c_broken, &
158 : elpa_c_kernel_ids, &
159 : elpa_c_kernel_names, &
160 : elpa_c_kernel_descriptions
161 :
162 : CONTAINS
163 :
164 : #if defined(__ELPA)
165 : ! **************************************************************************************************
166 : !> \brief Return a printable name for an ELPA complex kernel.
167 : !> \param kernel ELPA complex kernel id
168 : !> \return ...
169 : ! **************************************************************************************************
170 29 : FUNCTION get_elpa_c_kernel_name(kernel) RESULT(kernel_name)
171 : INTEGER, INTENT(IN) :: kernel
172 : CHARACTER(len=default_string_length) :: kernel_name
173 :
174 : INTEGER :: i
175 :
176 29 : kernel_name = "id: "//TRIM(ADJUSTL(cp_to_string(kernel)))
177 211 : DO i = 1, SIZE(elpa_c_kernel_ids)
178 211 : IF (elpa_c_kernel_ids(i) == kernel) THEN
179 29 : kernel_name = elpa_c_kernel_names(i)
180 29 : EXIT
181 : END IF
182 : END DO
183 29 : END FUNCTION get_elpa_c_kernel_name
184 : #endif
185 :
186 : ! **************************************************************************************************
187 : !> \brief Sets the active ELPA kernel for complex matrices.
188 : !> \param requested_kernel one of the elpa_c_kernel_ids
189 : ! **************************************************************************************************
190 10646 : SUBROUTINE set_elpa_c_kernel(requested_kernel)
191 : INTEGER, INTENT(IN) :: requested_kernel
192 :
193 : #if defined(__ELPA)
194 : INTEGER :: cpuid
195 :
196 10646 : elpa_c_kernel = requested_kernel
197 :
198 : ! Resolve AUTO kernel.
199 10646 : IF (elpa_c_kernel == ELPA_2STAGE_COMPLEX_INVALID) THEN
200 10640 : cpuid = m_cpuid_static()
201 0 : SELECT CASE (cpuid)
202 : CASE (MACHINE_CPU_GENERIC)
203 0 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_GENERIC
204 : CASE (MACHINE_X86_SSE4)
205 0 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_SSE_BLOCK2
206 : CASE (MACHINE_X86_AVX)
207 0 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_AVX_BLOCK2
208 : CASE (MACHINE_X86_AVX2)
209 10640 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_AVX2_BLOCK2
210 : CASE (MACHINE_X86_AVX512)
211 10640 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_AVX512_BLOCK2
212 : END SELECT
213 :
214 : ! Prefer GPU kernel if available.
215 : #if !defined(__NO_OFFLOAD_ELPA)
216 : #if defined(__OFFLOAD_CUDA)
217 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_NVIDIA_GPU
218 : #endif
219 : #if defined(__OFFLOAD_HIP)
220 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_AMD_GPU
221 : #endif
222 : #if defined(__OFFLOAD_OPENCL)
223 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_INTEL_GPU_SYCL
224 : #endif
225 : #endif
226 : ! If we could not find a suitable kernel then use ELPA_2STAGE_COMPLEX_DEFAULT.
227 10640 : IF (elpa_c_kernel == ELPA_2STAGE_COMPLEX_INVALID) THEN
228 0 : elpa_c_kernel = ELPA_2STAGE_COMPLEX_DEFAULT
229 : END IF
230 : END IF
231 : #else
232 : MARK_USED(requested_kernel)
233 : #endif
234 10646 : END SUBROUTINE set_elpa_c_kernel
235 :
236 : #if defined(__ELPA)
237 : ! **************************************************************************************************
238 : !> \brief Returns .TRUE. if the complex kernel is an SSE/AVX/AVX2/AVX512 BLOCK2 variant
239 : !> (these may be affected by the GCC 15.2 regression, marekandreas/elpa#77).
240 : !> \param kernel ...
241 : !> \return ...
242 : ! **************************************************************************************************
243 32 : FUNCTION is_block2_c_kernel(kernel) RESULT(is_block2)
244 : INTEGER, INTENT(IN) :: kernel
245 : LOGICAL :: is_block2
246 :
247 : is_block2 = (kernel == ELPA_2STAGE_COMPLEX_SSE_BLOCK2) .OR. &
248 : (kernel == ELPA_2STAGE_COMPLEX_AVX_BLOCK2) .OR. &
249 : (kernel == ELPA_2STAGE_COMPLEX_AVX2_BLOCK2) .OR. &
250 32 : (kernel == ELPA_2STAGE_COMPLEX_AVX512_BLOCK2)
251 32 : END FUNCTION is_block2_c_kernel
252 : #endif
253 :
254 : ! **************************************************************************************************
255 : !> \brief One-time runtime correctness check for ELPA complex BLOCK2 kernels.
256 : !> A small deterministic eigenproblem is solved on a single process with the
257 : !> same production routine used for the actual diagonalizations, and the
258 : !> eigenpair residual as well as the eigenvector orthogonality are verified.
259 : !> This guards against mis-compiled kernels (GCC 15.2 regression, marekandreas/elpa#77).
260 : !> If the check fails, a module flag is set such that cp_cfm_heevd falls back to ScaLAPACK.
261 : !> \param para_env communicator of the run, used to agree on the check result
262 : ! **************************************************************************************************
263 32 : SUBROUTINE check_elpa_c_kernel_correctness(para_env)
264 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
265 :
266 : #if defined(__ELPA)
267 : CHARACTER(LEN=*), PARAMETER :: routineN = 'check_elpa_c_kernel_correctness'
268 : CHARACTER(LEN=4*default_string_length) :: message
269 : INTEGER :: handle, i, is_broken, j
270 : ! Geometry aligned with the standalone reproducer of marekandreas/elpa#77.
271 : INTEGER, PARAMETER :: na_test = 20, nev_test = 18, nblk_test = 8
272 : REAL(KIND=dp) :: eps_ortho, eps_residual, rnd_im, rnd_re, test
273 : REAL(KIND=dp), PARAMETER :: th = 1.0E-8_dp
274 32 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues
275 : TYPE(cp_blacs_env_type), POINTER :: context
276 : TYPE(cp_cfm_type) :: eigenvectors, matrix, matrix_ref
277 : TYPE(cp_fm_redistribute_info) :: rdinfo
278 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
279 : TYPE(mp_para_env_type) :: para_env_self
280 : TYPE(rng_stream_type) :: rng_stream
281 :
282 32 : IF (elpa_c_correctness_checked) RETURN ! one-shot per process
283 :
284 32 : IF (.NOT. is_block2_c_kernel(elpa_c_kernel)) RETURN ! only BLOCK2 kernels are affected
285 30 : IF (elpa_one_stage) RETURN ! the 1-stage solver does not use kernels
286 24 : elpa_c_correctness_checked = .TRUE.
287 :
288 24 : CALL timeset(routineN, handle)
289 :
290 : ! Solve a small deterministic Hermitian eigenproblem redundantly on every process,
291 : ! reusing the production ELPA path (kernel setup, fallbacks and solver call).
292 24 : para_env_self = mp_comm_self
293 24 : CALL cp_blacs_env_create(context, para_env_self)
294 : CALL cp_fm_struct_create(fmstruct, para_env=para_env_self, context=context, &
295 : nrow_global=na_test, ncol_global=na_test, &
296 24 : nrow_block=nblk_test, ncol_block=nblk_test)
297 24 : CALL cp_cfm_create(matrix, fmstruct, name="elpa_check_mat")
298 24 : CALL cp_cfm_create(matrix_ref, fmstruct, name="elpa_check_ref")
299 24 : CALL cp_cfm_create(eigenvectors, fmstruct, name="elpa_check_vec")
300 :
301 : ! Fill with a dense random Hermitian matrix; the RNG stream starts from the
302 : ! library default seed, hence the matrix is identical on all processes.
303 24 : rng_stream = rng_stream_type(name="elpa_kernel_correctness", distribution_type=UNIFORM)
304 504 : DO j = 1, na_test
305 5040 : DO i = 1, j - 1
306 4560 : rnd_re = 2.0_dp*rng_stream%next() - 1.0_dp
307 4560 : rnd_im = 2.0_dp*rng_stream%next() - 1.0_dp
308 4560 : matrix%local_data(i, j) = CMPLX(rnd_re, rnd_im, KIND=dp)
309 5040 : matrix%local_data(j, i) = CONJG(matrix%local_data(i, j))
310 : END DO
311 480 : rnd_re = 2.0_dp*rng_stream%next() - 1.0_dp
312 504 : matrix%local_data(j, j) = CMPLX(rnd_re, 0.0_dp, KIND=dp)
313 : END DO
314 24 : CALL cp_cfm_to_cfm(matrix, matrix_ref) ! the solver destroys its input matrix
315 :
316 24 : ALLOCATE (eigenvalues(nev_test))
317 24 : CALL cp_cfm_diag_elpa_base(matrix, eigenvectors, eigenvalues, rdinfo)
318 :
319 : ! Check the orthogonality of the eigenvectors, |Z^H*Z - I|
320 : CALL cp_cfm_gemm("C", "N", nev_test, nev_test, na_test, (1.0_dp, 0.0_dp), &
321 24 : eigenvectors, eigenvectors, (0.0_dp, 0.0_dp), matrix)
322 24 : eps_ortho = 0.0_dp
323 456 : outer_ortho: DO i = 1, nev_test
324 8232 : DO j = 1, nev_test
325 7776 : test = ABS(matrix%local_data(i, j))
326 7776 : IF (i == j) test = ABS(matrix%local_data(i, j) - (1.0_dp, 0.0_dp))
327 7776 : IF (test > eps_ortho) eps_ortho = test
328 8208 : IF (eps_ortho > th) EXIT outer_ortho
329 : END DO
330 : END DO outer_ortho
331 :
332 : ! Check the eigenpair residuals, |A*Z - Z*diag(eigenvalues)|
333 : CALL cp_cfm_gemm("N", "N", na_test, nev_test, na_test, (1.0_dp, 0.0_dp), &
334 24 : matrix_ref, eigenvectors, (0.0_dp, 0.0_dp), matrix)
335 24 : eps_residual = 0.0_dp
336 504 : outer_res: DO i = 1, na_test
337 9144 : DO j = 1, nev_test
338 8640 : test = ABS(matrix%local_data(i, j) - eigenvectors%local_data(i, j)*eigenvalues(j))
339 8640 : IF (test > eps_residual) eps_residual = test
340 9120 : IF (eps_residual > th) EXIT outer_res
341 : END DO
342 : END DO outer_res
343 :
344 : ! Mis-compiled kernels produce grossly wrong results (deviations ~1e-1), not subtle rounding.
345 24 : elpa_c_broken = (eps_ortho > th) .OR. (eps_residual > th)
346 :
347 : ! Agree on the result to guarantee a consistent fallback on all processes.
348 24 : is_broken = 0
349 24 : IF (elpa_c_broken) is_broken = 1
350 24 : CALL para_env%max(is_broken)
351 24 : elpa_c_broken = is_broken == 1
352 24 : IF (elpa_c_broken) THEN
353 : message = "The ELPA complex kernel "//TRIM(get_elpa_c_kernel_name(elpa_c_kernel))// &
354 : " failed a runtime correctness check (orthogonality: "//TRIM(cp_to_string(eps_ortho))// &
355 : ", residual: "//TRIM(cp_to_string(eps_residual))//") and may return incorrect"// &
356 : " eigenvectors. Complex matrices fall back to ScaLAPACK. Consider the GENERIC"// &
357 0 : " kernel or rebuilding ELPA with -fno-tree-slp-vectorize."
358 0 : CALL cp_warn(__LOCATION__, TRIM(message))
359 : END IF
360 :
361 24 : DEALLOCATE (eigenvalues)
362 24 : CALL cp_cfm_release(eigenvectors)
363 24 : CALL cp_cfm_release(matrix_ref)
364 24 : CALL cp_cfm_release(matrix)
365 24 : CALL cp_fm_struct_release(fmstruct)
366 24 : CALL cp_blacs_env_release(context)
367 :
368 24 : CALL timestop(handle)
369 : #else
370 : MARK_USED(para_env)
371 : #endif
372 848 : END SUBROUTINE check_elpa_c_kernel_correctness
373 :
374 : ! **************************************************************************************************
375 : !> \brief Returns .TRUE. if the runtime check detected broken BLOCK2 kernels.
376 : !> \return ...
377 : ! **************************************************************************************************
378 72 : FUNCTION is_elpa_c_broken() RESULT(broken)
379 : LOGICAL :: broken
380 :
381 72 : broken = elpa_c_broken
382 72 : END FUNCTION is_elpa_c_broken
383 :
384 : ! **************************************************************************************************
385 : !> \brief Driver routine to diagonalize a CFM matrix with the ELPA library.
386 : !> \param matrix the matrix that is diagonalized
387 : !> \param eigenvectors eigenvectors of the input matrix
388 : !> \param eigenvalues eigenvalues of the input matrix
389 : ! **************************************************************************************************
390 72 : SUBROUTINE cp_cfm_diag_elpa(matrix, eigenvectors, eigenvalues)
391 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
392 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
393 :
394 : #if defined(__ELPA)
395 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_diag_elpa'
396 :
397 : INTEGER :: handle
398 : TYPE(cp_cfm_type) :: eigenvectors_new, matrix_new
399 : TYPE(cp_fm_redistribute_info) :: rdinfo
400 :
401 72 : CALL timeset(routineN, handle)
402 :
403 : ! Determine if the input matrix needs to be redistributed before diagonalization.
404 : ! Heuristics are used to determine the optimal number of CPUs for diagonalization.
405 : ! The redistributed matrix is stored in matrix_new, which is just a pointer
406 : ! to the original matrix if no redistribution is required.
407 : ! With ELPA, we have to make sure that all processor columns have nonzero width
408 : CALL cp_cfm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
409 72 : caller_is_elpa=.TRUE., redist_info=rdinfo)
410 :
411 : ! Call ELPA on CPUs that hold the new matrix
412 72 : IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
413 36 : CALL cp_cfm_diag_elpa_base(matrix_new, eigenvectors_new, eigenvalues, rdinfo)
414 : END IF
415 :
416 : ! Redistribute results and clean up
417 72 : CALL cp_cfm_redistribute_end(matrix, eigenvectors, eigenvalues, matrix_new, eigenvectors_new)
418 :
419 72 : CALL timestop(handle)
420 : #else
421 : eigenvalues = 0
422 : MARK_USED(matrix)
423 : MARK_USED(eigenvectors)
424 :
425 : CPABORT("CP2K compiled without the ELPA library.")
426 : #endif
427 72 : END SUBROUTINE cp_cfm_diag_elpa
428 :
429 : #if defined(__ELPA)
430 : ! **************************************************************************************************
431 : !> \brief Actual routine that calls ELPA to diagonalize a CFM matrix.
432 : !> \param matrix the matrix that is diagonalized
433 : !> \param eigenvectors eigenvectors of the input matrix
434 : !> \param eigenvalues eigenvalues of the input matrix
435 : !> \param rdinfo ...
436 : ! **************************************************************************************************
437 60 : SUBROUTINE cp_cfm_diag_elpa_base(matrix, eigenvectors, eigenvalues, rdinfo)
438 :
439 : TYPE(cp_cfm_type), INTENT(IN) :: matrix, eigenvectors
440 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
441 : TYPE(cp_fm_redistribute_info), INTENT(IN) :: rdinfo
442 :
443 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_diag_elpa_base'
444 :
445 : INTEGER :: handle
446 :
447 : CLASS(elpa_t), POINTER :: elpa_obj
448 : CHARACTER(len=default_string_length) :: kernel_name
449 : CHARACTER(len=2*default_string_length) :: message
450 : TYPE(mp_comm_type) :: group
451 : INTEGER :: fallback_kernel, &
452 : mypcol, myprow, n, &
453 : n_rows, n_cols, &
454 : nblk, neig, io_unit, &
455 : success
456 60 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eval
457 : TYPE(cp_blacs_env_type), POINTER :: context
458 : TYPE(cp_logger_type), POINTER :: logger
459 60 : INTEGER, DIMENSION(:), POINTER :: ncol_locals
460 : #if defined(__HAS_IEEE_EXCEPTIONS)
461 : LOGICAL, DIMENSION(5) :: halt
462 : #endif
463 :
464 60 : CALL timeset(routineN, handle)
465 60 : NULLIFY (logger)
466 60 : NULLIFY (ncol_locals)
467 :
468 60 : logger => cp_get_default_logger()
469 60 : io_unit = cp_logger_get_default_io_unit(logger)
470 :
471 60 : n = matrix%matrix_struct%nrow_global
472 60 : context => matrix%matrix_struct%context
473 60 : group = matrix%matrix_struct%para_env
474 :
475 60 : myprow = context%mepos(1)
476 60 : mypcol = context%mepos(2)
477 :
478 : ! elpa needs the full matrix
479 60 : CALL cp_cfm_uplo_to_full(matrix, eigenvectors)
480 :
481 : CALL cp_fm_struct_get(matrix%matrix_struct, &
482 : local_leading_dimension=n_rows, &
483 : ncol_local=n_cols, &
484 : nrow_block=nblk, &
485 60 : ncol_locals=ncol_locals)
486 :
487 : ! ELPA will fail in 'solve_tridi', with no useful error message, fail earlier
488 120 : IF (io_unit > 0 .AND. ANY(ncol_locals == 0)) THEN
489 0 : CALL rdinfo%write(io_unit)
490 0 : CPABORT("ELPA [pre-fail]: Problem contains processor column with zero width.")
491 : END IF
492 :
493 60 : neig = SIZE(eigenvalues, 1)
494 : ! ELPA's QR decomposition is only available for real matrices
495 :
496 60 : IF (io_unit > 0 .AND. elpa_print) THEN
497 : WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
498 29 : "ELPA| Matrix diagonalization information"
499 :
500 29 : kernel_name = get_elpa_c_kernel_name(elpa_c_kernel)
501 :
502 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
503 29 : "ELPA| Matrix order (NA) ", n, &
504 29 : "ELPA| Matrix block size (NBLK) ", nblk, &
505 29 : "ELPA| Number of eigenvectors (NEV) ", neig, &
506 29 : "ELPA| Local rows (LOCAL_NROWS) ", n_rows, &
507 58 : "ELPA| Local columns (LOCAL_NCOLS) ", n_cols
508 : WRITE (UNIT=io_unit, FMT="(T2,A,T61,A20)") &
509 29 : "ELPA| Kernel ", ADJUSTR(TRIM(kernel_name))
510 : END IF
511 :
512 : ! the full eigenvalues vector is needed
513 180 : ALLOCATE (eval(n))
514 :
515 60 : elpa_obj => elpa_allocate()
516 :
517 60 : CALL elpa_obj%set("na", n, success)
518 60 : CPASSERT(success == ELPA_OK)
519 :
520 60 : CALL elpa_obj%set("nev", neig, success)
521 60 : CPASSERT(success == ELPA_OK)
522 :
523 60 : CALL elpa_obj%set("local_nrows", n_rows, success)
524 60 : CPASSERT(success == ELPA_OK)
525 :
526 60 : CALL elpa_obj%set("local_ncols", n_cols, success)
527 60 : CPASSERT(success == ELPA_OK)
528 :
529 60 : CALL elpa_obj%set("nblk", nblk, success)
530 60 : CPASSERT(success == ELPA_OK)
531 :
532 60 : CALL elpa_obj%set("mpi_comm_parent", group%get_handle(), success)
533 60 : CPASSERT(success == ELPA_OK)
534 :
535 60 : CALL elpa_obj%set("process_row", myprow, success)
536 60 : CPASSERT(success == ELPA_OK)
537 :
538 60 : CALL elpa_obj%set("process_col", mypcol, success)
539 60 : CPASSERT(success == ELPA_OK)
540 :
541 60 : success = elpa_obj%setup()
542 60 : CPASSERT(success == ELPA_OK)
543 :
544 : CALL elpa_obj%set("solver", &
545 : MERGE(ELPA_SOLVER_1STAGE, ELPA_SOLVER_2STAGE, elpa_one_stage), &
546 108 : success)
547 60 : IF (success /= ELPA_OK) THEN
548 0 : CPABORT("Setting solver for ELPA failed")
549 : END IF
550 :
551 : ! enabling the GPU must happen before setting the kernel
552 0 : SELECT CASE (elpa_c_kernel)
553 : CASE (ELPA_2STAGE_COMPLEX_NVIDIA_GPU)
554 0 : CALL elpa_obj%set("nvidia-gpu", 1, success)
555 0 : CPASSERT(success == ELPA_OK)
556 : CASE (ELPA_2STAGE_COMPLEX_AMD_GPU)
557 0 : CALL elpa_obj%set("amd-gpu", 1, success)
558 0 : CPASSERT(success == ELPA_OK)
559 : CASE (ELPA_2STAGE_COMPLEX_INTEL_GPU_SYCL)
560 0 : CALL elpa_obj%set("intel-gpu", 1, success)
561 60 : CPASSERT(success == ELPA_OK)
562 : END SELECT
563 :
564 60 : IF (.NOT. elpa_one_stage) THEN
565 : ! Keep ELPA's configured default in case the requested kernel is unavailable.
566 48 : CALL elpa_obj%get("complex_kernel", fallback_kernel, success)
567 48 : CPASSERT(success == ELPA_OK)
568 :
569 48 : CALL elpa_obj%set("complex_kernel", elpa_c_kernel, success)
570 48 : IF (success /= ELPA_OK) THEN
571 : message = "Requested ELPA complex kernel "//TRIM(get_elpa_c_kernel_name(elpa_c_kernel))// &
572 0 : " is unavailable; falling back to "//TRIM(get_elpa_c_kernel_name(fallback_kernel))
573 0 : CALL cp_warn(__LOCATION__, TRIM(message))
574 0 : CALL elpa_obj%set("complex_kernel", fallback_kernel, success)
575 0 : CPASSERT(success == ELPA_OK)
576 : ! Avoid retrying the unavailable kernel for every subsequent diagonalization.
577 0 : elpa_c_kernel = fallback_kernel
578 : END IF
579 : END IF
580 :
581 : ! Set number of threads only when ELPA was built with OpenMP support.
582 60 : IF (elpa_obj%can_set("omp_threads", omp_get_max_threads()) == ELPA_OK) THEN
583 60 : CALL elpa_obj%set("omp_threads", omp_get_max_threads(), success)
584 60 : CPASSERT(success == ELPA_OK)
585 : END IF
586 :
587 : ! ELPA solver: calculate the Eigenvalues/vectors
588 : #if defined(__HAS_IEEE_EXCEPTIONS)
589 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
590 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
591 : #endif
592 60 : CALL elpa_obj%eigenvectors(matrix%local_data, eval, eigenvectors%local_data, success)
593 : #if defined(__HAS_IEEE_EXCEPTIONS)
594 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
595 : #endif
596 :
597 60 : IF (success /= ELPA_OK) THEN
598 0 : CPABORT("ELPA failed to diagonalize a matrix")
599 : END IF
600 :
601 60 : CALL elpa_deallocate(elpa_obj, success)
602 60 : CPASSERT(success == ELPA_OK)
603 :
604 1896 : eigenvalues(1:neig) = eval(1:neig)
605 60 : DEALLOCATE (eval)
606 :
607 60 : CALL timestop(handle)
608 :
609 180 : END SUBROUTINE cp_cfm_diag_elpa_base
610 : #endif
611 :
612 : END MODULE cp_cfm_elpa
|
