ELPA: avoid multiple initializations (cp2k#3786)

- Note: potentially ensuring the original "active device" is down to the library (ELPA),
        and CP2K is only responsible for toggling the active device for its own code.
        This note refers to a potential side-effect of (re-)initializing ELPA.
- Avoid potentially high cost of initialization (context creation, etc.)
- Initialization SAVE state does not require atomics (expected to be called from master).
- Harmonized CFM and FM interfaces.
- Code/format cleanup.

src/fm/cp_cfm_diag.F

@@ -26,13 +26,14 @@ MODULE cp_cfm_diag
 #if defined (__HAS_IEEE_EXCEPTIONS)
    USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
                               ieee_set_halting_mode, &
-                              ieee_all
+                              IEEE_ALL
 #endif
 
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
    PRIVATE
+
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
 
    PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
@@ -190,7 +191,7 @@ END SUBROUTINE cp_cfm_geeig
    SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
 
       TYPE(cp_cfm_type), INTENT(IN)                      :: amatrix, bmatrix, eigenvectors
-      REAL(KIND=dp), DIMENSION(:)                        :: eigenvalues
+      REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigenvalues
       TYPE(cp_cfm_type), INTENT(IN)                      :: work
       REAL(KIND=dp), INTENT(IN)                          :: epseig
 
@@ -238,7 +239,7 @@ SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work,
          ! Set small eigenvalues to a dummy save value
          evals(nc + 1:nao) = 1.0_dp
       END IF
-      ! calculate U*s**(-1/2)
+      ! Calculate U*s**(-1/2)
       cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
       CALL cp_cfm_column_scale(work, cevals)
       ! Reduce to get U^(-C) * H * U^(-1)

src/fm/cp_fm_diag.F

@@ -155,6 +155,8 @@ SUBROUTINE diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_neigvec_min_i
       INTEGER, INTENT(IN)                                :: dlaf_neigvec_min_input
       REAL(KIND=dp), INTENT(IN)                          :: eps_check_diag_input
 
+      LOGICAL, SAVE                                      :: initialized = .FALSE.
+
       fallback_applied = .FALSE.
 
       IF (diag_lib == "ScaLAPACK") THEN
@@ -180,12 +182,13 @@ SUBROUTINE diag_init(diag_lib, fallback_applied, elpa_kernel, elpa_neigvec_min_i
          CPABORT("ERROR in diag_init: Initialization of unknown diagonalization library requested")
       END IF
 
-      ! Initialization of requested diagonalization library:
-      IF (diag_type == FM_DIAG_TYPE_ELPA) THEN
+      ! Initialization of requested diagonalization library
+      IF (.NOT. initialized .AND. diag_type == FM_DIAG_TYPE_ELPA) THEN
          CALL initialize_elpa_library()
          CALL set_elpa_kernel(elpa_kernel)
          CALL set_elpa_qr(elpa_qr, elpa_qr_unsafe)
          CALL set_elpa_print(elpa_print)
+         initialized = .TRUE.
       END IF
 
       elpa_neigvec_min = elpa_neigvec_min_input
@@ -1449,22 +1452,22 @@ SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work,
 
       INTEGER                                            :: handle, i, icol, irow, nao, nc, ncol, &
                                                             nmo, nx
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: seigval
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
 
       CALL timeset(routineN, handle)
 
       ! Test sizees
       CALL cp_fm_get_info(amatrix, nrow_global=nao)
       nmo = SIZE(eigenvalues)
-      ALLOCATE (seigval(nao))
+      ALLOCATE (evals(nao))
 
       ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
       CALL cp_fm_scale(-1.0_dp, bmatrix)
-      CALL choose_eigv_solver(matrix=bmatrix, eigenvectors=work, eigenvalues=seigval)
-      seigval(:) = -seigval(:)
+      CALL choose_eigv_solver(matrix=bmatrix, eigenvectors=work, eigenvalues=evals)
+      evals(:) = -evals(:)
       nc = nao
       DO i = 1, nao
-         IF (seigval(i) < epseig) THEN
+         IF (evals(i) < epseig) THEN
             nc = i - 1
             EXIT
          END IF
@@ -1484,11 +1487,11 @@ SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work,
             END DO
          END DO
          ! Set small eigenvalues to a dummy save value
-         seigval(nc + 1:nao) = 1.0_dp
+         evals(nc + 1:nao) = 1.0_dp
       END IF
       ! Calculate U*s**(-1/2)
-      seigval(:) = 1.0_dp/SQRT(seigval(:))
-      CALL cp_fm_column_scale(work, seigval)
+      evals(:) = 1.0_dp/SQRT(evals(:))
+      CALL cp_fm_column_scale(work, evals)
       ! Reduce to get U^(-T) * H * U^(-1)
       CALL cp_fm_gemm("T", "N", nao, nao, nao, 1.0_dp, work, amatrix, 0.0_dp, bmatrix)
       CALL cp_fm_gemm("N", "N", nao, nao, nao, 1.0_dp, bmatrix, work, 0.0_dp, amatrix)
@@ -1504,7 +1507,7 @@ SUBROUTINE cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work,
       ! Restore vectors C = U^(-1) * C*
       CALL cp_fm_gemm("N", "N", nao, nx, nc, 1.0_dp, work, bmatrix, 0.0_dp, eigenvectors)
 
-      DEALLOCATE (seigval)
+      DEALLOCATE (evals)
 
       CALL timestop(handle)