(git:21ef868)
Loading...
Searching...
No Matches
cp_fm_elpa.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
182CONTAINS
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 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 kernel_name = "id: "//trim(adjustl(cp_to_string(kernel)))
197 DO i = 1, SIZE(elpa_kernel_ids)
198 IF (elpa_kernel_ids(i) == kernel) THEN
199 kernel_name = elpa_kernel_names(i)
200 EXIT
201 END IF
202 END DO
203 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 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 IF (elpa_init(20180525) /= elpa_ok) THEN
218 cpabort("The linked ELPA library does not support the required API version")
219 END IF
220 IF (PRESENT(one_stage)) elpa_one_stage = one_stage
221 IF (PRESENT(should_print)) elpa_print = should_print
222 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 END SUBROUTINE initialize_elpa_library
230
231! **************************************************************************************************
232!> \brief Finalize the ELPA library
233! **************************************************************************************************
234 SUBROUTINE finalize_elpa_library()
235#if defined(__ELPA)
236 CALL elpa_uninit()
237#else
238 cpabort("Finalization of ELPA library requested but not enabled during build")
239#endif
240 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 SUBROUTINE set_elpa_kernel(requested_kernel)
247 INTEGER, INTENT(IN) :: requested_kernel
248
249#if defined(__ELPA)
250 INTEGER :: cpuid
251
252 elpa_kernel = requested_kernel
253
254 ! Resolve AUTO kernel.
255 IF (elpa_kernel == elpa_2stage_real_invalid) THEN
256 cpuid = m_cpuid_static()
257 SELECT CASE (cpuid)
258 CASE (machine_cpu_generic)
259 elpa_kernel = elpa_2stage_real_generic
260 CASE (machine_x86_sse4)
261 elpa_kernel = elpa_2stage_real_sse_block4
262 CASE (machine_x86_avx)
263 elpa_kernel = elpa_2stage_real_avx_block4
264 CASE (machine_x86_avx2)
265 elpa_kernel = elpa_2stage_real_avx2_block4
266 CASE (machine_x86_avx512)
267 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 IF (elpa_kernel == elpa_2stage_real_invalid) THEN
284 elpa_kernel = elpa_2stage_real_default
285 END IF
286 END IF
287#else
288 mark_used(requested_kernel)
289#endif
290 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 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 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 caller_is_elpa=.true., redist_info=rdinfo)
318
319 ! Call ELPA on CPUs that hold the new matrix
320 IF (ASSOCIATED(matrix_new%matrix_struct)) THEN
321 CALL cp_fm_diag_elpa_base(matrix_new, eigenvectors_new, eigenvalues, rdinfo)
322 END IF
323
324 ! Redistribute results and clean up
325 CALL cp_fm_redistribute_end(matrix, eigenvectors, eigenvalues, matrix_new, eigenvectors_new)
326
327 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: ncol_locals
370#if defined(__HAS_IEEE_EXCEPTIONS)
371 LOGICAL, DIMENSION(5) :: halt
372#endif
373
374 CALL timeset(routinen, handle)
375 NULLIFY (logger)
376 NULLIFY (ncol_locals)
377
378 check_eigenvalues = .false.
379
380 logger => cp_get_default_logger()
381 io_unit = cp_logger_get_default_io_unit(logger)
382
383 n = matrix%matrix_struct%nrow_global
384 context => matrix%matrix_struct%context
385 group = matrix%matrix_struct%para_env
386
387 myprow = context%mepos(1)
388 mypcol = context%mepos(2)
389
390 ! elpa needs the full matrix
391 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 ncol_locals=ncol_locals)
398
399 ! ELPA will fail in 'solve_tridi', with no useful error message, fail earlier
400 IF (io_unit > 0 .AND. any(ncol_locals == 0)) THEN
401 CALL rdinfo%write(io_unit)
402 CALL cp_fm_write_info(matrix, io_unit)
403 cpabort("ELPA [pre-fail]: Problem contains processor column with zero width.")
404 END IF
405
406 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 use_qr = elpa_qr .AND. (modulo(n, 2) == 0)
415 ! Matrix order and block size must be greater than or equal to 64
416 IF (.NOT. elpa_qr_unsafe) THEN
417 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 IF (use_qr .AND. elpa_qr_unsafe .AND. elpa_print) THEN
422 check_eigenvalues = .true.
423 END IF
424
425 CALL matrix%matrix_struct%para_env%bcast(check_eigenvalues)
426
427 IF (check_eigenvalues) THEN
428 ! Allocate and initialize needed temporaries to compute eigenvalues without ELPA QR
429 ALLOCATE (eval_noqr(n))
430 CALL cp_fm_create(matrix=matrix_noqr, matrix_struct=matrix%matrix_struct)
431 CALL cp_fm_to_fm(matrix, matrix_noqr)
432 CALL cp_fm_create(matrix=eigenvectors_noqr, matrix_struct=eigenvectors%matrix_struct)
433 CALL cp_fm_uplo_to_full(matrix_noqr, eigenvectors_noqr)
434 END IF
435
436 IF (io_unit > 0 .AND. elpa_print) THEN
437 WRITE (unit=io_unit, fmt="(/,T2,A)") &
438 "ELPA| Matrix diagonalization information"
439
440 kernel_name = get_elpa_kernel_name(elpa_kernel)
441
442 WRITE (unit=io_unit, fmt="(T2,A,T71,I10)") &
443 "ELPA| Matrix order (NA) ", n, &
444 "ELPA| Matrix block size (NBLK) ", nblk, &
445 "ELPA| Number of eigenvectors (NEV) ", neig, &
446 "ELPA| Local rows (LOCAL_NROWS) ", n_rows, &
447 "ELPA| Local columns (LOCAL_NCOLS) ", n_cols
448 WRITE (unit=io_unit, fmt="(T2,A,T61,A20)") &
449 "ELPA| Kernel ", adjustr(trim(kernel_name))
450 IF (elpa_qr) THEN
451 WRITE (unit=io_unit, fmt="(T2,A,T78,A3)") &
452 "ELPA| QR step requested ", "YES"
453 ELSE
454 WRITE (unit=io_unit, fmt="(T2,A,T79,A2)") &
455 "ELPA| QR step requested ", "NO"
456 END IF
457
458 IF (elpa_qr) THEN
459 IF (use_qr) THEN
460 WRITE (unit=io_unit, fmt="(T2,A,T78,A3)") &
461 "ELPA| Matrix is suitable for QR ", "YES"
462 ELSE
463 WRITE (unit=io_unit, fmt="(T2,A,T79,A2)") &
464 "ELPA| Matrix is suitable for QR ", "NO"
465 END IF
466 IF (.NOT. use_qr) THEN
467 IF (modulo(n, 2) /= 0) THEN
468 WRITE (unit=io_unit, fmt="(T2,A)") &
469 "ELPA| Matrix order is NOT even"
470 END IF
471 IF ((nblk < 64) .AND. (.NOT. elpa_qr_unsafe)) THEN
472 WRITE (unit=io_unit, fmt="(T2,A)") &
473 "ELPA| Matrix block size is NOT 64 or greater"
474 END IF
475 ELSE
476 IF ((nblk < 64) .AND. elpa_qr_unsafe) THEN
477 WRITE (unit=io_unit, fmt="(T2,A)") &
478 "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 ALLOCATE (eval(n))
486
487 elpa_obj => elpa_allocate()
488
489 CALL elpa_obj%set("na", n, success)
490 cpassert(success == elpa_ok)
491
492 CALL elpa_obj%set("nev", neig, success)
493 cpassert(success == elpa_ok)
494
495 CALL elpa_obj%set("local_nrows", n_rows, success)
496 cpassert(success == elpa_ok)
497
498 CALL elpa_obj%set("local_ncols", n_cols, success)
499 cpassert(success == elpa_ok)
500
501 CALL elpa_obj%set("nblk", nblk, success)
502 cpassert(success == elpa_ok)
503
504 CALL elpa_obj%set("mpi_comm_parent", group%get_handle(), success)
505 cpassert(success == elpa_ok)
506
507 CALL elpa_obj%set("process_row", myprow, success)
508 cpassert(success == elpa_ok)
509
510 CALL elpa_obj%set("process_col", mypcol, success)
511 cpassert(success == elpa_ok)
512
513 success = elpa_obj%setup()
514 cpassert(success == elpa_ok)
515
516 CALL elpa_obj%set("solver", &
517 merge(elpa_solver_1stage, elpa_solver_2stage, elpa_one_stage), &
518 success)
519 IF (success /= elpa_ok) THEN
520 cpabort("Setting solver for ELPA failed")
521 END IF
522
523 ! enabling the GPU must happen before setting the kernel
524 SELECT CASE (elpa_kernel)
525 CASE (elpa_2stage_real_nvidia_gpu)
526 CALL elpa_obj%set("nvidia-gpu", 1, success)
527 cpassert(success == elpa_ok)
528 CASE (elpa_2stage_real_amd_gpu)
529 CALL elpa_obj%set("amd-gpu", 1, success)
530 cpassert(success == elpa_ok)
531 CASE (elpa_2stage_real_intel_gpu_sycl)
532 CALL elpa_obj%set("intel-gpu", 1, success)
533 cpassert(success == elpa_ok)
534 END SELECT
535
536 IF (.NOT. elpa_one_stage) THEN
537 ! Keep ELPA's configured default in case the requested kernel is unavailable.
538 CALL elpa_obj%get("real_kernel", fallback_kernel, success)
539 cpassert(success == elpa_ok)
540
541 CALL elpa_obj%set("real_kernel", elpa_kernel, success)
542 IF (success /= elpa_ok) THEN
543 message = "Requested ELPA real kernel "//trim(get_elpa_kernel_name(elpa_kernel))// &
544 " is unavailable; falling back to "//trim(get_elpa_kernel_name(fallback_kernel))
545 cpwarn(trim(message))
546 CALL elpa_obj%set("real_kernel", fallback_kernel, success)
547 cpassert(success == elpa_ok)
548 ! Avoid retrying the unavailable kernel for every subsequent diagonalization.
549 elpa_kernel = fallback_kernel
550 END IF
551
552 IF (use_qr) THEN
553 CALL elpa_obj%set("qr", 1, success)
554 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 IF (elpa_obj%can_set("omp_threads", omp_get_max_threads()) == elpa_ok) THEN
560 CALL elpa_obj%set("omp_threads", omp_get_max_threads(), success)
561 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 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 IF (success /= elpa_ok) THEN
575 cpabort("ELPA failed to diagonalize a matrix")
576 END IF
577
578 IF (check_eigenvalues) THEN
579 ! run again without QR
580 CALL elpa_obj%set("qr", 0, success)
581 cpassert(success == elpa_ok)
582
583 CALL elpa_obj%eigenvectors(matrix_noqr%local_data, eval_noqr, eigenvectors_noqr%local_data, success)
584 IF (success /= elpa_ok) THEN
585 cpabort("ELPA failed to diagonalize a matrix even without QR decomposition")
586 END IF
587
588 IF (any(abs(eval(1:neig) - eval_noqr(1:neig)) > th)) THEN
589 cpabort("ELPA failed to calculate Eigenvalues with ELPA's QR decomposition")
590 END IF
591
592 DEALLOCATE (eval_noqr)
593 CALL cp_fm_release(matrix_noqr)
594 CALL cp_fm_release(eigenvectors_noqr)
595 END IF
596
597 CALL elpa_deallocate(elpa_obj, success)
598 cpassert(success == elpa_ok)
599
600 eigenvalues(1:neig) = eval(1:neig)
601 DEALLOCATE (eval)
602
603 CALL timestop(handle)
604
605 END SUBROUTINE cp_fm_diag_elpa_base
606#endif
607
608END MODULE cp_fm_elpa
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
machine::m_cpuid_static
Target architecture or instruction set extension according to compiler target flags.
Definition machine.F:90
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_diag_utils
Auxiliary tools to redistribute cp_fm_type matrices before and after diagonalization....
Definition cp_fm_diag_utils.F:16
cp_fm_diag_utils::cp_fm_redistribute_end
subroutine, public cp_fm_redistribute_end(matrix, eigenvectors, eig, matrix_new, eigenvectors_new)
Redistributes eigenvectors and eigenvalues back to the original communicator group.
Definition cp_fm_diag_utils.F:463
cp_fm_diag_utils::cp_fm_redistribute_start
subroutine, public cp_fm_redistribute_start(matrix, eigenvectors, matrix_new, eigenvectors_new, caller_is_elpa, redist_info)
Determines the optimal number of CPUs for matrix diagonalization and redistributes the input matrices...
Definition cp_fm_diag_utils.F:270
cp_fm_elpa
Wrapper for ELPA.
Definition cp_fm_elpa.F:12
cp_fm_elpa::set_elpa_kernel
subroutine, public set_elpa_kernel(requested_kernel)
Sets the active ELPA kernel.
Definition cp_fm_elpa.F:247
cp_fm_elpa::elpa_qr
logical, save, public elpa_qr
Definition cp_fm_elpa.F:161
cp_fm_elpa::cp_fm_diag_elpa
subroutine, public cp_fm_diag_elpa(matrix, eigenvectors, eigenvalues)
Driver routine to diagonalize a FM matrix with the ELPA library.
Definition cp_fm_elpa.F:299
cp_fm_elpa::elpa_print
logical, save, public elpa_print
Definition cp_fm_elpa.F:161
cp_fm_elpa::finalize_elpa_library
subroutine, public finalize_elpa_library()
Finalize the ELPA library.
Definition cp_fm_elpa.F:235
cp_fm_elpa::elpa_one_stage
logical, save, public elpa_one_stage
Definition cp_fm_elpa.F:168
cp_fm_elpa::initialize_elpa_library
subroutine, public initialize_elpa_library(one_stage, qr, should_print)
Initialize the ELPA library.
Definition cp_fm_elpa.F:214
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_get
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
Definition cp_fm_struct.F:444
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_write_info
subroutine, public cp_fm_write_info(matrix, io_unit)
Write nicely formatted info about the FM to the given I/O unit (including the underlying FM struct)
Definition cp_fm_types.F:1070
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::machine_x86_avx
integer, parameter, public machine_x86_avx
Definition machine.F:57
machine::machine_x86_sse4
integer, parameter, public machine_x86_sse4
Definition machine.F:57
machine::machine_cpu_generic
integer, parameter, public machine_cpu_generic
Definition machine.F:57
machine::machine_x86_avx2
integer, parameter, public machine_x86_avx2
Definition machine.F:57
machine::machine_x86_avx512
integer, parameter, public machine_x86_avx512
Definition machine.F:57
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_fm_diag_utils::cp_fm_redistribute_info
Definition cp_fm_diag_utils.F:42
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:161