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
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE cp_fm_elpa
13 : USE cp_log_handling, ONLY: cp_to_string
14 : USE machine, ONLY: m_cpuid_static, &
15 : MACHINE_CPU_GENERIC, &
16 : MACHINE_X86_SSE4, &
17 : MACHINE_X86_AVX, &
18 : MACHINE_X86_AVX2, &
19 : MACHINE_X86_AVX512
20 : USE cp_blacs_env, ONLY: cp_blacs_env_type
21 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
22 : USE cp_fm_diag_utils, ONLY: cp_fm_redistribute_start, &
23 : cp_fm_redistribute_end, &
24 : cp_fm_redistribute_info
25 : USE cp_fm_struct, ONLY: cp_fm_struct_get
26 : USE cp_fm_types, ONLY: cp_fm_type, &
27 : cp_fm_to_fm, &
28 : cp_fm_release, &
29 : cp_fm_create, &
30 : cp_fm_write_info
31 : USE cp_log_handling, ONLY: cp_get_default_logger, &
32 : cp_logger_get_default_io_unit, &
33 : cp_logger_type
34 : USE kinds, ONLY: default_string_length, dp
35 : USE message_passing, ONLY: mp_comm_type
36 : USE OMP_LIB, ONLY: omp_get_max_threads
37 : #if defined(__HAS_IEEE_EXCEPTIONS)
38 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
39 : ieee_set_halting_mode, &
40 : IEEE_ALL
41 : #endif
42 : #include "../base/base_uses.f90"
43 :
44 : #if defined(__ELPA)
45 : USE elpa_constants, ONLY: ELPA_SOLVER_1STAGE, ELPA_SOLVER_2STAGE, ELPA_OK, &
46 : ELPA_2STAGE_REAL_INVALID, &
47 : ELPA_2STAGE_REAL_DEFAULT, &
48 : ELPA_2STAGE_REAL_GENERIC, &
49 : ELPA_2STAGE_REAL_GENERIC_SIMPLE, &
50 : ELPA_2STAGE_REAL_BGP, &
51 : ELPA_2STAGE_REAL_BGQ, &
52 : ELPA_2STAGE_REAL_SSE_ASSEMBLY, &
53 : ELPA_2STAGE_REAL_SSE_BLOCK2, &
54 : ELPA_2STAGE_REAL_SSE_BLOCK4, &
55 : ELPA_2STAGE_REAL_SSE_BLOCK6, &
56 : ELPA_2STAGE_REAL_AVX_BLOCK2, &
57 : ELPA_2STAGE_REAL_AVX_BLOCK4, &
58 : ELPA_2STAGE_REAL_AVX_BLOCK6, &
59 : ELPA_2STAGE_REAL_AVX2_BLOCK2, &
60 : ELPA_2STAGE_REAL_AVX2_BLOCK4, &
61 : ELPA_2STAGE_REAL_AVX2_BLOCK6, &
62 : ELPA_2STAGE_REAL_AVX512_BLOCK2, &
63 : ELPA_2STAGE_REAL_AVX512_BLOCK4, &
64 : ELPA_2STAGE_REAL_AVX512_BLOCK6, &
65 : ELPA_2STAGE_REAL_NVIDIA_GPU, &
66 : ELPA_2STAGE_REAL_AMD_GPU, &
67 : ELPA_2STAGE_REAL_INTEL_GPU_SYCL
68 :
69 : USE elpa, ONLY: elpa_t, elpa_init, elpa_uninit, &
70 : elpa_allocate, elpa_deallocate
71 : #endif
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_elpa'
78 :
79 : #if defined(__ELPA)
80 : INTEGER, DIMENSION(21), PARAMETER :: elpa_kernel_ids = [ &
81 : ELPA_2STAGE_REAL_INVALID, & ! auto
82 : ELPA_2STAGE_REAL_GENERIC, &
83 : ELPA_2STAGE_REAL_GENERIC_SIMPLE, &
84 : ELPA_2STAGE_REAL_BGP, &
85 : ELPA_2STAGE_REAL_BGQ, &
86 : ELPA_2STAGE_REAL_SSE_ASSEMBLY, &
87 : ELPA_2STAGE_REAL_SSE_BLOCK2, &
88 : ELPA_2STAGE_REAL_SSE_BLOCK4, &
89 : ELPA_2STAGE_REAL_SSE_BLOCK6, &
90 : ELPA_2STAGE_REAL_AVX_BLOCK2, &
91 : ELPA_2STAGE_REAL_AVX_BLOCK4, &
92 : ELPA_2STAGE_REAL_AVX_BLOCK6, &
93 : ELPA_2STAGE_REAL_AVX2_BLOCK2, &
94 : ELPA_2STAGE_REAL_AVX2_BLOCK4, &
95 : ELPA_2STAGE_REAL_AVX2_BLOCK6, &
96 : ELPA_2STAGE_REAL_AVX512_BLOCK2, &
97 : ELPA_2STAGE_REAL_AVX512_BLOCK4, &
98 : ELPA_2STAGE_REAL_AVX512_BLOCK6, &
99 : ELPA_2STAGE_REAL_NVIDIA_GPU, &
100 : ELPA_2STAGE_REAL_AMD_GPU, &
101 : ELPA_2STAGE_REAL_INTEL_GPU_SYCL]
102 :
103 : CHARACTER(len=14), DIMENSION(SIZE(elpa_kernel_ids)), PARAMETER :: &
104 : elpa_kernel_names = [CHARACTER(len=14) :: &
105 : "AUTO", &
106 : "GENERIC", &
107 : "GENERIC_SIMPLE", &
108 : "BGP", &
109 : "BGQ", &
110 : "SSE", &
111 : "SSE_BLOCK2", &
112 : "SSE_BLOCK4", &
113 : "SSE_BLOCK6", &
114 : "AVX_BLOCK2", &
115 : "AVX_BLOCK4", &
116 : "AVX_BLOCK6", &
117 : "AVX2_BLOCK2", &
118 : "AVX2_BLOCK4", &
119 : "AVX2_BLOCK6", &
120 : "AVX512_BLOCK2", &
121 : "AVX512_BLOCK4", &
122 : "AVX512_BLOCK6", &
123 : "NVIDIA_GPU", &
124 : "AMD_GPU", &
125 : "INTEL_GPU"]
126 :
127 : CHARACTER(len=44), DIMENSION(SIZE(elpa_kernel_ids)), PARAMETER :: &
128 : elpa_kernel_descriptions = [CHARACTER(len=44) :: &
129 : "Automatically selected kernel", &
130 : "Generic kernel", &
131 : "Simplified generic kernel", &
132 : "Kernel optimized for IBM BGP", &
133 : "Kernel optimized for IBM BGQ", &
134 : "Kernel optimized for x86_64/SSE", &
135 : "Kernel optimized for x86_64/SSE (block=2)", &
136 : "Kernel optimized for x86_64/SSE (block=4)", &
137 : "Kernel optimized for x86_64/SSE (block=6)", &
138 : "Kernel optimized for Intel AVX (block=2)", &
139 : "Kernel optimized for Intel AVX (block=4)", &
140 : "Kernel optimized for Intel AVX (block=6)", &
141 : "Kernel optimized for Intel AVX2 (block=2)", &
142 : "Kernel optimized for Intel AVX2 (block=4)", &
143 : "Kernel optimized for Intel AVX2 (block=6)", &
144 : "Kernel optimized for Intel AVX-512 (block=2)", &
145 : "Kernel optimized for Intel AVX-512 (block=4)", &
146 : "Kernel optimized for Intel AVX-512 (block=6)", &
147 : "Kernel targeting Nvidia GPUs", &
148 : "Kernel targeting AMD GPUs", &
149 : "Kernel targeting Intel GPUs"]
150 : #else
151 : INTEGER, DIMENSION(1), PARAMETER :: elpa_kernel_ids = [-1]
152 : CHARACTER(len=14), DIMENSION(1), PARAMETER :: elpa_kernel_names = ["AUTO"]
153 : CHARACTER(len=44), DIMENSION(1), PARAMETER :: elpa_kernel_descriptions = ["Automatically selected kernel"]
154 : #endif
155 :
156 : #if defined(__ELPA)
157 : INTEGER, SAVE :: elpa_kernel = elpa_kernel_ids(1) ! auto
158 : #endif
159 :
160 : ! elpa_qr_unsafe: disable block size limitations
161 : LOGICAL, SAVE :: elpa_qr_unsafe = .TRUE., &
162 : elpa_print = .FALSE., &
163 : elpa_qr = .FALSE.
164 :
165 : #if defined(__OFFLOAD_OPENCL)
166 : LOGICAL, SAVE :: elpa_one_stage = .TRUE.
167 : #else
168 : LOGICAL, SAVE :: elpa_one_stage = .FALSE.
169 : #endif
170 :
171 : PUBLIC :: cp_fm_diag_elpa, &
172 : set_elpa_kernel, &
173 : elpa_one_stage, &
174 : elpa_print, &
175 : elpa_qr, &
176 : elpa_kernel_ids, &
177 : elpa_kernel_names, &
178 : elpa_kernel_descriptions, &
179 : initialize_elpa_library, &
180 : finalize_elpa_library
181 :
182 : CONTAINS
183 :
184 : #if defined(__ELPA)
185 : ! **************************************************************************************************
186 : !> \brief Return a printable name for an ELPA real kernel.
187 : !> \param kernel ELPA real kernel id
188 : !> \return ...
189 : ! **************************************************************************************************
190 30 : FUNCTION get_elpa_kernel_name(kernel) RESULT(kernel_name)
191 : INTEGER, INTENT(IN) :: kernel
192 : CHARACTER(len=default_string_length) :: kernel_name
193 :
194 : INTEGER :: i
195 :
196 30 : kernel_name = "id: "//TRIM(ADJUSTL(cp_to_string(kernel)))
197 180 : DO i = 1, SIZE(elpa_kernel_ids)
198 180 : IF (elpa_kernel_ids(i) == kernel) THEN
199 30 : kernel_name = elpa_kernel_names(i)
200 30 : EXIT
201 : END IF
202 : END DO
203 30 : END FUNCTION get_elpa_kernel_name
204 : #endif
205 :
206 : ! **************************************************************************************************
207 : !> \brief Initialize the ELPA library
208 : !> \param one_stage ...
209 : !> \param qr ...
210 : !> \param should_print flag that determines if additional information
211 : !> is printed when the diagonalization routine is called.
212 : ! **************************************************************************************************
213 10426 : SUBROUTINE initialize_elpa_library(one_stage, qr, should_print)
214 : LOGICAL, INTENT(IN), OPTIONAL :: one_stage, qr, should_print
215 :
216 : #if defined(__ELPA)
217 10426 : IF (elpa_init(20180525) /= ELPA_OK) THEN
218 0 : CPABORT("The linked ELPA library does not support the required API version")
219 : END IF
220 10426 : IF (PRESENT(one_stage)) elpa_one_stage = one_stage
221 10426 : IF (PRESENT(should_print)) elpa_print = should_print
222 10426 : IF (PRESENT(qr)) elpa_qr = qr
223 : #else
224 : MARK_USED(one_stage)
225 : MARK_USED(qr)
226 : MARK_USED(should_print)
227 : CPABORT("Initialization of ELPA library requested but not enabled during build")
228 : #endif
229 10426 : END SUBROUTINE initialize_elpa_library
230 :
231 : ! **************************************************************************************************
232 : !> \brief Finalize the ELPA library
233 : ! **************************************************************************************************
234 10817 : SUBROUTINE finalize_elpa_library()
235 : #if defined(__ELPA)
236 10817 : CALL elpa_uninit()
237 : #else
238 : CPABORT("Finalization of ELPA library requested but not enabled during build")
239 : #endif
240 10817 : END SUBROUTINE finalize_elpa_library
241 :
242 : ! **************************************************************************************************
243 : !> \brief Sets the active ELPA kernel.
244 : !> \param requested_kernel one of the elpa_kernel_ids
245 : ! **************************************************************************************************
246 10426 : SUBROUTINE set_elpa_kernel(requested_kernel)
247 : INTEGER, INTENT(IN) :: requested_kernel
248 :
249 : #if defined(__ELPA)
250 : INTEGER :: cpuid
251 :
252 10426 : elpa_kernel = requested_kernel
253 :
254 : ! Resolve AUTO kernel.
255 10426 : IF (elpa_kernel == ELPA_2STAGE_REAL_INVALID) THEN
256 10422 : cpuid = m_cpuid_static()
257 0 : SELECT CASE (cpuid)
258 : CASE (MACHINE_CPU_GENERIC)
259 0 : elpa_kernel = ELPA_2STAGE_REAL_GENERIC
260 : CASE (MACHINE_X86_SSE4)
261 0 : elpa_kernel = ELPA_2STAGE_REAL_SSE_BLOCK4
262 : CASE (MACHINE_X86_AVX)
263 0 : elpa_kernel = ELPA_2STAGE_REAL_AVX_BLOCK4
264 : CASE (MACHINE_X86_AVX2)
265 10422 : elpa_kernel = ELPA_2STAGE_REAL_AVX2_BLOCK4
266 : CASE (MACHINE_X86_AVX512)
267 10422 : elpa_kernel = ELPA_2STAGE_REAL_AVX512_BLOCK4
268 : END SELECT
269 :
270 : ! Prefer GPU kernel if available.
271 : #if !defined(__NO_OFFLOAD_ELPA)
272 : #if defined(__OFFLOAD_CUDA)
273 : elpa_kernel = ELPA_2STAGE_REAL_NVIDIA_GPU
274 : #endif
275 : #if defined(__OFFLOAD_HIP)
276 : elpa_kernel = ELPA_2STAGE_REAL_AMD_GPU
277 : #endif
278 : #if defined(__OFFLOAD_OPENCL)
279 : elpa_kernel = ELPA_2STAGE_REAL_INTEL_GPU_SYCL
280 : #endif
281 : #endif
282 : ! If we could not find a suitable kernel then use ELPA_2STAGE_REAL_DEFAULT.
283 10422 : IF (elpa_kernel == ELPA_2STAGE_REAL_INVALID) THEN
284 0 : elpa_kernel = ELPA_2STAGE_REAL_DEFAULT
285 : END IF
286 : END IF
287 : #else
288 : MARK_USED(requested_kernel)
289 : #endif
290 10426 : END SUBROUTINE set_elpa_kernel
291 :
292 : ! **************************************************************************************************
293 : !> \brief Driver routine to diagonalize a FM matrix with the ELPA library.
294 : !> \param matrix the matrix that is diagonalized
295 : !> \param eigenvectors eigenvectors of the input matrix
296 : !> \param eigenvalues eigenvalues of the input matrix
297 : ! **************************************************************************************************
298 14420 : SUBROUTINE cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
299 : TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
300 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
301 :
302 : #if defined(__ELPA)
303 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_diag_elpa'
304 :
305 : INTEGER :: handle
306 : TYPE(cp_fm_type) :: eigenvectors_new, matrix_new
307 : TYPE(cp_fm_redistribute_info) :: rdinfo
308 :
309 14420 : CALL timeset(routineN, handle)
310 :
311 : ! Determine if the input matrix needs to be redistributed before diagonalization.
312 : ! Heuristics are used to determine the optimal number of CPUs for diagonalization.
313 : ! The redistributed matrix is stored in matrix_new, which is just a pointer
314 : ! to the original matrix if no redistribution is required.
315 : ! With ELPA, we have to make sure that all processor columns have nonzero width
316 : CALL cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, &
317 14420 : caller_is_elpa=.TRUE., redist_info=rdinfo)
318 :
319 : ! Call ELPA on CPUs that hold the new matrix
320 14420 : IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
321 9285 : CALL cp_fm_diag_elpa_base(matrix_new, eigenvectors_new, eigenvalues, rdinfo)
322 : END IF
323 :
324 : ! Redistribute results and clean up
325 14420 : CALL cp_fm_redistribute_end(matrix, eigenvectors, eigenvalues, matrix_new, eigenvectors_new)
326 :
327 14420 : CALL timestop(handle)
328 : #else
329 : eigenvalues = 0
330 : MARK_USED(matrix)
331 : MARK_USED(eigenvectors)
332 :
333 : CPABORT("CP2K compiled without the ELPA library.")
334 : #endif
335 14420 : END SUBROUTINE cp_fm_diag_elpa
336 :
337 : #if defined(__ELPA)
338 : ! **************************************************************************************************
339 : !> \brief Actual routine that calls ELPA to diagonalize a FM matrix.
340 : !> \param matrix the matrix that is diagonalized
341 : !> \param eigenvectors eigenvectors of the input matrix
342 : !> \param eigenvalues eigenvalues of the input matrix
343 : !> \param rdinfo ...
344 : ! **************************************************************************************************
345 9285 : SUBROUTINE cp_fm_diag_elpa_base(matrix, eigenvectors, eigenvalues, rdinfo)
346 :
347 : TYPE(cp_fm_type), INTENT(IN) :: matrix, eigenvectors
348 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
349 : TYPE(cp_fm_redistribute_info), INTENT(IN) :: rdinfo
350 :
351 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_diag_elpa_base'
352 :
353 : INTEGER :: handle
354 :
355 : CLASS(elpa_t), POINTER :: elpa_obj
356 : CHARACTER(len=default_string_length) :: kernel_name, message
357 : TYPE(mp_comm_type) :: group
358 : INTEGER :: fallback_kernel, &
359 : mypcol, myprow, n, &
360 : n_rows, n_cols, &
361 : nblk, neig, io_unit, &
362 : success
363 : LOGICAL :: use_qr, check_eigenvalues
364 9285 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: eval, eval_noqr
365 : TYPE(cp_blacs_env_type), POINTER :: context
366 : TYPE(cp_fm_type) :: matrix_noqr, eigenvectors_noqr
367 : TYPE(cp_logger_type), POINTER :: logger
368 : REAL(KIND=dp), PARAMETER :: th = 1.0E-14_dp
369 9285 : INTEGER, DIMENSION(:), POINTER :: ncol_locals
370 : #if defined(__HAS_IEEE_EXCEPTIONS)
371 : LOGICAL, DIMENSION(5) :: halt
372 : #endif
373 :
374 9285 : CALL timeset(routineN, handle)
375 9285 : NULLIFY (logger)
376 9285 : NULLIFY (ncol_locals)
377 :
378 9285 : check_eigenvalues = .FALSE.
379 :
380 9285 : logger => cp_get_default_logger()
381 9285 : io_unit = cp_logger_get_default_io_unit(logger)
382 :
383 9285 : n = matrix%matrix_struct%nrow_global
384 9285 : context => matrix%matrix_struct%context
385 9285 : group = matrix%matrix_struct%para_env
386 :
387 9285 : myprow = context%mepos(1)
388 9285 : mypcol = context%mepos(2)
389 :
390 : ! elpa needs the full matrix
391 9285 : CALL cp_fm_uplo_to_full(matrix, eigenvectors)
392 :
393 : CALL cp_fm_struct_get(matrix%matrix_struct, &
394 : local_leading_dimension=n_rows, &
395 : ncol_local=n_cols, &
396 : nrow_block=nblk, &
397 9285 : ncol_locals=ncol_locals)
398 :
399 : ! ELPA will fail in 'solve_tridi', with no useful error message, fail earlier
400 18576 : IF (io_unit > 0 .AND. ANY(ncol_locals == 0)) THEN
401 0 : CALL rdinfo%write(io_unit)
402 0 : CALL cp_fm_write_info(matrix, io_unit)
403 0 : CPABORT("ELPA [pre-fail]: Problem contains processor column with zero width.")
404 : END IF
405 :
406 9285 : neig = SIZE(eigenvalues, 1)
407 : ! Decide if matrix is suitable for ELPA to use QR
408 : ! The definition of what is considered a suitable matrix depends on the ELPA version
409 : ! The relevant ELPA files to check are
410 : ! - Proper matrix order: src/elpa2/elpa2_template.F90
411 : ! - Proper block size: test/Fortran/test.F90
412 : ! Note that the names of these files might change in different ELPA versions
413 : ! Matrix order must be even
414 9285 : use_qr = elpa_qr .AND. (MODULO(n, 2) == 0)
415 : ! Matrix order and block size must be greater than or equal to 64
416 9285 : IF (.NOT. elpa_qr_unsafe) THEN
417 0 : use_qr = use_qr .AND. (n >= 64) .AND. (nblk >= 64)
418 : END IF
419 :
420 : ! Check if eigenvalues computed with elpa_qr_unsafe should be verified
421 9285 : IF (use_qr .AND. elpa_qr_unsafe .AND. elpa_print) THEN
422 40 : check_eigenvalues = .TRUE.
423 : END IF
424 :
425 9285 : CALL matrix%matrix_struct%para_env%bcast(check_eigenvalues)
426 :
427 9285 : IF (check_eigenvalues) THEN
428 : ! Allocate and initialize needed temporaries to compute eigenvalues without ELPA QR
429 120 : ALLOCATE (eval_noqr(n))
430 40 : CALL cp_fm_create(matrix=matrix_noqr, matrix_struct=matrix%matrix_struct)
431 40 : CALL cp_fm_to_fm(matrix, matrix_noqr)
432 40 : CALL cp_fm_create(matrix=eigenvectors_noqr, matrix_struct=eigenvectors%matrix_struct)
433 80 : CALL cp_fm_uplo_to_full(matrix_noqr, eigenvectors_noqr)
434 : END IF
435 :
436 9285 : IF (io_unit > 0 .AND. elpa_print) THEN
437 : WRITE (UNIT=io_unit, FMT="(/,T2,A)") &
438 30 : "ELPA| Matrix diagonalization information"
439 :
440 30 : kernel_name = get_elpa_kernel_name(elpa_kernel)
441 :
442 : WRITE (UNIT=io_unit, FMT="(T2,A,T71,I10)") &
443 30 : "ELPA| Matrix order (NA) ", n, &
444 30 : "ELPA| Matrix block size (NBLK) ", nblk, &
445 30 : "ELPA| Number of eigenvectors (NEV) ", neig, &
446 30 : "ELPA| Local rows (LOCAL_NROWS) ", n_rows, &
447 60 : "ELPA| Local columns (LOCAL_NCOLS) ", n_cols
448 : WRITE (UNIT=io_unit, FMT="(T2,A,T61,A20)") &
449 30 : "ELPA| Kernel ", ADJUSTR(TRIM(kernel_name))
450 30 : IF (elpa_qr) THEN
451 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
452 20 : "ELPA| QR step requested ", "YES"
453 : ELSE
454 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
455 10 : "ELPA| QR step requested ", "NO"
456 : END IF
457 :
458 30 : IF (elpa_qr) THEN
459 20 : IF (use_qr) THEN
460 : WRITE (UNIT=io_unit, FMT="(T2,A,T78,A3)") &
461 20 : "ELPA| Matrix is suitable for QR ", "YES"
462 : ELSE
463 : WRITE (UNIT=io_unit, FMT="(T2,A,T79,A2)") &
464 0 : "ELPA| Matrix is suitable for QR ", "NO"
465 : END IF
466 20 : IF (.NOT. use_qr) THEN
467 0 : IF (MODULO(n, 2) /= 0) THEN
468 : WRITE (UNIT=io_unit, FMT="(T2,A)") &
469 0 : "ELPA| Matrix order is NOT even"
470 : END IF
471 0 : IF ((nblk < 64) .AND. (.NOT. elpa_qr_unsafe)) THEN
472 : WRITE (UNIT=io_unit, FMT="(T2,A)") &
473 0 : "ELPA| Matrix block size is NOT 64 or greater"
474 : END IF
475 : ELSE
476 20 : IF ((nblk < 64) .AND. elpa_qr_unsafe) THEN
477 : WRITE (UNIT=io_unit, FMT="(T2,A)") &
478 10 : "ELPA| Matrix block size check was bypassed"
479 : END IF
480 : END IF
481 : END IF
482 : END IF
483 :
484 : ! the full eigenvalues vector is needed
485 27855 : ALLOCATE (eval(n))
486 :
487 9285 : elpa_obj => elpa_allocate()
488 :
489 9285 : CALL elpa_obj%set("na", n, success)
490 9285 : CPASSERT(success == ELPA_OK)
491 :
492 9285 : CALL elpa_obj%set("nev", neig, success)
493 9285 : CPASSERT(success == ELPA_OK)
494 :
495 9285 : CALL elpa_obj%set("local_nrows", n_rows, success)
496 9285 : CPASSERT(success == ELPA_OK)
497 :
498 9285 : CALL elpa_obj%set("local_ncols", n_cols, success)
499 9285 : CPASSERT(success == ELPA_OK)
500 :
501 9285 : CALL elpa_obj%set("nblk", nblk, success)
502 9285 : CPASSERT(success == ELPA_OK)
503 :
504 9285 : CALL elpa_obj%set("mpi_comm_parent", group%get_handle(), success)
505 9285 : CPASSERT(success == ELPA_OK)
506 :
507 9285 : CALL elpa_obj%set("process_row", myprow, success)
508 9285 : CPASSERT(success == ELPA_OK)
509 :
510 9285 : CALL elpa_obj%set("process_col", mypcol, success)
511 9285 : CPASSERT(success == ELPA_OK)
512 :
513 9285 : success = elpa_obj%setup()
514 9285 : CPASSERT(success == ELPA_OK)
515 :
516 : CALL elpa_obj%set("solver", &
517 : MERGE(ELPA_SOLVER_1STAGE, ELPA_SOLVER_2STAGE, elpa_one_stage), &
518 18550 : success)
519 9285 : IF (success /= ELPA_OK) THEN
520 0 : CPABORT("Setting solver for ELPA failed")
521 : END IF
522 :
523 : ! enabling the GPU must happen before setting the kernel
524 0 : SELECT CASE (elpa_kernel)
525 : CASE (ELPA_2STAGE_REAL_NVIDIA_GPU)
526 0 : CALL elpa_obj%set("nvidia-gpu", 1, success)
527 0 : CPASSERT(success == ELPA_OK)
528 : CASE (ELPA_2STAGE_REAL_AMD_GPU)
529 0 : CALL elpa_obj%set("amd-gpu", 1, success)
530 0 : CPASSERT(success == ELPA_OK)
531 : CASE (ELPA_2STAGE_REAL_INTEL_GPU_SYCL)
532 0 : CALL elpa_obj%set("intel-gpu", 1, success)
533 9285 : CPASSERT(success == ELPA_OK)
534 : END SELECT
535 :
536 9285 : IF (.NOT. elpa_one_stage) THEN
537 : ! Keep ELPA's configured default in case the requested kernel is unavailable.
538 9265 : CALL elpa_obj%get("real_kernel", fallback_kernel, success)
539 9265 : CPASSERT(success == ELPA_OK)
540 :
541 9265 : CALL elpa_obj%set("real_kernel", elpa_kernel, success)
542 9265 : IF (success /= ELPA_OK) THEN
543 : message = "Requested ELPA real kernel "//TRIM(get_elpa_kernel_name(elpa_kernel))// &
544 0 : " is unavailable; falling back to "//TRIM(get_elpa_kernel_name(fallback_kernel))
545 0 : CPWARN(TRIM(message))
546 0 : CALL elpa_obj%set("real_kernel", fallback_kernel, success)
547 0 : CPASSERT(success == ELPA_OK)
548 : ! Avoid retrying the unavailable kernel for every subsequent diagonalization.
549 0 : elpa_kernel = fallback_kernel
550 : END IF
551 :
552 9265 : IF (use_qr) THEN
553 40 : CALL elpa_obj%set("qr", 1, success)
554 40 : CPASSERT(success == ELPA_OK)
555 : END IF
556 : END IF
557 :
558 : ! Set number of threads only when ELPA was built with OpenMP support.
559 9285 : IF (elpa_obj%can_set("omp_threads", omp_get_max_threads()) == ELPA_OK) THEN
560 9285 : CALL elpa_obj%set("omp_threads", omp_get_max_threads(), success)
561 9285 : CPASSERT(success == ELPA_OK)
562 : END IF
563 :
564 : ! ELPA solver: calculate the Eigenvalues/vectors
565 : #if defined(__HAS_IEEE_EXCEPTIONS)
566 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
567 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
568 : #endif
569 9285 : CALL elpa_obj%eigenvectors(matrix%local_data, eval, eigenvectors%local_data, success)
570 : #if defined(__HAS_IEEE_EXCEPTIONS)
571 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
572 : #endif
573 :
574 9285 : IF (success /= ELPA_OK) THEN
575 0 : CPABORT("ELPA failed to diagonalize a matrix")
576 : END IF
577 :
578 9285 : IF (check_eigenvalues) THEN
579 : ! run again without QR
580 40 : CALL elpa_obj%set("qr", 0, success)
581 40 : CPASSERT(success == ELPA_OK)
582 :
583 40 : CALL elpa_obj%eigenvectors(matrix_noqr%local_data, eval_noqr, eigenvectors_noqr%local_data, success)
584 40 : IF (success /= ELPA_OK) THEN
585 0 : CPABORT("ELPA failed to diagonalize a matrix even without QR decomposition")
586 : END IF
587 :
588 3400 : IF (ANY(ABS(eval(1:neig) - eval_noqr(1:neig)) > th)) THEN
589 0 : CPABORT("ELPA failed to calculate Eigenvalues with ELPA's QR decomposition")
590 : END IF
591 :
592 40 : DEALLOCATE (eval_noqr)
593 40 : CALL cp_fm_release(matrix_noqr)
594 40 : CALL cp_fm_release(eigenvectors_noqr)
595 : END IF
596 :
597 9285 : CALL elpa_deallocate(elpa_obj, success)
598 9285 : CPASSERT(success == ELPA_OK)
599 :
600 770780 : eigenvalues(1:neig) = eval(1:neig)
601 9285 : DEALLOCATE (eval)
602 :
603 9285 : CALL timestop(handle)
604 :
605 37140 : END SUBROUTINE cp_fm_diag_elpa_base
606 : #endif
607 :
608 : END MODULE cp_fm_elpa
|
