LCOV - code coverage report
Current view: top level - src/fm - cp_cfm_diag.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.5 % 85 82
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief used for collecting diagonalization schemes available for cp_cfm_type
      10              : !> \note
      11              : !>      first version : only one routine right now
      12              : !> \author Joost VandeVondele (2003-09)
      13              : ! **************************************************************************************************
      14              : MODULE cp_cfm_diag
      15              :    USE cp_cfm_cholesky, ONLY: cp_cfm_cholesky_decompose
      16              :    USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm, &
      17              :                                   cp_cfm_column_scale, &
      18              :                                   cp_cfm_scale, &
      19              :                                   cp_cfm_triangular_invert, &
      20              :                                   cp_cfm_triangular_multiply
      21              :    USE cp_cfm_types, ONLY: cp_cfm_get_info, &
      22              :                            cp_cfm_set_element, &
      23              :                            cp_cfm_to_cfm, &
      24              :                            cp_cfm_type
      25              : #if defined(__DLAF)
      26              :    USE cp_cfm_dlaf_api, ONLY: cp_cfm_diag_gen_dlaf, &
      27              :                               cp_cfm_diag_dlaf
      28              :    USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
      29              :    USE cp_fm_diag, ONLY: diag_type, dlaf_neigvec_min, FM_DIAG_TYPE_DLAF
      30              : #endif
      31              :    USE kinds, ONLY: dp
      32              : #if defined (__HAS_IEEE_EXCEPTIONS)
      33              :    USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
      34              :                               ieee_set_halting_mode, &
      35              :                               IEEE_ALL
      36              : #endif
      37              : 
      38              : #include "../base/base_uses.f90"
      39              : 
      40              :    IMPLICIT NONE
      41              :    PRIVATE
      42              : 
      43              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_diag'
      44              : 
      45              :    PUBLIC :: cp_cfm_heevd, cp_cfm_geeig, cp_cfm_geeig_canon
      46              : 
      47              : CONTAINS
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief Perform a diagonalisation of a complex matrix
      51              : !> \param matrix ...
      52              : !> \param eigenvectors ...
      53              : !> \param eigenvalues ...
      54              : !> \par History
      55              : !>      12.2024 Added DLA-Future support [Rocco Meli]
      56              : !> \author Joost VandeVondele
      57              : ! **************************************************************************************************
      58        26256 :    SUBROUTINE cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
      59              : 
      60              :       TYPE(cp_cfm_type), INTENT(IN)                      :: matrix, eigenvectors
      61              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigenvalues
      62              : 
      63              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_cfm_heevd'
      64              : 
      65              :       INTEGER                                            :: handle
      66              : 
      67        26256 :       CALL timeset(routineN, handle)
      68              : 
      69              : #if defined(__DLAF)
      70              :       IF (diag_type == FM_DIAG_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
      71              :          ! Initialize DLA-Future on-demand; if already initialized, does nothing
      72              :          CALL cp_dlaf_initialize()
      73              : 
      74              :          ! Create DLAF grid from BLACS context; if already present, does nothing
      75              :          CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
      76              : 
      77              :          CALL cp_cfm_diag_dlaf(matrix, eigenvectors, eigenvalues)
      78              :       ELSE
      79              : #endif
      80        26256 :          CALL cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
      81              : #if defined(__DLAF)
      82              :       END IF
      83              : #endif
      84              : 
      85        26256 :       CALL timestop(handle)
      86              : 
      87        26256 :    END SUBROUTINE cp_cfm_heevd
      88              : 
      89              : ! **************************************************************************************************
      90              : !> \brief Perform a diagonalisation of a complex matrix
      91              : !> \param matrix ...
      92              : !> \param eigenvectors ...
      93              : !> \param eigenvalues ...
      94              : !> \par History
      95              : !>      - (De)Allocation checks updated (15.02.2011,MK)
      96              : !> \author Joost VandeVondele
      97              : ! **************************************************************************************************
      98        26256 :    SUBROUTINE cp_cfm_heevd_base(matrix, eigenvectors, eigenvalues)
      99              : 
     100              :       TYPE(cp_cfm_type), INTENT(IN)            :: matrix, eigenvectors
     101              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
     102              : 
     103              :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_heevd_base'
     104              : 
     105        26256 :       COMPLEX(KIND=dp), DIMENSION(:), POINTER  :: work
     106              :       COMPLEX(KIND=dp), DIMENSION(:, :), &
     107        26256 :          POINTER                               :: m
     108              :       INTEGER                                  :: handle, info, liwork, &
     109              :                                                   lrwork, lwork, n
     110        26256 :       INTEGER, DIMENSION(:), POINTER           :: iwork
     111        26256 :       REAL(KIND=dp), DIMENSION(:), POINTER     :: rwork
     112              : #if defined(__parallel)
     113              :       INTEGER, DIMENSION(9)                    :: descm, descv
     114              :       COMPLEX(KIND=dp), DIMENSION(:, :), &
     115        26256 :          POINTER                               :: v
     116              : #if defined (__HAS_IEEE_EXCEPTIONS)
     117              :       LOGICAL, DIMENSION(5)                    :: halt
     118              : #endif
     119              : #endif
     120              : 
     121        26256 :       CALL timeset(routineN, handle)
     122              : 
     123        26256 :       n = matrix%matrix_struct%nrow_global
     124        26256 :       m => matrix%local_data
     125        26256 :       ALLOCATE (iwork(1), rwork(1), work(1))
     126              :       ! work space query
     127        26256 :       lwork = -1
     128        26256 :       lrwork = -1
     129        26256 :       liwork = -1
     130              : 
     131              : #if defined(__parallel)
     132        26256 :       v => eigenvectors%local_data
     133       262560 :       descm(:) = matrix%matrix_struct%descriptor(:)
     134       262560 :       descv(:) = eigenvectors%matrix_struct%descriptor(:)
     135              :       CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
     136        26256 :                    work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
     137              :       ! The work space query for lwork does not return always sufficiently large values.
     138              :       ! Let's add some margin to avoid crashes.
     139        26256 :       lwork = CEILING(REAL(work(1), KIND=dp)) + 1000
     140              :       ! needed to correct for a bug in scalapack, unclear how much the right number is
     141        26256 :       lrwork = CEILING(rwork(1)) + 1000000
     142        26256 :       liwork = iwork(1)
     143              : #else
     144              :       CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
     145              :                   work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
     146              :       lwork = CEILING(REAL(work(1), KIND=dp))
     147              :       lrwork = CEILING(rwork(1))
     148              :       liwork = iwork(1)
     149              : #endif
     150              : 
     151        26256 :       DEALLOCATE (iwork, rwork, work)
     152       183792 :       ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
     153              : 
     154              : #if defined(__parallel)
     155              : ! Scalapack takes advantage of IEEE754 exceptions for speedup.
     156              : ! Therefore, we disable floating point traps temporarily.
     157              : #if defined (__HAS_IEEE_EXCEPTIONS)
     158              :       CALL ieee_get_halting_mode(IEEE_ALL, halt)
     159              :       CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
     160              : #endif
     161              : 
     162              :       CALL pzheevd('V', 'U', n, m(1, 1), 1, 1, descm, eigenvalues(1), v(1, 1), 1, 1, descv, &
     163        26256 :                    work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
     164              : 
     165              : #if defined (__HAS_IEEE_EXCEPTIONS)
     166              :       CALL ieee_set_halting_mode(IEEE_ALL, halt)
     167              : #endif
     168              : #else
     169              :       CALL zheevd('V', 'U', n, m(1, 1), SIZE(m, 1), eigenvalues(1), &
     170              :                   work(1), lwork, rwork(1), lrwork, iwork(1), liwork, info)
     171              :       eigenvectors%local_data = matrix%local_data
     172              : #endif
     173              : 
     174        26256 :       DEALLOCATE (iwork, rwork, work)
     175        26256 :       IF (info /= 0) &
     176            0 :          CPABORT("Diagonalisation of a complex matrix failed")
     177              : 
     178        26256 :       CALL timestop(handle)
     179              : 
     180        26256 :    END SUBROUTINE cp_cfm_heevd_base
     181              : 
     182              : ! **************************************************************************************************
     183              : !> \brief General Eigenvalue Problem  AX = BXE
     184              : !>        Single option version: Cholesky decomposition of B
     185              : !> \param amatrix ...
     186              : !> \param bmatrix ...
     187              : !> \param eigenvectors ...
     188              : !> \param eigenvalues ...
     189              : !> \param work ...
     190              : !> \par History
     191              : !>      12.2024 Added DLA-Future support [Rocco Meli]
     192              : ! **************************************************************************************************
     193        16710 :    SUBROUTINE cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
     194              : 
     195              :       TYPE(cp_cfm_type), INTENT(IN)                      :: amatrix, bmatrix, eigenvectors
     196              :       REAL(KIND=dp), DIMENSION(:)                        :: eigenvalues
     197              :       TYPE(cp_cfm_type), INTENT(IN)                      :: work
     198              : 
     199              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_cfm_geeig'
     200              : 
     201              :       INTEGER                                            :: handle, nao, nmo
     202              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     203              : 
     204        16710 :       CALL timeset(routineN, handle)
     205              : 
     206        16710 :       CALL cp_cfm_get_info(amatrix, nrow_global=nao)
     207        50130 :       ALLOCATE (evals(nao))
     208        16710 :       nmo = SIZE(eigenvalues)
     209              : 
     210              : #if defined(__DLAF)
     211              :       IF (diag_type == FM_DIAG_TYPE_DLAF .AND. amatrix%matrix_struct%nrow_global >= dlaf_neigvec_min) THEN
     212              :          ! Initialize DLA-Future on-demand; if already initialized, does nothing
     213              :          CALL cp_dlaf_initialize()
     214              : 
     215              :          ! Create DLAF grid from BLACS context; if already present, does nothing
     216              :          CALL cp_dlaf_create_grid(amatrix%matrix_struct%context%get_handle())
     217              :          CALL cp_dlaf_create_grid(bmatrix%matrix_struct%context%get_handle())
     218              :          CALL cp_dlaf_create_grid(eigenvectors%matrix_struct%context%get_handle())
     219              : 
     220              :          ! Use DLA-Future generalized eigenvalue solver for large matrices
     221              :          CALL cp_cfm_diag_gen_dlaf(amatrix, bmatrix, work, evals)
     222              :       ELSE
     223              : #endif
     224              :          ! Cholesky decompose S=U(T)U
     225        16710 :          CALL cp_cfm_cholesky_decompose(bmatrix)
     226              :          ! Invert to get U^(-1)
     227        16710 :          CALL cp_cfm_triangular_invert(bmatrix)
     228              :          ! Reduce to get U^(-T) * H * U^(-1)
     229        16710 :          CALL cp_cfm_triangular_multiply(bmatrix, amatrix, side="R")
     230        16710 :          CALL cp_cfm_triangular_multiply(bmatrix, amatrix, transa_tr="C")
     231              :          ! Diagonalize
     232        16710 :          CALL cp_cfm_heevd(matrix=amatrix, eigenvectors=work, eigenvalues=evals)
     233              :          ! Restore vectors C = U^(-1) * C*
     234        16710 :          CALL cp_cfm_triangular_multiply(bmatrix, work)
     235              : #if defined(__DLAF)
     236              :       END IF
     237              : #endif
     238              : 
     239        16710 :       CALL cp_cfm_to_cfm(work, eigenvectors, nmo)
     240       333414 :       eigenvalues(1:nmo) = evals(1:nmo)
     241              : 
     242        16710 :       DEALLOCATE (evals)
     243              : 
     244        16710 :       CALL timestop(handle)
     245              : 
     246        16710 :    END SUBROUTINE cp_cfm_geeig
     247              : 
     248              : ! **************************************************************************************************
     249              : !> \brief General Eigenvalue Problem  AX = BXE
     250              : !>        Use canonical orthogonalization
     251              : !> \param amatrix ...
     252              : !> \param bmatrix ...
     253              : !> \param eigenvectors ...
     254              : !> \param eigenvalues ...
     255              : !> \param work ...
     256              : !> \param epseig ...
     257              : ! **************************************************************************************************
     258          432 :    SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
     259              : 
     260              :       TYPE(cp_cfm_type), INTENT(IN)                      :: amatrix, bmatrix, eigenvectors
     261              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigenvalues
     262              :       TYPE(cp_cfm_type), INTENT(IN)                      :: work
     263              :       REAL(KIND=dp), INTENT(IN)                          :: epseig
     264              : 
     265              :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_geeig_canon'
     266              :       COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
     267              :          czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
     268              : 
     269              :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: cevals
     270              :       INTEGER                                            :: handle, i, icol, irow, nao, nc, ncol, &
     271              :                                                             nmo, nx
     272              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
     273              : 
     274          432 :       CALL timeset(routineN, handle)
     275              : 
     276              :       ! Test sizes
     277          432 :       CALL cp_cfm_get_info(amatrix, nrow_global=nao)
     278          432 :       nmo = SIZE(eigenvalues)
     279         2160 :       ALLOCATE (evals(nao), cevals(nao))
     280              : 
     281              :       ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
     282          432 :       CALL cp_cfm_scale(-cone, bmatrix)
     283          432 :       CALL cp_cfm_heevd(bmatrix, work, evals)
     284        19502 :       evals(:) = -evals(:)
     285          432 :       nc = nao
     286        18938 :       DO i = 1, nao
     287        18938 :          IF (evals(i) < epseig) THEN
     288           72 :             nc = i - 1
     289           72 :             EXIT
     290              :          END IF
     291              :       END DO
     292          432 :       CPASSERT(nc /= 0)
     293              : 
     294          432 :       IF (nc /= nao) THEN
     295           72 :          IF (nc < nmo) THEN
     296              :             ! Copy NULL space definition to last vectors of eigenvectors (if needed)
     297            0 :             ncol = nmo - nc
     298            0 :             CALL cp_cfm_to_cfm(work, eigenvectors, ncol, nc + 1, nc + 1)
     299              :          END IF
     300              :          ! Set NULL space in eigenvector matrix of S to zero
     301          636 :          DO icol = nc + 1, nao
     302        80028 :             DO irow = 1, nao
     303        79956 :                CALL cp_cfm_set_element(work, irow, icol, czero)
     304              :             END DO
     305              :          END DO
     306              :          ! Set small eigenvalues to a dummy save value
     307          636 :          evals(nc + 1:nao) = 1.0_dp
     308              :       END IF
     309              :       ! Calculate U*s**(-1/2)
     310        19502 :       cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
     311          432 :       CALL cp_cfm_column_scale(work, cevals)
     312              :       ! Reduce to get U^(-C) * H * U^(-1)
     313          432 :       CALL cp_cfm_gemm("C", "N", nao, nao, nao, cone, work, amatrix, czero, bmatrix)
     314          432 :       CALL cp_cfm_gemm("N", "N", nao, nao, nao, cone, bmatrix, work, czero, amatrix)
     315          432 :       IF (nc /= nao) THEN
     316              :          ! set diagonal values to save large value
     317          636 :          DO icol = nc + 1, nao
     318          636 :             CALL cp_cfm_set_element(amatrix, icol, icol, CMPLX(10000.0_dp, 0.0_dp, KIND=dp))
     319              :          END DO
     320              :       END IF
     321              :       ! Diagonalize
     322          432 :       CALL cp_cfm_heevd(amatrix, bmatrix, evals)
     323         6678 :       eigenvalues(1:nmo) = evals(1:nmo)
     324          432 :       nx = MIN(nc, nmo)
     325              :       ! Restore vectors C = U^(-1) * C*
     326          432 :       CALL cp_cfm_gemm("N", "N", nao, nx, nc, cone, work, bmatrix, czero, eigenvectors)
     327              : 
     328          432 :       DEALLOCATE (evals)
     329              : 
     330          432 :       CALL timestop(handle)
     331              : 
     332          864 :    END SUBROUTINE cp_cfm_geeig_canon
     333              : 
     334              : END MODULE cp_cfm_diag
        

Generated by: LCOV version 2.0-1