Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Provides interfaces to LAPACK eigenvalue/SVD routines
      10             : !> \note
      11             : !>      We are using LAPACK interfaces, so please make sure in IBM/AIX you have
      12             : !>      the lapack library before essl: "xlf90 ... -llapack -lessl" !!!
      13             : !> \par History
      14             : !>      JGH (26-5-2001): delay D/S C/Z problem to the lapack library call
      15             : !> \author APSI
      16             : ! **************************************************************************************************
      17             : MODULE eigenvalueproblems
      18             : 
      19             :    USE kinds,                           ONLY: dp
      20             :    USE lapack,                          ONLY: lapack_cgesvd,&
      21             :                                               lapack_chpev,&
      22             :                                               lapack_sgesvd,&
      23             :                                               lapack_ssyev
      24             : #include "../base/base_uses.f90"
      25             : 
      26             :    IMPLICIT NONE
      27             : 
      28             :    PRIVATE
      29             : 
      30             :    PUBLIC :: diagonalise
      31             : 
      32             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eigenvalueproblems'
      33             : 
      34             :    INTERFACE diagonalise
      35             :       MODULE PROCEDURE diagonalise_ssyev
      36             :       MODULE PROCEDURE diagonalise_chpev
      37             :    END INTERFACE
      38             : 
      39             :    INTERFACE singular_values
      40             :       MODULE PROCEDURE cp2k_sgesvd
      41             :       MODULE PROCEDURE cp2k_cgesvd
      42             :    END INTERFACE
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief ...
      48             : !> \param matrix ...
      49             : !> \param mysize ...
      50             : !> \param storageform ...
      51             : !> \param eigenvalues ...
      52             : !> \param eigenvectors ...
      53             : ! **************************************************************************************************
      54         988 :    SUBROUTINE diagonalise_ssyev(matrix, mysize, storageform, eigenvalues, &
      55         988 :                                 eigenvectors)
      56             : 
      57             :       REAL(KIND=dp), INTENT(IN)                          :: matrix(:, :)
      58             :       INTEGER, INTENT(IN)                                :: mysize
      59             :       CHARACTER(LEN=*), INTENT(IN)                       :: storageform
      60             :       REAL(KIND=dp), INTENT(OUT)                         :: eigenvalues(:), eigenvectors(:, :)
      61             : 
      62             :       CHARACTER, PARAMETER                               :: jobz = "V"
      63             : 
      64             :       CHARACTER                                          :: uplo
      65             :       INTEGER                                            :: info, lda, lwork
      66        1976 :       REAL(KIND=dp)                                      :: work(3*mysize - 1)
      67             : 
      68             :       IF (storageform(1:5) == "Lower" .OR. &
      69         988 :           storageform(1:5) == "LOWER" .OR. &
      70             :           storageform(1:5) == "lower") THEN
      71           0 :          uplo = "L"
      72             :       ELSE IF (storageform(1:5) == "Upper" .OR. &
      73         988 :                storageform(1:5) == "upper" .OR. &
      74             :                storageform(1:5) == "UPPER") THEN
      75         988 :          uplo = "U"
      76             :       ELSE
      77           0 :          CPABORT("Unknown form of storage")
      78             :       END IF
      79             : 
      80         988 :       lda = SIZE(matrix, 1)
      81         988 :       lwork = 3*mysize - 1
      82             : 
      83       12844 :       eigenvectors = matrix
      84             : 
      85             :       CALL lapack_ssyev(jobz, uplo, mysize, eigenvectors, lda, eigenvalues, &
      86         988 :                         work, lwork, info)
      87         988 :       IF (info /= 0) THEN
      88           0 :          CPABORT("Error in diagonalisation")
      89             :       END IF
      90             : 
      91         988 :    END SUBROUTINE diagonalise_ssyev
      92             : 
      93             : ! **************************************************************************************************
      94             : !> \brief ...
      95             : !> \param matrix ...
      96             : !> \param mysize ...
      97             : !> \param storageform ...
      98             : !> \param eigenvalues ...
      99             : !> \param eigenvectors ...
     100             : ! **************************************************************************************************
     101           0 :    SUBROUTINE diagonalise_chpev(matrix, mysize, storageform, eigenvalues, &
     102           0 :                                 eigenvectors)
     103             : 
     104             :       COMPLEX(KIND=dp), INTENT(INOUT)                    :: matrix(:)
     105             :       INTEGER, INTENT(IN)                                :: mysize
     106             :       CHARACTER(LEN=*), INTENT(IN)                       :: storageform
     107             :       REAL(KIND=dp), INTENT(OUT)                         :: eigenvalues(:)
     108             :       COMPLEX(KIND=dp), INTENT(OUT)                      :: eigenvectors(:, :)
     109             : 
     110             :       CHARACTER, PARAMETER                               :: jobz = "V"
     111             : 
     112             :       CHARACTER                                          :: uplo
     113             :       INTEGER                                            :: info
     114           0 :       COMPLEX(KIND=dp)                                   :: work(2*mysize - 1)
     115           0 :       REAL(KIND=dp)                                      :: rwork(3*mysize - 2)
     116             : 
     117             :       IF (storageform(1:5) == "Lower" .OR. &
     118           0 :           storageform(1:5) == "LOWER" .OR. &
     119             :           storageform(1:5) == "lower") THEN
     120           0 :          uplo = "L"
     121             :       ELSE IF (storageform(1:5) == "Upper" .OR. &
     122           0 :                storageform(1:5) == "upper" .OR. &
     123             :                storageform(1:5) == "UPPER") THEN
     124           0 :          uplo = "U"
     125             :       ELSE
     126           0 :          CPABORT("Unknown form of storage")
     127             :       END IF
     128             : 
     129             :       CALL lapack_chpev(jobz, uplo, mysize, matrix, eigenvalues, &
     130           0 :                         eigenvectors, mysize, work, rwork, info)
     131           0 :       IF (info /= 0) THEN
     132           0 :          CPABORT("Error in diagonalisation")
     133             :       END IF
     134             : 
     135           0 :    END SUBROUTINE diagonalise_chpev
     136             : 
     137             : ! **************************************************************************************************
     138             : !> \brief ...
     139             : !> \param matrix ...
     140             : !> \param svalues ...
     141             : !> \param mrow ...
     142             : !> \param ncol ...
     143             : !> \param uvec ...
     144             : !> \param vtvec ...
     145             : ! **************************************************************************************************
     146           0 :    SUBROUTINE cp2k_sgesvd(matrix, svalues, mrow, ncol, uvec, vtvec)
     147             : 
     148             :       REAL(KIND=dp), INTENT(IN)                          :: matrix(:, :)
     149             :       REAL(KIND=dp), INTENT(OUT)                         :: svalues(:)
     150             :       INTEGER, INTENT(IN)                                :: mrow, ncol
     151             :       REAL(KIND=dp), INTENT(OUT)                         :: uvec(:, :), vtvec(:, :)
     152             : 
     153             :       CHARACTER, PARAMETER                               :: jobu = "A", jobvt = "A"
     154             : 
     155             :       INTEGER                                            :: info, lda, ldu, ldvt, lwork
     156           0 :       REAL(KIND=dp)                                      :: work(25*(mrow + ncol))
     157             : 
     158           0 :       lwork = 25*(mrow + ncol)
     159           0 :       lda = SIZE(matrix, 1)
     160           0 :       ldu = SIZE(uvec, 1)
     161           0 :       ldvt = SIZE(vtvec, 1)
     162             : 
     163             :       CALL lapack_sgesvd(jobu, jobvt, mrow, ncol, matrix, lda, svalues, &
     164           0 :                          uvec, ldu, vtvec, ldvt, work, lwork, info)
     165           0 :       IF (info /= 0) THEN
     166           0 :          CPABORT("Error in singular value decomposition.")
     167             :       END IF
     168             : 
     169           0 :    END SUBROUTINE cp2k_sgesvd
     170             : 
     171             : ! **************************************************************************************************
     172             : !> \brief ...
     173             : !> \param matrix ...
     174             : !> \param svalues ...
     175             : !> \param mrow ...
     176             : !> \param ncol ...
     177             : !> \param uvec ...
     178             : !> \param vtvec ...
     179             : ! **************************************************************************************************
     180           0 :    SUBROUTINE cp2k_cgesvd(matrix, svalues, mrow, ncol, uvec, vtvec)
     181             : 
     182             :       COMPLEX(KIND=dp), INTENT(IN)                       :: matrix(:, :)
     183             :       REAL(KIND=dp), INTENT(OUT)                         :: svalues(:)
     184             :       INTEGER, INTENT(IN)                                :: mrow, ncol
     185             :       COMPLEX(KIND=dp), INTENT(OUT)                      :: uvec(:, :), vtvec(:, :)
     186             : 
     187             :       CHARACTER, PARAMETER                               :: jobu = "A", jobvt = "A"
     188             : 
     189             :       INTEGER                                            :: info, lda, ldu, ldvt, lwork
     190           0 :       COMPLEX(KIND=dp)                                   :: work(25*(mrow + ncol))
     191           0 :       REAL(KIND=dp)                                      :: rwork(25*(mrow + ncol))
     192             : 
     193           0 :       lwork = 25*(mrow + ncol)
     194           0 :       lda = SIZE(matrix, 1)
     195           0 :       ldu = SIZE(uvec, 1)
     196           0 :       ldvt = SIZE(vtvec, 1)
     197             : 
     198             :       CALL lapack_cgesvd(jobu, jobvt, mrow, ncol, matrix, lda, svalues, &
     199           0 :                          uvec, ldu, vtvec, ldvt, work, lwork, rwork, info)
     200           0 :       IF (info /= 0) THEN
     201           0 :          CPABORT("Error in singular value decomposition.")
     202             :       END IF
     203             : 
     204           0 :    END SUBROUTINE cp2k_cgesvd
     205             : 
     206             : END MODULE eigenvalueproblems
     207             :