Disable IEEE exceptions for LAPACK

- FP-exceptions for ScaLAPACK are disabled since a long time.
- __HAS_IEEE_EXCEPTIONS denotes intrinsic module available.
- Harmonized some CFM/FM code format/comments.
This commit is contained in:
Hans Pabst 2025-10-24 15:12:29 +02:00
parent a6dc286aa2
commit 9c0f83189b
4 changed files with 35 additions and 36 deletions

View file

@ -114,9 +114,9 @@ CONTAINS
INTEGER, DIMENSION(9) :: descm, descv
COMPLEX(KIND=dp), DIMENSION(:, :), &
POINTER :: v
#endif
#if defined (__HAS_IEEE_EXCEPTIONS)
LOGICAL, DIMENSION(5) :: halt
#endif
#endif
CALL timeset(routineN, handle)
@ -152,29 +152,26 @@ CONTAINS
DEALLOCATE (iwork, rwork, work)
ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
#if defined(__parallel)
! Scalapack takes advantage of IEEE754 exceptions for speedup.
! (Sca-)LAPACK takes advantage of IEEE754 exceptions for speedup.
! Therefore, we disable floating point traps temporarily.
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_get_halting_mode(IEEE_ALL, halt)
CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
#endif
#if defined(__parallel)
CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_set_halting_mode(IEEE_ALL, halt)
#endif
#else
CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
eigenvectors%local_data = matrix%local_data
#endif
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_set_halting_mode(IEEE_ALL, halt)
#endif
DEALLOCATE (iwork, rwork, work)
IF (info /= 0) &
CPABORT("Diagonalisation of a complex matrix failed")
IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
CALL timestop(handle)

View file

@ -108,7 +108,7 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE cp_cfm_pzpotri_dlaf
! **************************************************************************************************
! **************************************************************************************************
!> \brief DLA-Future eigensolver for complex Hermitian matrices
!> \param matrix ...
!> \param eigenvectors ...
@ -210,6 +210,7 @@ CONTAINS
! Use eigenvectors matrix as workspace
CALL cp_cfm_uplo_to_full(matrix, eigenvectors)
! Create DLAF grid from BLACS context (if already present, does nothing)
blacs_context = matrix%matrix_struct%context%get_handle()
CALL cp_dlaf_create_grid(blacs_context)
@ -282,6 +283,7 @@ CONTAINS
CALL cp_cfm_uplo_to_full(amatrix, eigenvectors)
CALL cp_cfm_uplo_to_full(bmatrix, eigenvectors)
! Create DLAF grid from BLACS context (if already present, does nothing)
blacs_context = amatrix%matrix_struct%context%get_handle()
CALL cp_dlaf_create_grid(blacs_context)

View file

@ -86,10 +86,10 @@ MODULE cp_fm_diag
! The following saved variables are diagonalization global
! Stores the default library for diagonalization
INTEGER, SAVE, PUBLIC :: diag_type = 0
INTEGER, SAVE, PUBLIC :: diag_type = 0
! Minimum number of eigenvectors for the use of the ELPA eigensolver.
! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
INTEGER, SAVE :: elpa_neigvec_min = 0
INTEGER, SAVE :: elpa_neigvec_min = 0
#if defined(__DLAF)
! Minimum number of eigenvectors for the use of the DLAF eigensolver.
! The ScaLAPACK eigensolver is used as fallback for all smaller cases.
@ -701,13 +701,13 @@ CONTAINS
INTEGER, DIMENSION(:), ALLOCATABLE :: iclustr
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: gap
INTEGER, EXTERNAL :: iceil, numroc
#if defined (__HAS_IEEE_EXCEPTIONS)
LOGICAL, DIMENSION(5) :: halt
#endif
#else
INTEGER :: nla, nlz
INTEGER, EXTERNAL :: ilaenv
#endif
#if defined (__HAS_IEEE_EXCEPTIONS)
LOGICAL, DIMENSION(5) :: halt
#endif
! by default all
n = matrix%matrix_struct%nrow_global
@ -786,22 +786,19 @@ CONTAINS
ALLOCATE (iwork(liwork))
ALLOCATE (work(lwork))
! Scalapack takes advantage of IEEE754 exceptions for speedup.
! ScaLAPACK takes advantage of IEEE754 exceptions for speedup.
! Therefore, we disable floating point traps temporarily.
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_get_halting_mode(IEEE_ALL, halt)
CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
#endif
CALL pdsyevx(job_type, "I", "U", n, a(1, 1), 1, 1, desca, vl, vu, 1, neig_local, abstol, m, nz, w(1), orfac, &
z(1, 1), 1, 1, descz, work(1), lwork, iwork(1), liwork, ifail(1), iclustr(1), gap, info)
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_set_halting_mode(IEEE_ALL, halt)
#endif
! Error handling
IF (info /= 0) THEN
IF (ionode) THEN
output_unit = cp_logger_get_unit_nr(logger, local=.FALSE.)
@ -823,12 +820,8 @@ CONTAINS
END IF
! Release work storage
DEALLOCATE (gap)
DEALLOCATE (iclustr)
DEALLOCATE (ifail)
DEALLOCATE (iwork)
DEALLOCATE (work)
#else
@ -856,11 +849,19 @@ CONTAINS
nla = SIZE(a, 1)
nlz = SIZE(z, 1)
! LAPACK takes advantage of IEEE754 exceptions for speedup.
! Therefore, we disable floating point traps temporarily.
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_get_halting_mode(IEEE_ALL, halt)
CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
#endif
CALL dsyevx(job_type, "I", "U", n, a(1, 1), nla, vl, vu, 1, neig_local, &
abstol, m, w, z(1, 1), nlz, work(1), lwork, iwork(1), ifail(1), info)
#if defined (__HAS_IEEE_EXCEPTIONS)
CALL ieee_set_halting_mode(IEEE_ALL, halt)
#endif
! Error handling
IF (info /= 0) THEN
output_unit = cp_logger_get_unit_nr(logger, local=.FALSE.)
WRITE (unit=output_unit, FMT="(/,(T3,A,T12,1X,I10))") &
@ -872,13 +873,11 @@ CONTAINS
CPABORT("Error in DSYEVX (ScaLAPACK)")
END IF
#endif
! Release work storage
DEALLOCATE (ifail)
DEALLOCATE (iwork)
DEALLOCATE (work)
#endif
eigenvalues(1:neig_local) = w(1:neig_local)
DEALLOCATE (w)

