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 provides a unified interface to lapack geev routines
      10             : !> \par History
      11             : !>       2014.09 created [Florian Schiffmann]
      12             : !>       2023.12 Removed support for single-precision [Ole Schuett]
      13             : !>       2024.12 Removed support for complex input matrices [Ole Schuett]
      14             : !> \author Florian Schiffmann
      15             : ! **************************************************************************************************
      16             : MODULE arnoldi_geev
      17             :    USE kinds,                           ONLY: dp
      18             : #include "../base/base_uses.f90"
      19             : 
      20             :    IMPLICIT NONE
      21             : 
      22             :    PRIVATE
      23             : 
      24             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_geev'
      25             : 
      26             :    PUBLIC :: arnoldi_general_local_diag, arnoldi_tridiag_local_diag, arnoldi_symm_local_diag
      27             : 
      28             : CONTAINS
      29             : 
      30             : ! **************************************************************************************************
      31             : !> \brief ...
      32             : !> \param jobvr ...
      33             : !> \param matrix ...
      34             : !> \param ndim ...
      35             : !> \param evals ...
      36             : !> \param revec ...
      37             : ! **************************************************************************************************
      38        4463 :    SUBROUTINE arnoldi_symm_local_diag(jobvr, matrix, ndim, evals, revec)
      39             :       CHARACTER(1)                                       :: jobvr
      40             :       REAL(dp), DIMENSION(:, :)                          :: matrix
      41             :       INTEGER                                            :: ndim
      42             :       COMPLEX(dp), DIMENSION(:)                          :: evals
      43             :       COMPLEX(dp), DIMENSION(:, :)                       :: revec
      44             : 
      45        8926 :       INTEGER                                            :: i, info, liwork, lwork, iwork(3 + 5*ndim)
      46        8926 :       REAL(dp)                                           :: tmp_array(ndim, ndim), &
      47        8926 :                                                             work(1 + 6*ndim + 2*ndim**2)
      48        8926 :       REAL(dp), DIMENSION(ndim)                          :: eval
      49             : 
      50        4463 :       lwork = 1 + 6*ndim + 2*ndim**2
      51        4463 :       liwork = 3 + 5*ndim
      52             : 
      53     1374755 :       tmp_array(:, :) = matrix(:, :)
      54        4463 :       CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info)
      55             : 
      56       78384 :       DO i = 1, ndim
      57     1370292 :          revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, dp), dp)
      58       78384 :          evals(i) = CMPLX(eval(i), 0.0, dp)
      59             :       END DO
      60             : 
      61        4463 :    END SUBROUTINE arnoldi_symm_local_diag
      62             : 
      63             : ! **************************************************************************************************
      64             : !> \brief ...
      65             : !> \param jobvl ...
      66             : !> \param jobvr ...
      67             : !> \param matrix ...
      68             : !> \param ndim ...
      69             : !> \param evals ...
      70             : !> \param revec ...
      71             : !> \param levec ...
      72             : ! **************************************************************************************************
      73        8652 :    SUBROUTINE arnoldi_tridiag_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
      74             :       CHARACTER(1)                                       :: jobvl, jobvr
      75             :       REAL(dp), DIMENSION(:, :)                          :: matrix
      76             :       INTEGER                                            :: ndim
      77             :       COMPLEX(dp), DIMENSION(:)                          :: evals
      78             :       COMPLEX(dp), DIMENSION(:, :)                       :: revec, levec
      79             : 
      80             :       INTEGER                                            :: i, info
      81       17304 :       REAL(dp)                                           :: work(20*ndim)
      82       17304 :       REAL(dp), DIMENSION(ndim)                          :: diag, offdiag
      83       17304 :       REAL(dp), DIMENSION(ndim, ndim)                    :: evec_r
      84             : 
      85             :       MARK_USED(jobvl) !the argument has to be here for the template to work
      86             : 
      87        8652 :       levec(1, 1) = CMPLX(0.0, 0.0, dp)
      88        8652 :       info = 0
      89        8652 :       diag(ndim) = matrix(ndim, ndim)
      90       95377 :       DO i = 1, ndim - 1
      91       86725 :          diag(i) = matrix(i, i)
      92       95377 :          offdiag(i) = matrix(i + 1, i)
      93             : 
      94             :       END DO
      95             : 
      96        8652 :       CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info)
      97             : 
      98      104029 :       DO i = 1, ndim
      99     2214198 :          revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
     100      104029 :          evals(i) = CMPLX(diag(i), 0.0, dp)
     101             :       END DO
     102        8652 :    END SUBROUTINE arnoldi_tridiag_local_diag
     103             : 
     104             : ! **************************************************************************************************
     105             : !> \brief ...
     106             : !> \param jobvl ...
     107             : !> \param jobvr ...
     108             : !> \param matrix ...
     109             : !> \param ndim ...
     110             : !> \param evals ...
     111             : !> \param revec ...
     112             : !> \param levec ...
     113             : ! **************************************************************************************************
     114      116673 :    SUBROUTINE arnoldi_general_local_diag(jobvl, jobvr, matrix, ndim, evals, revec, levec)
     115             :       CHARACTER(1)                                       :: jobvl, jobvr
     116             :       REAL(dp), DIMENSION(:, :)                          :: matrix
     117             :       INTEGER                                            :: ndim
     118             :       COMPLEX(dp), DIMENSION(:)                          :: evals
     119             :       COMPLEX(dp), DIMENSION(:, :)                       :: revec, levec
     120             : 
     121             :       INTEGER                                            :: i, info, lwork
     122      233346 :       LOGICAL                                            :: selects(ndim)
     123      233346 :       REAL(dp)                                           :: norm, tmp_array(ndim, ndim), &
     124      233346 :                                                             work(20*ndim)
     125      233346 :       REAL(dp), DIMENSION(ndim)                          :: eval1, eval2
     126      116673 :       REAL(dp), DIMENSION(ndim, ndim)                    :: evec_l, evec_r
     127             : 
     128             :       MARK_USED(jobvr) !the argument has to be here for the template to work
     129             :       MARK_USED(jobvl) !the argument has to be here for the template to work
     130             : 
     131     1115087 :       eval1 = REAL(0.0, dp); eval2 = REAL(0.0, dp)
     132     6134419 :       tmp_array(:, :) = matrix(:, :)
     133             :       ! ask lapack how much space it would like in the work vector, don't ask me why
     134      116673 :       lwork = -1
     135      116673 :       CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
     136             : 
     137      116673 :       lwork = MIN(20*ndim, INT(work(1)))
     138      116673 :       CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info)
     139      116673 :       CALL dtrevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info)
     140             : 
     141             :       ! compose the eigenvectors, lapacks way of storing them is a pain
     142             :       ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec
     143             :       ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized
     144      116673 :       i = 1
     145      615880 :       DO WHILE (i .LE. ndim)
     146      615880 :          IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, dp))) THEN
     147    11536285 :             evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i)))
     148     6017746 :             revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, dp), dp)
     149     6017746 :             levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, dp), dp)
     150      499207 :             i = i + 1
     151           0 :          ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, dp))) THEN
     152           0 :             norm = SQRT(SUM(evec_r(:, i)**2.0_dp) + SUM(evec_r(:, i + 1)**2.0_dp))
     153           0 :             revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), dp)/norm
     154           0 :             revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), dp)/norm
     155           0 :             levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), dp)
     156           0 :             levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), dp)
     157           0 :             i = i + 2
     158             :          ELSE
     159           0 :             CPABORT('something went wrong while sorting the EV in arnoldi_geev')
     160             :          END IF
     161             :       END DO
     162             : 
     163             :       ! this is to keep the interface consistent with complex geev
     164      615880 :       DO i = 1, ndim
     165      615880 :          evals(i) = CMPLX(eval1(i), eval2(i), dp)
     166             :       END DO
     167             : 
     168      116673 :    END SUBROUTINE arnoldi_general_local_diag
     169             : 
     170             : END MODULE arnoldi_geev