View file

@ -101,7 +101,6 @@ CONTAINS
INTEGER :: handle
CALL timeset(routineN, handle)
#if defined(__DLAF)
CALL dlaf_pspotrf(uplo, n, a, ia, ja, desca, info)
#else
@ -247,7 +246,7 @@ CONTAINS
CHARACTER, PARAMETER :: uplo = 'L'
CHARACTER(LEN=100) :: message
INTEGER :: dlaf_handle, handle, n
INTEGER :: blacs_context, dlaf_handle, handle, n
INTEGER, DIMENSION(9) :: desca, descz
INTEGER, TARGET :: info
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, z
@ -259,8 +258,9 @@ CONTAINS
! Use eigenvectors matrix as workspace
CALL cp_fm_uplo_to_full(matrix, eigenvectors)
! Create DLAF grid from BLACS context; if already present, does nothing
CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
! Create DLAF grid from BLACS context (if already present, does nothing)
blacs_context = matrix%matrix_struct%context%get_handle()
CALL cp_dlaf_create_grid(blacs_context)
n = matrix%matrix_struct%nrow_global
@ -272,9 +272,7 @@ CONTAINS
info = -1
CALL timeset(dlaf_name, dlaf_handle)
CALL dlaf_pdsyevd(uplo, n, a, 1, 1, desca, eigenvalues, z, 1, 1, descz, info)
CALL timestop(dlaf_handle)
IF (info /= 0) THEN
@ -295,6 +293,7 @@ CONTAINS
MARK_USED(dlaf_handle)
MARK_USED(dlaf_name)
MARK_USED(message)
MARK_USED(blacs_context)
CPABORT("CP2K compiled without DLA-Future-Fortran library.")
#endif
@ -357,7 +356,7 @@ CONTAINS
CHARACTER, PARAMETER :: uplo = 'L'
CHARACTER(LEN=100) :: message
INTEGER :: dlaf_handle, handle, n
INTEGER :: blacs_context, dlaf_handle, handle, n
INTEGER, DIMENSION(9) :: desca, descb, descz
INTEGER, TARGET :: info
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, z
@ -371,7 +370,8 @@ CONTAINS
CALL cp_fm_uplo_to_full(b_matrix, eigenvectors)
! Create DLAF grid from BLACS context; if already present, does nothing
CALL cp_dlaf_create_grid(a_matrix%matrix_struct%context%get_handle())
blacs_context = a_matrix%matrix_struct%context%get_handle()
CALL cp_dlaf_create_grid(blacs_context)
n = a_matrix%matrix_struct%nrow_global
@ -406,6 +406,7 @@ CONTAINS
MARK_USED(uplo)
MARK_USED(n)
MARK_USED(info)
MARK_USED(blacs_context)
MARK_USED(dlaf_handle)
MARK_USED(dlaf_name)
MARK_USED(message)