LCOV - code coverage report
Current view: top level - src/common - mathlib.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 473 608 77.8 %
Date: 2025-02-18 08:24:35 Functions: 37 41 90.2 %

          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 Collection of simple mathematical functions and subroutines
      10             : !> \par History
      11             : !>      FUNCTION angle updated and FUNCTION dihedral angle added; cleaned
      12             : !>      (13.03.2004,MK)
      13             : !> \author MK (15.11.1998)
      14             : ! **************************************************************************************************
      15             : MODULE mathlib
      16             : 
      17             :    USE kinds,                           ONLY: default_string_length,&
      18             :                                               dp
      19             :    USE mathconstants,                   ONLY: euler,&
      20             :                                               fac,&
      21             :                                               oorootpi
      22             : #include "../base/base_uses.f90"
      23             : 
      24             :    IMPLICIT NONE
      25             : 
      26             :    PRIVATE
      27             : 
      28             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mathlib'
      29             :    REAL(KIND=dp), PARAMETER             :: eps_geo = 1.0E-6_dp
      30             : 
      31             :    ! Public subroutines
      32             : 
      33             :    PUBLIC :: build_rotmat, &
      34             :              jacobi, &
      35             :              diamat_all, &
      36             :              invmat, &
      37             :              invmat_symm, &
      38             :              invert_matrix, &
      39             :              get_pseudo_inverse_svd, &
      40             :              get_pseudo_inverse_diag, &
      41             :              symmetrize_matrix, &
      42             :              unit_matrix, diag, &
      43             :              erfc_cutoff, &
      44             :              diag_antisym, &
      45             :              diag_complex
      46             : 
      47             :    ! Public functions
      48             : 
      49             :    PUBLIC :: angle, &
      50             :              binomial, &
      51             :              binomial_gen, &
      52             :              multinomial, &
      53             :              det_3x3, &
      54             :              dihedral_angle, &
      55             :              gcd, &
      56             :              inv_3x3, &
      57             :              lcm, &
      58             :              vector_product, &
      59             :              pswitch, &
      60             :              rotate_vector, &
      61             :              reflect_vector, &
      62             :              expint, abnormal_value, &
      63             :              get_diag, &
      64             :              set_diag
      65             : 
      66             :    INTERFACE det_3x3
      67             :       MODULE PROCEDURE det_3x3_1, det_3x3_2
      68             :    END INTERFACE
      69             : 
      70             :    INTERFACE invert_matrix
      71             :       MODULE PROCEDURE invert_matrix_d, invert_matrix_z
      72             :    END INTERFACE
      73             : 
      74             :    INTERFACE set_diag
      75             :       MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
      76             :    END INTERFACE
      77             : 
      78             :    INTERFACE swap
      79             :       MODULE PROCEDURE swap_scalar, swap_vector
      80             :    END INTERFACE
      81             : 
      82             :    INTERFACE unit_matrix
      83             :       MODULE PROCEDURE unit_matrix_d, unit_matrix_z
      84             :    END INTERFACE
      85             : 
      86             : CONTAINS
      87             : 
      88             : ! **************************************************************************************************
      89             : !> \brief Polynomial (5th degree) switching function
      90             : !>        f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) = 0
      91             : !> \param x ...
      92             : !> \param a ...
      93             : !> \param b ...
      94             : !> \param order ...
      95             : !> \return =0 : f(x)
      96             : !> \return =1 : f'(x)
      97             : !> \return =2 : f"(x)
      98             : ! **************************************************************************************************
      99       16030 :    FUNCTION pswitch(x, a, b, order) RESULT(fx)
     100             :       REAL(KIND=dp)                                      :: x, a, b
     101             :       INTEGER                                            :: order
     102             :       REAL(KIND=dp)                                      :: fx
     103             : 
     104             :       REAL(KIND=dp)                                      :: u, u2, u3
     105             : 
     106       16030 :       CPASSERT(b > a)
     107       16030 :       IF (x < a .OR. x > b) THEN
     108             :          ! outside switching intervall
     109       14990 :          IF (order > 0) THEN
     110             :             ! derivatives are 0
     111             :             fx = 0.0_dp
     112             :          ELSE
     113        7495 :             IF (x < a) THEN
     114             :                ! x < a => f(x) = 1
     115             :                fx = 1.0_dp
     116             :             ELSE
     117             :                ! x > b => f(x) = 0
     118        7288 :                fx = 0.0_dp
     119             :             END IF
     120             :          END IF
     121             :       ELSE
     122             :          ! renormalized coordinate
     123        1040 :          u = (x - a)/(b - a)
     124        1560 :          SELECT CASE (order)
     125             :          CASE (0)
     126         520 :             u2 = u*u
     127         520 :             u3 = u2*u
     128         520 :             fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
     129             :          CASE (1)
     130         520 :             u2 = u*u
     131         520 :             fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
     132         520 :             fx = fx/(b - a)
     133             :          CASE (2)
     134           0 :             u2 = u*u
     135           0 :             fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
     136           0 :             fx = fx/(b - a)**2
     137             :          CASE DEFAULT
     138        1040 :             CPABORT('order not defined')
     139             :          END SELECT
     140             :       END IF
     141             : 
     142       16030 :    END FUNCTION pswitch
     143             : 
     144             : ! **************************************************************************************************
     145             : !> \brief determines if a value is not normal (e.g. for Inf and Nan)
     146             : !>        based on IO to work also under optimization.
     147             : !> \param a input value
     148             : !> \return TRUE for NaN and Inf
     149             : ! **************************************************************************************************
     150      344444 :    LOGICAL FUNCTION abnormal_value(a)
     151             :       REAL(KIND=dp)                                      :: a
     152             : 
     153             :       CHARACTER(LEN=32)                                  :: buffer
     154             : 
     155      344444 :       abnormal_value = .FALSE.
     156             :       ! the function should work when compiled with -ffast-math and similar
     157             :       ! unfortunately, that option asserts that all numbers are normals,
     158             :       ! which the compiler uses to optimize the function to .FALSE. if based on the IEEE module
     159             :       ! therefore, pass this to the Fortran runtime/printf, if things are NaN or Inf, error out.
     160      344444 :       WRITE (buffer, *) a
     161      344444 :       IF (INDEX(buffer, "N") .NE. 0 .OR. INDEX(buffer, "n") .NE. 0) abnormal_value = .TRUE.
     162             : 
     163      344444 :    END FUNCTION abnormal_value
     164             : 
     165             : ! **************************************************************************************************
     166             : !> \brief  Calculation of the angle between the vectors a and b.
     167             : !>         The angle is returned in radians.
     168             : !> \param a ...
     169             : !> \param b ...
     170             : !> \return ...
     171             : !> \date    14.10.1998
     172             : !> \author  MK
     173             : !> \version 1.0
     174             : ! **************************************************************************************************
     175      606672 :    PURE FUNCTION angle(a, b) RESULT(angle_ab)
     176             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: a, b
     177             :       REAL(KIND=dp)                                      :: angle_ab
     178             : 
     179             :       REAL(KIND=dp)                                      :: length_of_a, length_of_b
     180      606672 :       REAL(KIND=dp), DIMENSION(SIZE(a, 1))               :: a_norm, b_norm
     181             : 
     182     2426688 :       length_of_a = SQRT(DOT_PRODUCT(a, a))
     183     2426688 :       length_of_b = SQRT(DOT_PRODUCT(b, b))
     184             : 
     185      606672 :       IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
     186     2426688 :          a_norm(:) = a(:)/length_of_a
     187     2426688 :          b_norm(:) = b(:)/length_of_b
     188     2426688 :          angle_ab = ACOS(MIN(MAX(DOT_PRODUCT(a_norm, b_norm), -1.0_dp), 1.0_dp))
     189             :       ELSE
     190             :          angle_ab = 0.0_dp
     191             :       END IF
     192             : 
     193      606672 :    END FUNCTION angle
     194             : 
     195             : ! **************************************************************************************************
     196             : !> \brief   The binomial coefficient n over k for 0 <= k <= n is calculated,
     197             : !>            otherwise zero is returned.
     198             : !> \param n ...
     199             : !> \param k ...
     200             : !> \return ...
     201             : !> \date    08.03.1999
     202             : !> \author  MK
     203             : !> \version 1.0
     204             : ! **************************************************************************************************
     205   344248841 :    ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
     206             :       INTEGER, INTENT(IN)                                :: n, k
     207             :       REAL(KIND=dp)                                      :: n_over_k
     208             : 
     209   344248841 :       IF ((k >= 0) .AND. (k <= n)) THEN
     210   341590495 :          n_over_k = fac(n)/(fac(n - k)*fac(k))
     211             :       ELSE
     212             :          n_over_k = 0.0_dp
     213             :       END IF
     214             : 
     215   344248841 :    END FUNCTION binomial
     216             : 
     217             : ! **************************************************************************************************
     218             : !> \brief   The generalized binomial coefficient z over k for 0 <= k <= n is calculated.
     219             : !>            (z)   z*(z-1)*...*(z-k+2)*(z-k+1)
     220             : !>            ( ) = ---------------------------
     221             : !>            (k)                 k!
     222             : !> \param z ...
     223             : !> \param k ...
     224             : !> \return ...
     225             : !> \date    11.11.2019
     226             : !> \author  FS
     227             : !> \version 1.0
     228             : ! **************************************************************************************************
     229      171640 :    ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
     230             :       REAL(KIND=dp), INTENT(IN)                          :: z
     231             :       INTEGER, INTENT(IN)                                :: k
     232             :       REAL(KIND=dp)                                      :: z_over_k
     233             : 
     234             :       INTEGER                                            :: i
     235             : 
     236      171640 :       IF (k >= 0) THEN
     237             :          z_over_k = 1.0_dp
     238      257460 :          DO i = 1, k
     239      257460 :             z_over_k = z_over_k*(z - i + 1)/REAL(i, dp)
     240             :          END DO
     241             :       ELSE
     242             :          z_over_k = 0.0_dp
     243             :       END IF
     244             : 
     245      171640 :    END FUNCTION binomial_gen
     246             : 
     247             : ! **************************************************************************************************
     248             : !> \brief Calculates the multinomial coefficients
     249             : !> \param n ...
     250             : !> \param k ...
     251             : !> \return ...
     252             : !> \author Ole Schuett
     253             : ! **************************************************************************************************
     254        8982 :    PURE FUNCTION multinomial(n, k) RESULT(res)
     255             :       INTEGER, INTENT(IN)                                :: n
     256             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: k
     257             :       REAL(KIND=dp)                                      :: res
     258             : 
     259             :       INTEGER                                            :: i
     260             :       REAL(KIND=dp)                                      :: denom
     261             : 
     262       71856 :       IF (ALL(k >= 0) .AND. SUM(k) == n) THEN
     263        5280 :          denom = 1.0_dp
     264       21120 :          DO i = 1, SIZE(k)
     265       21120 :             denom = denom*fac(k(i))
     266             :          END DO
     267        5280 :          res = fac(n)/denom
     268             :       ELSE
     269             :          res = 0.0_dp
     270             :       END IF
     271             : 
     272        8982 :    END FUNCTION multinomial
     273             : 
     274             : ! **************************************************************************************************
     275             : !> \brief   The rotation matrix rotmat which rotates a vector about a
     276             : !>          rotation axis defined by the vector a is build up.
     277             : !>          The rotation angle is phi (radians).
     278             : !> \param phi ...
     279             : !> \param a ...
     280             : !> \param rotmat ...
     281             : !> \date    16.10.1998
     282             : !> \author  MK
     283             : !> \version 1.0
     284             : ! **************************************************************************************************
     285       15759 :    PURE SUBROUTINE build_rotmat(phi, a, rotmat)
     286             :       REAL(KIND=dp), INTENT(IN)                          :: phi
     287             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
     288             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: rotmat
     289             : 
     290             :       REAL(KIND=dp)                                      :: cosp, cost, length_of_a, sinp
     291             :       REAL(KIND=dp), DIMENSION(3)                        :: d
     292             : 
     293       15759 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
     294             :       ! Check the length of the vector a
     295       15759 :       IF (length_of_a > eps_geo) THEN
     296             : 
     297       63036 :          d(:) = a(:)/length_of_a
     298             : 
     299       15759 :          cosp = COS(phi)
     300       15759 :          sinp = SIN(phi)
     301       15759 :          cost = 1.0_dp - cosp
     302             : 
     303       15759 :          rotmat(1, 1) = d(1)*d(1)*cost + cosp
     304       15759 :          rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
     305       15759 :          rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
     306       15759 :          rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
     307       15759 :          rotmat(2, 2) = d(2)*d(2)*cost + cosp
     308       15759 :          rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
     309       15759 :          rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
     310       15759 :          rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
     311       15759 :          rotmat(3, 3) = d(3)*d(3)*cost + cosp
     312             :       ELSE
     313           0 :          CALL unit_matrix(rotmat)
     314             :       END IF
     315             : 
     316       15759 :    END SUBROUTINE build_rotmat
     317             : 
     318             : ! **************************************************************************************************
     319             : !> \brief   Returns the determinante of the 3x3 matrix a.
     320             : !> \param a ...
     321             : !> \return ...
     322             : !> \date    13.03.2004
     323             : !> \author  MK
     324             : !> \version 1.0
     325             : ! **************************************************************************************************
     326    16465682 :    PURE FUNCTION det_3x3_1(a) RESULT(det_a)
     327             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: a
     328             :       REAL(KIND=dp)                                      :: det_a
     329             : 
     330             :       det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
     331             :               a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
     332    16465682 :               a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
     333             : 
     334    16465682 :    END FUNCTION det_3x3_1
     335             : 
     336             : ! **************************************************************************************************
     337             : !> \brief   Returns the determinante of the 3x3 matrix a given by its columns.
     338             : !> \param a1 ...
     339             : !> \param a2 ...
     340             : !> \param a3 ...
     341             : !> \return ...
     342             : !> \date    13.03.2004
     343             : !> \author  MK
     344             : !> \version 1.0
     345             : ! **************************************************************************************************
     346        6201 :    PURE FUNCTION det_3x3_2(a1, a2, a3) RESULT(det_a)
     347             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a1, a2, a3
     348             :       REAL(KIND=dp)                                      :: det_a
     349             : 
     350             :       det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
     351             :               a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
     352        6201 :               a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
     353             : 
     354        6201 :    END FUNCTION det_3x3_2
     355             : 
     356             : ! **************************************************************************************************
     357             : !> \brief Diagonalize the symmetric n by n matrix a using the LAPACK
     358             : !>        library. Only the upper triangle of matrix a is used.
     359             : !>        Externals (LAPACK 3.0)
     360             : !> \param a ...
     361             : !> \param eigval ...
     362             : !> \param dac ...
     363             : !> \date    29.03.1999
     364             : !> \par Variables
     365             : !>      - a       : Symmetric matrix to be diagonalized (input; upper triangle) ->
     366             : !>      -           eigenvectors of the matrix a (output).
     367             : !>      - dac     : If true, then the divide-and-conquer algorithm is applied.
     368             : !>      - eigval  : Eigenvalues of the matrix a (output).
     369             : !> \author  MK
     370             : !> \version 1.0
     371             : ! **************************************************************************************************
     372       74776 :    SUBROUTINE diamat_all(a, eigval, dac)
     373             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
     374             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigval
     375             :       LOGICAL, INTENT(IN), OPTIONAL                      :: dac
     376             : 
     377             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'diamat_all'
     378             : 
     379             :       INTEGER                                            :: handle, info, liwork, lwork, n, nb
     380       74776 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
     381             :       INTEGER, EXTERNAL                                  :: ilaenv
     382             :       LOGICAL                                            :: divide_and_conquer
     383       74776 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work
     384             : 
     385             :       EXTERNAL dsyev, dsyevd
     386             : 
     387       74776 :       CALL timeset(routineN, handle)
     388             : 
     389             :       ! Get the size of the matrix a
     390       74776 :       n = SIZE(a, 1)
     391             : 
     392             :       ! Check the size of matrix a
     393       74776 :       IF (SIZE(a, 2) /= n) THEN
     394           0 :          CPABORT("Check the size of matrix a (parameter #1)")
     395             :       END IF
     396             : 
     397             :       ! Check the size of vector eigval
     398       74776 :       IF (SIZE(eigval) /= n) THEN
     399           0 :          CPABORT("The dimension of vector eigval is too small")
     400             :       END IF
     401             : 
     402             :       ! Check, if the divide-and-conquer algorithm is requested
     403             : 
     404       74776 :       IF (PRESENT(dac)) THEN
     405         205 :          divide_and_conquer = dac
     406             :       ELSE
     407             :          divide_and_conquer = .FALSE.
     408             :       END IF
     409             : 
     410             :       ! Get the optimal work storage size
     411             : 
     412         205 :       IF (divide_and_conquer) THEN
     413         205 :          lwork = 2*n**2 + 6*n + 1
     414         205 :          liwork = 5*n + 3
     415             :       ELSE
     416       74571 :          nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
     417       74571 :          lwork = (nb + 2)*n
     418             :       END IF
     419             : 
     420             :       ! Allocate work storage
     421             : 
     422      224328 :       ALLOCATE (work(lwork))
     423       74776 :       IF (divide_and_conquer) THEN
     424         615 :          ALLOCATE (iwork(liwork))
     425             :       END IF
     426             : 
     427             :       ! Diagonalize the matrix a
     428             : 
     429       74776 :       info = 0
     430             :       IF (divide_and_conquer) THEN
     431         205 :          CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
     432             :       ELSE
     433       74767 :          CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
     434             :       END IF
     435             : 
     436       74776 :       IF (info /= 0) THEN
     437           0 :          IF (divide_and_conquer) THEN
     438           0 :             CPABORT("The matrix diagonalization with dsyevd failed")
     439             :          ELSE
     440           0 :             CPABORT("The matrix diagonalization with dsyev failed")
     441             :          END IF
     442             :       END IF
     443             : 
     444             :       ! Release work storage
     445       74776 :       DEALLOCATE (work)
     446             : 
     447       74776 :       IF (divide_and_conquer) THEN
     448         205 :          DEALLOCATE (iwork)
     449             :       END IF
     450             : 
     451       74776 :       CALL timestop(handle)
     452             : 
     453      149552 :    END SUBROUTINE diamat_all
     454             : 
     455             : ! **************************************************************************************************
     456             : !> \brief   Returns the dihedral angle, i.e. the angle between the planes
     457             : !>          defined by the vectors (-ab,bc) and (cd,-bc).
     458             : !>          The dihedral angle is returned in radians.
     459             : !> \param ab ...
     460             : !> \param bc ...
     461             : !> \param cd ...
     462             : !> \return ...
     463             : !> \date    13.03.2004
     464             : !> \author  MK
     465             : !> \version 1.0
     466             : ! **************************************************************************************************
     467         116 :    PURE FUNCTION dihedral_angle(ab, bc, cd) RESULT(dihedral_angle_abcd)
     468             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: ab, bc, cd
     469             :       REAL(KIND=dp)                                      :: dihedral_angle_abcd
     470             : 
     471             :       REAL(KIND=dp)                                      :: det_abcd
     472             :       REAL(KIND=dp), DIMENSION(3)                        :: abc, bcd
     473             : 
     474         464 :       abc = vector_product(bc, -ab)
     475         464 :       bcd = vector_product(cd, -bc)
     476             :       ! Calculate the normal vectors of the planes
     477             :       ! defined by the points a,b,c and b,c,d
     478             : 
     479         464 :       det_abcd = det_3x3(abc, bcd, -bc)
     480         116 :       dihedral_angle_abcd = SIGN(1.0_dp, det_abcd)*angle(abc, bcd)
     481             : 
     482         116 :    END FUNCTION dihedral_angle
     483             : 
     484             : ! **************************************************************************************************
     485             : !> \brief   Return the diagonal elements of matrix a as a vector.
     486             : !> \param a ...
     487             : !> \return ...
     488             : !> \date    20.11.1998
     489             : !> \author  MK
     490             : !> \version 1.0
     491             : ! **************************************************************************************************
     492       35680 :    PURE FUNCTION get_diag(a) RESULT(a_diag)
     493             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: a
     494             :       REAL(KIND=dp), &
     495             :          DIMENSION(MIN(SIZE(a, 1), SIZE(a, 2)))          :: a_diag
     496             : 
     497             :       INTEGER                                            :: i, n
     498             : 
     499       35680 :       n = MIN(SIZE(a, 1), SIZE(a, 2))
     500             : 
     501     3612682 :       DO i = 1, n
     502     3612682 :          a_diag(i) = a(i, i)
     503             :       END DO
     504             : 
     505       35680 :    END FUNCTION get_diag
     506             : 
     507             : ! **************************************************************************************************
     508             : !> \brief   Returns the inverse of the 3 x 3 matrix a.
     509             : !> \param a ...
     510             : !> \return ...
     511             : !> \date    13.03.2004
     512             : !> \author  MK
     513             : !> \version 1.0
     514             : ! **************************************************************************************************
     515    15705655 :    PURE FUNCTION inv_3x3(a) RESULT(a_inv)
     516             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: a
     517             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: a_inv
     518             : 
     519             :       REAL(KIND=dp)                                      :: det_a
     520             : 
     521    15705655 :       det_a = 1.0_dp/det_3x3(a)
     522             : 
     523    15705655 :       a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
     524    15705655 :       a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
     525    15705655 :       a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
     526             : 
     527    15705655 :       a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
     528    15705655 :       a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
     529    15705655 :       a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
     530             : 
     531    15705655 :       a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
     532    15705655 :       a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
     533    15705655 :       a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
     534             : 
     535    15705655 :    END FUNCTION inv_3x3
     536             : 
     537             : ! **************************************************************************************************
     538             : !> \brief returns inverse of matrix using the lapack routines DGETRF and DGETRI
     539             : !> \param a ...
     540             : !> \param info ...
     541             : ! **************************************************************************************************
     542        6062 :    SUBROUTINE invmat(a, info)
     543             :       REAL(KIND=dp), INTENT(INOUT)                       :: a(:, :)
     544             :       INTEGER, INTENT(OUT)                               :: info
     545             : 
     546             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'invmat'
     547             : 
     548             :       INTEGER                                            :: handle, lwork, n
     549        6062 :       INTEGER, ALLOCATABLE                               :: ipiv(:)
     550        6062 :       REAL(KIND=dp), ALLOCATABLE                         :: work(:)
     551             : 
     552        6062 :       CALL timeset(routineN, handle)
     553             : 
     554        6062 :       n = SIZE(a, 1)
     555        6062 :       lwork = 20*n
     556       18186 :       ALLOCATE (ipiv(n))
     557       18186 :       ALLOCATE (work(lwork))
     558       57762 :       ipiv = 0
     559     1040062 :       work = 0._dp
     560        6062 :       info = 0
     561      438574 :       CALL dgetrf(n, n, a, n, ipiv, info)
     562        6062 :       IF (info == 0) THEN
     563      438544 :          CALL dgetri(n, a, n, ipiv, work, lwork, info)
     564             :       END IF
     565        6062 :       DEALLOCATE (ipiv, work)
     566             : 
     567        6062 :       CALL timestop(handle)
     568             : 
     569        6062 :    END SUBROUTINE invmat
     570             : 
     571             : ! **************************************************************************************************
     572             : !> \brief returns inverse of real symmetric, positive definite matrix
     573             : !> \param a matrix
     574             : !> \param potrf if cholesky decomposition of a was already done using dpotrf.
     575             : !>        If not given, cholesky decomposition of a will be done before inversion.
     576             : !> \param uplo indicating if the upper or lower triangle of a is stored.
     577             : !> \author Dorothea Golze [02.2015]
     578             : ! **************************************************************************************************
     579        2401 :    SUBROUTINE invmat_symm(a, potrf, uplo)
     580             :       REAL(KIND=dp), INTENT(INOUT)                       :: a(:, :)
     581             :       LOGICAL, INTENT(IN), OPTIONAL                      :: potrf
     582             :       CHARACTER(LEN=1), INTENT(IN), OPTIONAL             :: uplo
     583             : 
     584             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'invmat_symm'
     585             : 
     586             :       CHARACTER(LEN=1)                                   :: myuplo
     587             :       INTEGER                                            :: handle, info, n
     588             :       LOGICAL                                            :: do_potrf
     589             : 
     590        2401 :       CALL timeset(routineN, handle)
     591             : 
     592        2401 :       do_potrf = .TRUE.
     593        2401 :       IF (PRESENT(potrf)) do_potrf = potrf
     594             : 
     595        2401 :       myuplo = 'U'
     596        2401 :       IF (PRESENT(uplo)) myuplo = uplo
     597             : 
     598        2401 :       n = SIZE(a, 1)
     599        2401 :       info = 0
     600             : 
     601             :       ! do cholesky decomposition
     602        2401 :       IF (do_potrf) THEN
     603        3384 :          CALL dpotrf(myuplo, n, a, n, info)
     604        2184 :          IF (info /= 0) CPABORT("DPOTRF failed")
     605             :       END IF
     606             : 
     607             :       ! do inversion using the cholesky decomposition
     608        3601 :       CALL dpotri(myuplo, n, a, n, info)
     609        2401 :       IF (info /= 0) CPABORT("Matrix inversion failed")
     610             : 
     611             :       ! complete the matrix
     612        2401 :       IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
     613        2401 :          CALL symmetrize_matrix(a, "upper_to_lower")
     614             :       ELSE
     615           0 :          CALL symmetrize_matrix(a, "lower_to_upper")
     616             :       END IF
     617             : 
     618        2401 :       CALL timestop(handle)
     619             : 
     620        2401 :    END SUBROUTINE invmat_symm
     621             : 
     622             : ! **************************************************************************************************
     623             : !> \brief  Compute the inverse of the n by n real matrix a using the LAPACK
     624             : !>         library
     625             : !> \param a ...
     626             : !> \param a_inverse ...
     627             : !> \param eval_error ...
     628             : !> \param option ...
     629             : !> \param improve ...
     630             : !> \date   23.03.1999
     631             : !> \par Variables
     632             : !>       - a        : Real matrix to be inverted (input).
     633             : !>       - a_inverse: Inverse of the matrix a (output).
     634             : !>       - a_lu     : LU factorization of matrix a.
     635             : !>       - a_norm   : Norm of matrix a.
     636             : !>       - error    : Estimated error of the inversion.
     637             : !>       - r_cond   : Reciprocal condition number of the matrix a.
     638             : !>       - trans    : "N" => invert a
     639             : !>       -            "T" => invert transpose(a)
     640             : !> \author MK
     641             : !> \version 1.0
     642             : !> \note NB add improve argument, used to disable call to dgerfs
     643             : ! **************************************************************************************************
     644        2066 :    SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
     645             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: a
     646             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: a_inverse
     647             :       REAL(KIND=dp), INTENT(OUT)                         :: eval_error
     648             :       CHARACTER(LEN=1), INTENT(IN), OPTIONAL             :: option
     649             :       LOGICAL, INTENT(IN), OPTIONAL                      :: improve
     650             : 
     651             :       CHARACTER(LEN=1)                                   :: norm, trans
     652             :       CHARACTER(LEN=default_string_length)               :: message
     653             :       INTEGER                                            :: info, iter, n
     654        2066 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ipiv, iwork
     655             :       LOGICAL                                            :: do_improve
     656             :       REAL(KIND=dp)                                      :: a_norm, old_eval_error, r_cond
     657        2066 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: berr, ferr, work
     658        2066 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: a_lu, b
     659             :       REAL(KIND=dp), EXTERNAL                            :: dlange
     660             : 
     661             :       EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
     662             : 
     663             :       ! Check for optional parameter
     664        2066 :       IF (PRESENT(option)) THEN
     665         468 :          trans = option
     666             :       ELSE
     667        1598 :          trans = "N"
     668             :       END IF
     669             : 
     670        2066 :       IF (PRESENT(improve)) THEN
     671         534 :          do_improve = improve
     672             :       ELSE
     673             :          do_improve = .TRUE.
     674             :       END IF
     675             : 
     676             :       ! Get the dimension of matrix a
     677        2066 :       n = SIZE(a, 1)
     678             : 
     679             :       ! Check array dimensions
     680        2066 :       IF (n == 0) THEN
     681           0 :          CPABORT("Matrix to be inverted of zero size")
     682             :       END IF
     683             : 
     684        2066 :       IF (n /= SIZE(a, 2)) THEN
     685           0 :          CPABORT("Check the array bounds of parameter #1")
     686             :       END IF
     687             : 
     688        2066 :       IF ((n /= SIZE(a_inverse, 1)) .OR. &
     689             :           (n /= SIZE(a_inverse, 2))) THEN
     690           0 :          CPABORT("Check the array bounds of parameter #2")
     691             :       END IF
     692             : 
     693             :       ! Allocate work storage
     694        8264 :       ALLOCATE (a_lu(n, n))
     695        6198 :       ALLOCATE (b(n, n))
     696        6198 :       ALLOCATE (berr(n))
     697        4132 :       ALLOCATE (ferr(n))
     698        6198 :       ALLOCATE (ipiv(n))
     699        4132 :       ALLOCATE (iwork(n))
     700        6198 :       ALLOCATE (work(4*n))
     701             : 
     702     2451814 :       a_lu(1:n, 1:n) = a(1:n, 1:n)
     703             : 
     704             :       ! Compute the LU factorization of the matrix a
     705        2066 :       CALL dgetrf(n, n, a_lu, n, ipiv, info)
     706             : 
     707        2066 :       IF (info /= 0) THEN
     708           0 :          CPABORT("The LU factorization in dgetrf failed")
     709             :       END IF
     710             : 
     711             :       ! Compute the norm of the matrix a
     712             : 
     713        2066 :       IF (trans == "N") THEN
     714        1844 :          norm = '1'
     715             :       ELSE
     716         222 :          norm = 'I'
     717             :       END IF
     718             : 
     719        2066 :       a_norm = dlange(norm, n, n, a, n, work)
     720             : 
     721             :       ! Compute the reciprocal of the condition number of a
     722             : 
     723        2066 :       CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
     724             : 
     725        2066 :       IF (info /= 0) THEN
     726           0 :          CPABORT("The computation of the condition number in dgecon failed")
     727             :       END IF
     728             : 
     729        2066 :       IF (r_cond < EPSILON(0.0_dp)) THEN
     730           0 :          WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
     731             :          CALL cp_abort(__LOCATION__, &
     732             :                        "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
     733           0 :                        "working precision)")
     734             :       END IF
     735             : 
     736             :       ! Solve a system of linear equations using the LU factorization computed by dgetrf
     737             : 
     738        2066 :       CALL unit_matrix(a_inverse)
     739             : 
     740        2066 :       CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
     741             : 
     742        2066 :       IF (info /= 0) THEN
     743           0 :          CPABORT("Solving the system of linear equations in dgetrs failed")
     744             :       END IF
     745             : 
     746             :       ! Improve the computed solution iteratively
     747        2066 :       CALL unit_matrix(b) ! Initialize right-hand sides
     748             : 
     749        2066 :       eval_error = 0.0_dp
     750             : 
     751        2066 :       IF (do_improve) THEN
     752        5932 :          DO iter = 1, 10
     753             : 
     754             :             CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
     755        5706 :                         work, iwork, info)
     756             : 
     757        5706 :             IF (info /= 0) THEN
     758           0 :                CPABORT("Improving the computed solution in dgerfs failed")
     759             :             END IF
     760             : 
     761        5706 :             old_eval_error = eval_error
     762      167334 :             eval_error = MAXVAL(ferr)
     763             : 
     764        5932 :             IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
     765             : 
     766             :          END DO
     767             :       END IF
     768             : 
     769             :       ! Release work storage
     770        2066 :       DEALLOCATE (work)
     771        2066 :       DEALLOCATE (iwork)
     772        2066 :       DEALLOCATE (ipiv)
     773        2066 :       DEALLOCATE (ferr)
     774        2066 :       DEALLOCATE (berr)
     775        2066 :       DEALLOCATE (b)
     776        2066 :       DEALLOCATE (a_lu)
     777             : 
     778        2066 :    END SUBROUTINE invert_matrix_d
     779             : 
     780             : ! **************************************************************************************************
     781             : !> \brief  Compute the inverse of the n by n complex matrix a using the LAPACK
     782             : !>         library
     783             : !> \param a ...
     784             : !> \param a_inverse ...
     785             : !> \param eval_error ...
     786             : !> \param option ...
     787             : !> \date   08.06.2009
     788             : !> \par Variables
     789             : !>       - a        : Complex matrix to be inverted (input).
     790             : !>       - a_inverse: Inverse of the matrix a (output).
     791             : !>       - a_lu     : LU factorization of matrix a.
     792             : !>       - a_norm   : Norm of matrix a.
     793             : !>       - error    : Estimated error of the inversion.
     794             : !>       - r_cond   : Reciprocal condition number of the matrix a.
     795             : !>       - trans    : "N" => invert a
     796             : !>       -            "T" => invert transpose(a)
     797             : !> \author MK
     798             : !> \version 1.0
     799             : ! **************************************************************************************************
     800           0 :    SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
     801             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN)      :: a
     802             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT)     :: a_inverse
     803             :       REAL(KIND=dp), INTENT(OUT)                         :: eval_error
     804             :       CHARACTER(LEN=1), INTENT(IN), OPTIONAL             :: option
     805             : 
     806             :       CHARACTER(LEN=1)                                   :: norm, trans
     807             :       CHARACTER(LEN=default_string_length)               :: message
     808           0 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: work
     809           0 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :)     :: a_lu, b
     810             :       INTEGER                                            :: info, iter, n
     811           0 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ipiv
     812             :       REAL(KIND=dp)                                      :: a_norm, old_eval_error, r_cond
     813           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: berr, ferr, rwork
     814             :       REAL(KIND=dp), EXTERNAL                            :: zlange
     815             : 
     816             :       EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
     817             : 
     818             :       ! Check for optional parameter
     819           0 :       IF (PRESENT(option)) THEN
     820           0 :          trans = option
     821             :       ELSE
     822           0 :          trans = "N"
     823             :       END IF
     824             : 
     825             :       ! Get the dimension of matrix a
     826           0 :       n = SIZE(a, 1)
     827             : 
     828             :       ! Check array dimensions
     829           0 :       IF (n == 0) THEN
     830           0 :          CPABORT("Matrix to be inverted of zero size")
     831             :       END IF
     832             : 
     833           0 :       IF (n /= SIZE(a, 2)) THEN
     834           0 :          CPABORT("Check the array bounds of parameter #1")
     835             :       END IF
     836             : 
     837           0 :       IF ((n /= SIZE(a_inverse, 1)) .OR. &
     838             :           (n /= SIZE(a_inverse, 2))) THEN
     839           0 :          CPABORT("Check the array bounds of parameter #2")
     840             :       END IF
     841             : 
     842             :       ! Allocate work storage
     843           0 :       ALLOCATE (a_lu(n, n))
     844           0 :       ALLOCATE (b(n, n))
     845           0 :       ALLOCATE (berr(n))
     846           0 :       ALLOCATE (ferr(n))
     847           0 :       ALLOCATE (ipiv(n))
     848           0 :       ALLOCATE (rwork(2*n))
     849           0 :       ALLOCATE (work(2*n))
     850             : 
     851           0 :       a_lu(1:n, 1:n) = a(1:n, 1:n)
     852             : 
     853             :       ! Compute the LU factorization of the matrix a
     854           0 :       CALL zgetrf(n, n, a_lu, n, ipiv, info)
     855             : 
     856           0 :       IF (info /= 0) THEN
     857           0 :          CPABORT("The LU factorization in dgetrf failed")
     858             :       END IF
     859             : 
     860             :       ! Compute the norm of the matrix a
     861             : 
     862           0 :       IF (trans == "N") THEN
     863           0 :          norm = '1'
     864             :       ELSE
     865           0 :          norm = 'I'
     866             :       END IF
     867             : 
     868           0 :       a_norm = zlange(norm, n, n, a, n, work)
     869             : 
     870             :       ! Compute the reciprocal of the condition number of a
     871             : 
     872           0 :       CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
     873             : 
     874           0 :       IF (info /= 0) THEN
     875           0 :          CPABORT("The computation of the condition number in dgecon failed")
     876             :       END IF
     877             : 
     878           0 :       IF (r_cond < EPSILON(0.0_dp)) THEN
     879           0 :          WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
     880             :          CALL cp_abort(__LOCATION__, &
     881             :                        "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
     882           0 :                        "working precision)")
     883             :       END IF
     884             : 
     885             :       ! Solve a system of linear equations using the LU factorization computed by dgetrf
     886             : 
     887           0 :       CALL unit_matrix(a_inverse)
     888             : 
     889           0 :       CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
     890             : 
     891           0 :       IF (info /= 0) THEN
     892           0 :          CPABORT("Solving the system of linear equations in dgetrs failed")
     893             :       END IF
     894             : 
     895             :       ! Improve the computed solution iteratively
     896           0 :       CALL unit_matrix(b) ! Initialize right-hand sides
     897             : 
     898           0 :       eval_error = 0.0_dp
     899             : 
     900           0 :       DO iter = 1, 10
     901             : 
     902             :          CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
     903           0 :                      work, rwork, info)
     904             : 
     905           0 :          IF (info /= 0) THEN
     906           0 :             CPABORT("Improving the computed solution in dgerfs failed")
     907             :          END IF
     908             : 
     909           0 :          old_eval_error = eval_error
     910           0 :          eval_error = MAXVAL(ferr)
     911             : 
     912           0 :          IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
     913             : 
     914             :       END DO
     915             : 
     916             :       ! Release work storage
     917           0 :       DEALLOCATE (work)
     918           0 :       DEALLOCATE (rwork)
     919           0 :       DEALLOCATE (ipiv)
     920           0 :       DEALLOCATE (ferr)
     921           0 :       DEALLOCATE (berr)
     922           0 :       DEALLOCATE (b)
     923           0 :       DEALLOCATE (a_lu)
     924             : 
     925           0 :    END SUBROUTINE invert_matrix_z
     926             : 
     927             : ! **************************************************************************************************
     928             : !> \brief returns the pseudoinverse of a real, square matrix using singular
     929             : !>         value decomposition
     930             : !> \param a matrix a
     931             : !> \param a_pinverse pseudoinverse of matrix a
     932             : !> \param rskip parameter for setting small singular values to zero
     933             : !> \param determinant determinant of matrix a (optional output)
     934             : !> \param sval array holding singular values of matrix a (optional output)
     935             : !> \author Dorothea Golze [02.2015]
     936             : ! **************************************************************************************************
     937         487 :    SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
     938             :       REAL(KIND=dp), DIMENSION(:, :)                     :: a, a_pinverse
     939             :       REAL(KIND=dp), INTENT(IN)                          :: rskip
     940             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: determinant
     941             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
     942             :          OPTIONAL, POINTER                               :: sval
     943             : 
     944             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_svd'
     945             : 
     946             :       INTEGER                                            :: handle, i, info, lwork, n
     947         487 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
     948         487 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: sig, work
     949         487 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: sig_plus, temp_mat, u, vt
     950             : 
     951         487 :       CALL timeset(routineN, handle)
     952             : 
     953         487 :       n = SIZE(a, 1)
     954        6818 :       ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
     955       48311 :       u(:, :) = 0.0_dp
     956       48311 :       vt(:, :) = 0.0_dp
     957        2677 :       sig(:) = 0.0_dp
     958       48311 :       sig_plus = 0.0_dp
     959         974 :       work = 0.0_dp
     960       18007 :       iwork = 0
     961         487 :       IF (PRESENT(determinant)) determinant = 1.0_dp
     962             : 
     963             :       ! work size query
     964         487 :       lwork = -1
     965             :       CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
     966         487 :                   lwork, iwork(1), info)
     967             : 
     968         487 :       IF (info /= 0) THEN
     969           0 :          CPABORT("ERROR in DGESDD: Could not retrieve work array sizes")
     970             :       END IF
     971         487 :       lwork = INT(work(1))
     972         487 :       DEALLOCATE (work)
     973        1461 :       ALLOCATE (work(lwork))
     974             : 
     975             :       ! do SVD
     976             :       CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
     977         487 :                   lwork, iwork(1), info)
     978             : 
     979         487 :       IF (info /= 0) THEN
     980           0 :          CPABORT("SVD failed")
     981             :       END IF
     982             : 
     983         487 :       IF (PRESENT(sval)) THEN
     984           0 :          CPASSERT(.NOT. ASSOCIATED(sval))
     985           0 :          ALLOCATE (sval(n))
     986           0 :          sval(:) = sig
     987             :       END IF
     988             : 
     989             :       ! set singular values that are too small to zero
     990        2677 :       DO i = 1, n
     991       50501 :          IF (sig(i) > rskip*MAXVAL(sig)) THEN
     992        2154 :             IF (PRESENT(determinant)) &
     993           0 :                determinant = determinant*sig(i)
     994        2154 :             sig_plus(i, i) = 1._dp/sig(i)
     995             :          ELSE
     996          36 :             sig_plus(i, i) = 0.0_dp
     997             :          END IF
     998             :       END DO
     999             : 
    1000             :       ! build pseudoinverse: V*sig_plus*UT
    1001         487 :       CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
    1002         487 :       CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
    1003             : 
    1004         487 :       DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
    1005             : 
    1006         487 :       CALL timestop(handle)
    1007             : 
    1008         487 :    END SUBROUTINE get_pseudo_inverse_svd
    1009             : 
    1010             : ! **************************************************************************************************
    1011             : !> \brief returns the pseudoinverse of a real, symmetric and positive definite
    1012             : !>        matrix using diagonalization.
    1013             : !> \param a matrix a
    1014             : !> \param a_pinverse pseudoinverse of matrix a
    1015             : !> \param rskip parameter for setting small eigenvalues to zero
    1016             : !> \author Dorothea Golze [02.2015]
    1017             : ! **************************************************************************************************
    1018        1161 :    SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
    1019             :       REAL(KIND=dp), DIMENSION(:, :)                     :: a, a_pinverse
    1020             :       REAL(KIND=dp), INTENT(IN)                          :: rskip
    1021             : 
    1022             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_diag'
    1023             : 
    1024             :       INTEGER                                            :: handle, i, info, lwork, n
    1025        1161 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eig, work
    1026        1161 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dinv, temp_mat
    1027             : 
    1028        1161 :       CALL timeset(routineN, handle)
    1029             : 
    1030        1161 :       info = 0
    1031        1161 :       n = SIZE(a, 1)
    1032        9288 :       ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
    1033      186519 :       dinv(:, :) = 0.0_dp
    1034       15075 :       eig(:) = 0.0_dp
    1035        2322 :       work(:) = 0.0_dp
    1036      186519 :       temp_mat = 0.0_dp
    1037             : 
    1038             :       ! work size query
    1039        1161 :       lwork = -1
    1040        1161 :       CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
    1041        1161 :       IF (info /= 0) THEN
    1042           0 :          CPABORT("ERROR in DSYEV: Could not retrieve work array sizes")
    1043             :       END IF
    1044        1161 :       lwork = INT(work(1))
    1045        1161 :       DEALLOCATE (work)
    1046        3483 :       ALLOCATE (work(lwork))
    1047      474237 :       work = 0.0_dp
    1048             : 
    1049             :       ! get eigenvalues and eigenvectors
    1050        1161 :       CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
    1051             : 
    1052        1161 :       IF (info /= 0) THEN
    1053           0 :          CPABORT("Matrix diagonalization failed")
    1054             :       END IF
    1055             : 
    1056             :       ! set eigenvalues that are too small to zero
    1057       15075 :       DO i = 1, n
    1058      200433 :          IF (eig(i) > rskip*MAXVAL(eig)) THEN
    1059       12454 :             dinv(i, i) = 1.0_dp/eig(i)
    1060             :          ELSE
    1061        1460 :             dinv(i, i) = 0._dp
    1062             :          END IF
    1063             :       END DO
    1064             : 
    1065             :       ! build pseudoinverse: U*dinv*UT
    1066        1161 :       CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
    1067        1161 :       CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
    1068             : 
    1069        1161 :       DEALLOCATE (eig, work, dinv, temp_mat)
    1070             : 
    1071        1161 :       CALL timestop(handle)
    1072             : 
    1073        1161 :    END SUBROUTINE get_pseudo_inverse_diag
    1074             : 
    1075             : ! **************************************************************************************************
    1076             : !> \brief  Reflection of the vector a through a mirror plane defined by the
    1077             : !>         normal vector b. The reflected vector a is stored in a_mirror.
    1078             : !> \param a ...
    1079             : !> \param b ...
    1080             : !> \return ...
    1081             : !> \date    16.10.1998
    1082             : !> \author  MK
    1083             : !> \version 1.0
    1084             : ! **************************************************************************************************
    1085       19487 :    PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
    1086             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a, b
    1087             :       REAL(KIND=dp), DIMENSION(3)                        :: a_mirror
    1088             : 
    1089             :       REAL(KIND=dp)                                      :: length_of_b, scapro
    1090             :       REAL(KIND=dp), DIMENSION(3)                        :: d
    1091             : 
    1092       19487 :       length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
    1093             : 
    1094       19487 :       IF (length_of_b > eps_geo) THEN
    1095             : 
    1096       77948 :          d(:) = b(:)/length_of_b
    1097             : 
    1098             :          ! Calculate the mirror image a_mirror of the vector a
    1099       19487 :          scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
    1100             : 
    1101       77948 :          a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
    1102             : 
    1103             :       ELSE
    1104             : 
    1105           0 :          a_mirror(:) = 0.0_dp
    1106             : 
    1107             :       END IF
    1108             : 
    1109       19487 :    END FUNCTION reflect_vector
    1110             : 
    1111             : ! **************************************************************************************************
    1112             : !> \brief   Rotation of the vector a about an rotation axis defined by the
    1113             : !>          vector b. The rotation angle is phi (radians). The rotated vector
    1114             : !>          a is stored in a_rot.
    1115             : !> \param a ...
    1116             : !> \param phi ...
    1117             : !> \param b ...
    1118             : !> \return ...
    1119             : !> \date    16.10.1998
    1120             : !> \author  MK
    1121             : !> \version 1.0
    1122             : ! **************************************************************************************************
    1123        4098 :    PURE FUNCTION rotate_vector(a, phi, b) RESULT(a_rot)
    1124             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1125             :       REAL(KIND=dp), INTENT(IN)                          :: phi
    1126             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: b
    1127             :       REAL(KIND=dp), DIMENSION(3)                        :: a_rot
    1128             : 
    1129             :       REAL(KIND=dp)                                      :: length_of_b
    1130             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rotmat
    1131             : 
    1132        4098 :       length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
    1133        4098 :       IF (length_of_b > eps_geo) THEN
    1134             : 
    1135             :          ! Build up the rotation matrix rotmat
    1136        3964 :          CALL build_rotmat(phi, b, rotmat)
    1137             : 
    1138             :          ! Rotate the vector a by phi about the axis defined by vector b
    1139       63424 :          a_rot(:) = MATMUL(rotmat, a)
    1140             : 
    1141             :       ELSE
    1142             : 
    1143         536 :          a_rot(:) = 0.0_dp
    1144             : 
    1145             :       END IF
    1146             : 
    1147        4098 :    END FUNCTION rotate_vector
    1148             : 
    1149             : ! **************************************************************************************************
    1150             : !> \brief   Set the diagonal elements of matrix a to b.
    1151             : !> \param a ...
    1152             : !> \param b ...
    1153             : !> \date    20.11.1998
    1154             : !> \author  MK
    1155             : !> \version 1.0
    1156             : ! **************************************************************************************************
    1157       36239 :    PURE SUBROUTINE set_diag_scalar_d(a, b)
    1158             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
    1159             :       REAL(KIND=dp), INTENT(IN)                          :: b
    1160             : 
    1161             :       INTEGER                                            :: i, n
    1162             : 
    1163       36239 :       n = MIN(SIZE(a, 1), SIZE(a, 2))
    1164      200835 :       DO i = 1, n
    1165      200835 :          a(i, i) = b
    1166             :       END DO
    1167             : 
    1168       36239 :    END SUBROUTINE set_diag_scalar_d
    1169             : 
    1170             : ! **************************************************************************************************
    1171             : !> \brief ...
    1172             : !> \param a ...
    1173             : !> \param b ...
    1174             : ! **************************************************************************************************
    1175           0 :    PURE SUBROUTINE set_diag_scalar_z(a, b)
    1176             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT)   :: a
    1177             :       COMPLEX(KIND=dp), INTENT(IN)                       :: b
    1178             : 
    1179             :       INTEGER                                            :: i, n
    1180             : 
    1181           0 :       n = MIN(SIZE(a, 1), SIZE(a, 2))
    1182           0 :       DO i = 1, n
    1183           0 :          a(i, i) = b
    1184             :       END DO
    1185             : 
    1186           0 :    END SUBROUTINE set_diag_scalar_z
    1187             : 
    1188             : ! **************************************************************************************************
    1189             : !> \brief   Symmetrize the matrix a.
    1190             : !> \param a ...
    1191             : !> \param option ...
    1192             : !> \date    16.10.1998
    1193             : !> \author  MK
    1194             : !> \version 1.0
    1195             : ! **************************************************************************************************
    1196       16473 :    SUBROUTINE symmetrize_matrix(a, option)
    1197             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
    1198             :       CHARACTER(LEN=*), INTENT(IN)                       :: option
    1199             : 
    1200             :       INTEGER                                            :: i, n
    1201             : 
    1202       16473 :       n = MIN(SIZE(a, 1), SIZE(a, 2))
    1203             : 
    1204       16473 :       IF (option == "lower_to_upper") THEN
    1205           0 :          DO i = 1, n - 1
    1206           0 :             a(i, i + 1:n) = a(i + 1:n, i)
    1207             :          END DO
    1208       16473 :       ELSE IF (option == "upper_to_lower") THEN
    1209      403937 :          DO i = 1, n - 1
    1210    22586006 :             a(i + 1:n, i) = a(i, i + 1:n)
    1211             :          END DO
    1212          36 :       ELSE IF (option == "anti_lower_to_upper") THEN
    1213           0 :          DO i = 1, n - 1
    1214           0 :             a(i, i + 1:n) = -a(i + 1:n, i)
    1215             :          END DO
    1216          36 :       ELSE IF (option == "anti_upper_to_lower") THEN
    1217         564 :          DO i = 1, n - 1
    1218        4716 :             a(i + 1:n, i) = -a(i, i + 1:n)
    1219             :          END DO
    1220             :       ELSE
    1221           0 :          CPABORT("Invalid option <"//TRIM(option)//"> was specified for parameter #2")
    1222             :       END IF
    1223             : 
    1224       16473 :    END SUBROUTINE symmetrize_matrix
    1225             : 
    1226             : ! **************************************************************************************************
    1227             : !> \brief   Set the matrix a to be a unit matrix.
    1228             : !> \param a ...
    1229             : !> \date    16.10.1998
    1230             : !> \author  MK
    1231             : !> \version 1.0
    1232             : ! **************************************************************************************************
    1233       36239 :    PURE SUBROUTINE unit_matrix_d(a)
    1234             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
    1235             : 
    1236     6702115 :       a(:, :) = 0.0_dp
    1237       36239 :       CALL set_diag(a, 1.0_dp)
    1238             : 
    1239       36239 :    END SUBROUTINE unit_matrix_d
    1240             : 
    1241             : ! **************************************************************************************************
    1242             : !> \brief ...
    1243             : !> \param a ...
    1244             : ! **************************************************************************************************
    1245           0 :    PURE SUBROUTINE unit_matrix_z(a)
    1246             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT)   :: a
    1247             : 
    1248           0 :       a(:, :) = (0.0_dp, 0.0_dp)
    1249           0 :       CALL set_diag(a, (1.0_dp, 0.0_dp))
    1250             : 
    1251           0 :    END SUBROUTINE unit_matrix_z
    1252             : 
    1253             : ! **************************************************************************************************
    1254             : !> \brief   Calculation of the vector product c = a x b.
    1255             : !> \param a ...
    1256             : !> \param b ...
    1257             : !> \return ...
    1258             : !> \date    16.10.1998
    1259             : !> \author  MK
    1260             : !> \version 1.0
    1261             : ! **************************************************************************************************
    1262        7659 :    PURE FUNCTION vector_product(a, b) RESULT(c)
    1263             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a, b
    1264             :       REAL(KIND=dp), DIMENSION(3)                        :: c
    1265             : 
    1266        7659 :       c(1) = a(2)*b(3) - a(3)*b(2)
    1267        7659 :       c(2) = a(3)*b(1) - a(1)*b(3)
    1268        7659 :       c(3) = a(1)*b(2) - a(2)*b(1)
    1269             : 
    1270        7659 :    END FUNCTION vector_product
    1271             : 
    1272             : ! **************************************************************************************************
    1273             : !> \brief computes the greatest common divisor of two number
    1274             : !> \param a ...
    1275             : !> \param b ...
    1276             : !> \return ...
    1277             : !> \author Joost VandeVondele
    1278             : ! **************************************************************************************************
    1279      958424 :    ELEMENTAL FUNCTION gcd(a, b)
    1280             :       INTEGER, INTENT(IN)                                :: a, b
    1281             :       INTEGER                                            :: gcd
    1282             : 
    1283             :       INTEGER                                            :: aa, ab, l, rem, s
    1284             : 
    1285      958424 :       aa = ABS(a)
    1286      958424 :       ab = ABS(b)
    1287      958424 :       IF (aa < ab) THEN
    1288             :          s = aa
    1289             :          l = ab
    1290             :       ELSE
    1291      938634 :          s = ab
    1292      938634 :          l = aa
    1293             :       END IF
    1294      958424 :       IF (s .NE. 0) THEN
    1295             :          DO
    1296      958424 :             rem = MOD(l, s)
    1297      958424 :             IF (rem == 0) EXIT
    1298             :             l = s
    1299      958424 :             s = rem
    1300             :          END DO
    1301             :          GCD = s
    1302             :       ELSE
    1303             :          GCD = l
    1304             :       END IF
    1305      958424 :    END FUNCTION gcd
    1306             : 
    1307             : ! **************************************************************************************************
    1308             : !> \brief computes the least common multiplier of two numbers
    1309             : !> \param a ...
    1310             : !> \param b ...
    1311             : !> \return ...
    1312             : !> \author Joost VandeVondele
    1313             : ! **************************************************************************************************
    1314      459367 :    ELEMENTAL FUNCTION lcm(a, b)
    1315             :       INTEGER, INTENT(IN)                                :: a, b
    1316             :       INTEGER                                            :: lcm
    1317             : 
    1318             :       INTEGER                                            :: tmp
    1319             : 
    1320      459367 :       tmp = gcd(a, b)
    1321      459367 :       IF (tmp == 0) THEN
    1322             :          lcm = 0
    1323             :       ELSE
    1324             :          ! could still overflow if the true lcm is larger than maxint
    1325      459367 :          lcm = ABS((a/tmp)*b)
    1326             :       END IF
    1327      459367 :    END FUNCTION lcm
    1328             : 
    1329             : ! **************************************************************************************************
    1330             : !> \brief computes the exponential integral
    1331             : !>      Ei(x) = Int(exp(-x*t)/t,t=1..infinity)  x>0
    1332             : !> \param x ...
    1333             : !> \return ...
    1334             : !> \author JGH (adapted from Numerical recipies)
    1335             : ! **************************************************************************************************
    1336           0 :    FUNCTION ei(x)
    1337             :       REAL(dp)                                           :: x, ei
    1338             : 
    1339             :       INTEGER, PARAMETER                                 :: maxit = 100
    1340             :       REAL(dp), PARAMETER                                :: eps = EPSILON(0.0_dp), &
    1341             :                                                             fpmin = TINY(0.0_dp)
    1342             : 
    1343             :       INTEGER                                            :: k
    1344             :       REAL(dp)                                           :: fact, prev, sum1, term
    1345             : 
    1346           0 :       IF (x <= 0._dp) THEN
    1347           0 :          CPABORT("Invalid argument")
    1348             :       END IF
    1349             : 
    1350           0 :       IF (x < fpmin) THEN
    1351           0 :          ei = LOG(x) + euler
    1352           0 :       ELSE IF (x <= -LOG(EPS)) THEN
    1353             :          sum1 = 0._dp
    1354             :          fact = 1._dp
    1355           0 :          DO k = 1, maxit
    1356           0 :             fact = fact*x/REAL(k, dp)
    1357           0 :             term = fact/REAL(k, dp)
    1358           0 :             sum1 = sum1 + term
    1359           0 :             IF (term < eps*sum1) EXIT
    1360             :          END DO
    1361           0 :          ei = sum1 + LOG(x) + euler
    1362             :       ELSE
    1363             :          sum1 = 0._dp
    1364             :          term = 1._dp
    1365           0 :          DO k = 1, maxit
    1366           0 :             prev = term
    1367           0 :             term = term*REAL(k, dp)/x
    1368           0 :             IF (term < eps) EXIT
    1369           0 :             IF (term < prev) THEN
    1370           0 :                sum1 = sum1 + term
    1371             :             ELSE
    1372           0 :                sum1 = sum1 - prev
    1373           0 :                EXIT
    1374             :             END IF
    1375             :          END DO
    1376           0 :          ei = EXP(x)*(1._dp + sum1)/x
    1377             :       END IF
    1378             : 
    1379           0 :    END FUNCTION ei
    1380             : 
    1381             : ! **************************************************************************************************
    1382             : !> \brief computes the exponential integral
    1383             : !>      En(x) = Int(exp(-x*t)/t^n,t=1..infinity)  x>0, n=0,1,..
    1384             : !>      Note: Ei(-x) = -E1(x)
    1385             : !> \param n ...
    1386             : !> \param x ...
    1387             : !> \return ...
    1388             : !> \par History
    1389             : !>      05.2007 Created
    1390             : !> \author Manuel Guidon (adapted from Numerical recipies)
    1391             : ! **************************************************************************************************
    1392   215454540 :    ELEMENTAL IMPURE FUNCTION expint(n, x)
    1393             :       INTEGER, INTENT(IN)                                :: n
    1394             :       REAL(dp), INTENT(IN)                               :: x
    1395             :       REAL(dp)                                           :: expint
    1396             : 
    1397             :       INTEGER, PARAMETER                                 :: maxit = 100
    1398             :       REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
    1399             :          fpmin = TINY(0.0_dp)
    1400             : 
    1401             :       INTEGER                                            :: i, ii, nm1
    1402             :       REAL(dp)                                           :: a, b, c, d, del, fact, h, psi
    1403             : 
    1404   215454540 :       nm1 = n - 1
    1405             : 
    1406   215454540 :       IF (n .LT. 0 .OR. x .LT. 0.0_dp .OR. (x .EQ. 0.0_dp .AND. (n .EQ. 0 .OR. n .EQ. 1))) THEN
    1407           0 :          CPABORT("Invalid argument")
    1408   215454540 :       ELSE IF (n .EQ. 0) THEN !Special case.
    1409           0 :          expint = EXP(-x)/x
    1410   215454540 :       ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
    1411           0 :          expint = 1.0_dp/nm1
    1412   215454540 :       ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
    1413   178113139 :          b = x + n
    1414   178113139 :          c = 1.0_dp/FPMIN
    1415   178113139 :          d = 1.0_dp/b
    1416   178113139 :          h = d
    1417  5315048305 :          DO i = 1, MAXIT
    1418  5315048305 :             a = -i*(nm1 + i)
    1419  5315048305 :             b = b + 2.0_dp
    1420  5315048305 :             d = 1.0_dp/(a*d + b)
    1421  5315048305 :             c = b + a/c
    1422  5315048305 :             del = c*d
    1423  5315048305 :             h = h*del
    1424  5315048305 :             IF (ABS(del - 1.0_dp) .LT. EPS) THEN
    1425   178113139 :                expint = h*EXP(-x)
    1426   178113139 :                RETURN
    1427             :             END IF
    1428             :          END DO
    1429           0 :          CPABORT("continued fraction failed in expint")
    1430             :       ELSE !Evaluate series.
    1431    37341401 :          IF (nm1 .NE. 0) THEN !Set first term.
    1432           0 :             expint = 1.0_dp/nm1
    1433             :          ELSE
    1434    37341401 :             expint = -LOG(x) - euler
    1435             :          END IF
    1436             :          fact = 1.0_dp
    1437   339896924 :          DO i = 1, MAXIT
    1438   339896924 :             fact = -fact*x/i
    1439   339896924 :             IF (i .NE. nm1) THEN
    1440   339896924 :                del = -fact/(i - nm1)
    1441             :             ELSE
    1442             :                psi = -euler !Compute I(n).
    1443           0 :                DO ii = 1, nm1
    1444           0 :                   psi = psi + 1.0_dp/ii
    1445             :                END DO
    1446           0 :                del = fact*(-LOG(x) + psi)
    1447             :             END IF
    1448   339896924 :             expint = expint + del
    1449   339896924 :             IF (ABS(del) .LT. ABS(expint)*EPS) RETURN
    1450             :          END DO
    1451           0 :          CPABORT("series failed in expint")
    1452             :       END IF
    1453             : 
    1454             :    END FUNCTION expint
    1455             : 
    1456             : ! **************************************************************************************************
    1457             : !> \brief  Jacobi matrix diagonalization. The eigenvalues are returned in
    1458             : !>         vector d and the eigenvectors are returned in matrix v in ascending
    1459             : !>         order.
    1460             : !>
    1461             : !> \param a ...
    1462             : !> \param d ...
    1463             : !> \param v ...
    1464             : !> \par History
    1465             : !>         - Creation (20.11.98, Matthias Krack)
    1466             : ! **************************************************************************************************
    1467       32049 :    SUBROUTINE jacobi(a, d, v)
    1468             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
    1469             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: d
    1470             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: v
    1471             : 
    1472             :       INTEGER                                            :: n
    1473             : 
    1474       32049 :       n = SIZE(d(:))
    1475             : 
    1476             :       ! Diagonalize matrix a
    1477       32049 :       CALL diag(n, a, d, v)
    1478             : 
    1479             :       ! Sort eigenvalues and eigenvector in ascending order
    1480       32049 :       CALL eigsrt(n, d, v)
    1481             : 
    1482       32049 :    END SUBROUTINE jacobi
    1483             : 
    1484             : ! **************************************************************************************************
    1485             : !> \brief  Diagonalize matrix a. The eigenvalues are returned in vector d
    1486             : !>         and the eigenvectors are returned in matrix v.
    1487             : !>
    1488             : !> \param n matrix/vector extent (problem size)
    1489             : !> \param a matrix to be diagonalised
    1490             : !> \param d vector of eigenvalues
    1491             : !> \param v matrix of eigenvectors
    1492             : !> \par History
    1493             : !>         - Creation (20.11.98, Matthias Krack)
    1494             : ! **************************************************************************************************
    1495       32049 :    SUBROUTINE diag(n, a, d, v)
    1496             :       INTEGER, INTENT(IN)                                :: n
    1497             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: a
    1498             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: d
    1499             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: v
    1500             : 
    1501             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'diag'
    1502             :       REAL(KIND=dp), PARAMETER                           :: a_eps = 1.0E-10_dp, d_eps = 1.0E-3_dp
    1503             : 
    1504             :       INTEGER                                            :: handle, i, ip, iq
    1505             :       REAL(KIND=dp)                                      :: a_max, apq, c, d_min, dip, diq, g, h, s, &
    1506             :                                                             t, tau, theta, tresh
    1507       32049 :       REAL(KIND=dp), DIMENSION(n)                        :: b, z
    1508             : 
    1509       32049 :       CALL timeset(routineN, handle)
    1510             : 
    1511       32049 :       a_max = 0.0_dp
    1512      117162 :       DO ip = 1, n - 1
    1513      935906 :          a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
    1514      117162 :          b(ip) = a(ip, ip) ! get_diag(a)
    1515             :       END DO
    1516       32049 :       b(n) = a(n, n)
    1517             : 
    1518       32049 :       CALL unit_matrix(v)
    1519             : 
    1520             :       ! Go for 50 iterations
    1521      155595 :       DO i = 1, 50
    1522      953984 :          d = b
    1523     1109579 :          d_min = MAX(d_eps, MINVAL(ABS(b)))
    1524      155595 :          IF (a_max < a_eps*d_min) THEN
    1525       32049 :             CALL timestop(handle)
    1526             :             RETURN
    1527             :          END IF
    1528      194091 :          tresh = MERGE(a_max, 0.0_dp, (i < 4))
    1529      804773 :          z = 0.0_dp
    1530      681227 :          DO ip = 1, n - 1
    1531     6203136 :             DO iq = ip + 1, n
    1532     5521909 :                dip = d(ip)
    1533     5521909 :                diq = d(iq)
    1534     5521909 :                apq = a(ip, iq)
    1535     5521909 :                g = 100.0_dp*ABS(apq)
    1536     6079590 :                IF (tresh < ABS(apq)) THEN
    1537     3228562 :                   h = diq - dip
    1538     3228562 :                   IF ((ABS(h) + g) .NE. ABS(h)) THEN
    1539     2446652 :                      theta = 0.5_dp*h/apq
    1540     2446652 :                      t = 1.0_dp/(ABS(theta) + SQRT(1.0_dp + theta**2))
    1541     2446652 :                      IF (theta < 0.0_dp) t = -t
    1542             :                   ELSE
    1543      781910 :                      t = apq/h
    1544             :                   END IF
    1545     3228562 :                   c = 1.0_dp/SQRT(1.0_dp + t**2)
    1546     3228562 :                   s = t*c
    1547     3228562 :                   tau = s/(1.0_dp + c)
    1548     3228562 :                   h = t*apq
    1549     3228562 :                   z(ip) = z(ip) - h
    1550     3228562 :                   z(iq) = z(iq) + h
    1551     3228562 :                   d(ip) = dip - h
    1552     3228562 :                   d(iq) = diq + h
    1553     3228562 :                   a(ip, iq) = 0.0_dp
    1554     3228562 :                   CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
    1555     3228562 :                   CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
    1556     3228562 :                   CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
    1557     3228562 :                   CALL jrotate(v(:, ip), v(:, iq), s, tau)
    1558             :                ELSE IF ((4 < i) .AND. &
    1559           0 :                         ((ABS(dip) + g) == ABS(dip)) .AND. &
    1560     2293347 :                         ((ABS(diq) + g) == ABS(diq))) THEN
    1561           0 :                   a(ip, iq) = 0.0_dp
    1562             :                END IF
    1563             :             END DO
    1564             :          END DO
    1565      804773 :          b = b + z
    1566             :          a_max = 0.0_dp
    1567      681227 :          DO ip = 1, n - 1
    1568     6760817 :             a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
    1569             :          END DO
    1570             :       END DO
    1571           0 :       WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
    1572             : 
    1573           0 :       CALL timestop(handle)
    1574             : 
    1575             :    END SUBROUTINE diag
    1576             : 
    1577             : ! **************************************************************************************************
    1578             : !> \brief  Perform a Jacobi rotation of the vectors a and b.
    1579             : !>
    1580             : !> \param a ...
    1581             : !> \param b ...
    1582             : !> \param ss ...
    1583             : !> \param tt ...
    1584             : !> \par History
    1585             : !>         - Creation (20.11.98, Matthias Krack)
    1586             : ! **************************************************************************************************
    1587    12914248 :    PURE SUBROUTINE jrotate(a, b, ss, tt)
    1588             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: a, b
    1589             :       REAL(KIND=dp), INTENT(IN)                          :: ss, tt
    1590             : 
    1591             :       REAL(KIND=dp)                                      :: u, v
    1592             : 
    1593    12914248 :       u = 1.0_dp - ss*tt
    1594    12914248 :       v = ss/u
    1595             : 
    1596   220663620 :       a = a*u - b*ss
    1597   220663620 :       b = b*(u + ss*v) + a*v
    1598             : 
    1599    12914248 :    END SUBROUTINE jrotate
    1600             : 
    1601             : ! **************************************************************************************************
    1602             : !> \brief Sort the values in vector d in ascending order and swap the
    1603             : !>        corresponding columns of matrix v.
    1604             : !>
    1605             : !> \param n ...
    1606             : !> \param d ...
    1607             : !> \param v ...
    1608             : !> \par History
    1609             : !>         - Creation (20.11.98, Matthias Krack)
    1610             : ! **************************************************************************************************
    1611       32049 :    SUBROUTINE eigsrt(n, d, v)
    1612             :       INTEGER, INTENT(IN)                                :: n
    1613             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: d
    1614             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: v
    1615             : 
    1616             :       INTEGER                                            :: i, j
    1617             : 
    1618      117162 :       DO i = 1, n - 1
    1619      935906 :          j = SUM(MINLOC(d(i:n))) + i - 1
    1620      117162 :          IF (j /= i) THEN
    1621       43155 :             CALL swap(d(i), d(j))
    1622       43155 :             CALL swap(v(:, i), v(:, j))
    1623             :          END IF
    1624             :       END DO
    1625             : 
    1626       32049 :    END SUBROUTINE eigsrt
    1627             : 
    1628             : ! **************************************************************************
    1629             : !> \brief Swap two scalars
    1630             : !>
    1631             : !> \param a ...
    1632             : !> \param b ...
    1633             : !> \par History
    1634             : !>         - Creation (20.11.98, Matthias Krack)
    1635             : ! **************************************************************************************************
    1636       43155 :    ELEMENTAL SUBROUTINE swap_scalar(a, b)
    1637             :       REAL(KIND=dp), INTENT(INOUT)                       :: a, b
    1638             : 
    1639             :       REAL(KIND=dp)                                      :: c
    1640             : 
    1641       43155 :       c = a
    1642       43155 :       a = b
    1643       43155 :       b = c
    1644             : 
    1645       43155 :    END SUBROUTINE swap_scalar
    1646             : 
    1647             : ! **************************************************************************
    1648             : !> \brief Swap two vectors
    1649             : !>
    1650             : !> \param a ...
    1651             : !> \param b ...
    1652             : !> \par History
    1653             : !>         - Creation (20.11.98, Matthias Krack)
    1654             : ! **************************************************************************************************
    1655       43155 :    SUBROUTINE swap_vector(a, b)
    1656             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: a, b
    1657             : 
    1658             :       INTEGER                                            :: i, n
    1659             :       REAL(KIND=dp)                                      :: c
    1660             : 
    1661       43155 :       n = SIZE(a)
    1662             : 
    1663       43155 :       IF (n /= SIZE(b)) THEN
    1664           0 :          CPABORT("Check the array bounds of the parameters")
    1665             :       END IF
    1666             : 
    1667      929589 :       DO i = 1, n
    1668      886434 :          c = a(i)
    1669      886434 :          a(i) = b(i)
    1670      929589 :          b(i) = c
    1671             :       END DO
    1672             : 
    1673       43155 :    END SUBROUTINE swap_vector
    1674             : 
    1675             : ! **************************************************************************************************
    1676             : !> \brief - compute a truncation radius for the shortrange operator
    1677             : !> \param eps target accuracy!> \param omg screening parameter
    1678             : !> \param omg ...
    1679             : !> \param r_cutoff cutoff radius
    1680             : !> \par History
    1681             : !>      10.2012 created [Hossein Banihashemian]
    1682             : !>      05.2019 moved here from hfx_types (A. Bussy)
    1683             : !> \author Hossein Banihashemian
    1684             : ! **************************************************************************************************
    1685          64 :    SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
    1686             :       IMPLICIT NONE
    1687             : 
    1688             :       REAL(dp), INTENT(in)  :: eps, omg
    1689             :       REAL(dp), INTENT(out) :: r_cutoff
    1690             : 
    1691             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'erfc_cutoff'
    1692             : 
    1693             :       REAL(dp), PARAMETER :: abstol = 1E-10_dp, soltol = 1E-16_dp
    1694             :       REAL(dp) :: r0, f0, fprime0, delta_r
    1695             :       INTEGER :: iter, handle
    1696             :       INTEGER, PARAMETER :: iterMAX = 1000
    1697             : 
    1698          64 :       CALL timeset(routineN, handle)
    1699             : 
    1700             :       ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
    1701          64 :       r0 = SQRT(-LOG(eps*omg*10**2))/omg
    1702          64 :       CALL eval_transc_func(r0, eps, omg, f0, fprime0)
    1703             : 
    1704         630 :       DO iter = 1, iterMAX
    1705         630 :          delta_r = f0/fprime0
    1706         630 :          r0 = r0 - delta_r
    1707         630 :          CALL eval_transc_func(r0, eps, omg, f0, fprime0)
    1708         630 :          IF (ABS(delta_r) .LT. abstol .OR. ABS(f0) .LT. soltol) EXIT
    1709             :       END DO
    1710          64 :       CPASSERT(iter <= itermax)
    1711          64 :       r_cutoff = r0
    1712             : 
    1713          64 :       CALL timestop(handle)
    1714             :    CONTAINS
    1715             : ! **************************************************************************************************
    1716             : !> \brief ...
    1717             : !> \param r ...
    1718             : !> \param eps ...
    1719             : !> \param omega ...
    1720             : !> \param fn ...
    1721             : !> \param df ...
    1722             : ! **************************************************************************************************
    1723         694 :       ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
    1724             :       REAL(dp), INTENT(in)                               :: r, eps, omega
    1725             :       REAL(dp), INTENT(out)                              :: fn, df
    1726             : 
    1727             :       REAL(dp)                                           :: qr
    1728             : 
    1729         694 :          qr = omega*r
    1730         694 :          fn = ERFC(qr) - r*eps
    1731         694 :          df = -2.0_dp*oorootpi*omega*EXP(-qr**2) - eps
    1732         694 :       END SUBROUTINE eval_transc_func
    1733             :    END SUBROUTINE erfc_cutoff
    1734             : 
    1735             : ! **************************************************************************************************
    1736             : !> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
    1737             : !> \param matrix Hermitian matrix is preserved
    1738             : !> \param eigenvectors ...
    1739             : !> \param eigenvalues ...
    1740             : !> \author A. Bussy
    1741             : ! **************************************************************************************************
    1742       55195 :    SUBROUTINE diag_complex(matrix, eigenvectors, eigenvalues)
    1743             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN)      :: matrix
    1744             :       COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT)     :: eigenvectors
    1745             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigenvalues
    1746             : 
    1747             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'diag_complex'
    1748             : 
    1749       55195 :       COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE        :: work
    1750             :       INTEGER                                            :: handle, info, liwork, lrwork, lwork, n
    1751       55195 :       INTEGER, DIMENSION(:), ALLOCATABLE                 :: iwork
    1752       55195 :       REAL(KIND=dp), DIMENSION(:), ALLOCATABLE           :: rwork
    1753             : 
    1754       55195 :       CALL timeset(routineN, handle)
    1755             : 
    1756       55195 :       IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
    1757             :       ! IF (MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) > 1e-14_dp) CPABORT("Expected hermitian matrix")
    1758             : 
    1759       55195 :       n = SIZE(matrix, 1)
    1760       55195 :       ALLOCATE (iwork(1), rwork(1), work(1))
    1761             : 
    1762             :       ! work space query
    1763       55195 :       lwork = -1
    1764       55195 :       lrwork = -1
    1765       55195 :       liwork = -1
    1766             : 
    1767       55195 :       CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
    1768             : 
    1769       55195 :       lwork = CEILING(REAL(work(1), KIND=dp))
    1770       55195 :       lrwork = CEILING(rwork(1))
    1771       55195 :       liwork = iwork(1)
    1772             : 
    1773       55195 :       DEALLOCATE (iwork, rwork, work)
    1774      386365 :       ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
    1775     2062481 :       eigenvectors(:, :) = matrix(:, :)
    1776             : 
    1777             :       ! final diagonalization
    1778       55195 :       CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
    1779             : 
    1780       55195 :       DEALLOCATE (iwork, rwork, work)
    1781             : 
    1782       55195 :       IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
    1783             : 
    1784       55195 :       CALL timestop(handle)
    1785             : 
    1786       55195 :    END SUBROUTINE diag_complex
    1787             : 
    1788             : ! **************************************************************************************************
    1789             : !> \brief Helper routine for diagonalizing anti symmetric matrices
    1790             : !> \param matrix ...
    1791             : !> \param evecs ...
    1792             : !> \param evals ...
    1793             : ! **************************************************************************************************
    1794        3539 :    SUBROUTINE diag_antisym(matrix, evecs, evals)
    1795             :       REAL(dp), DIMENSION(:, :)                          :: matrix
    1796             :       COMPLEX(dp), DIMENSION(:, :)                       :: evecs
    1797             :       COMPLEX(dp), DIMENSION(:)                          :: evals
    1798             : 
    1799             :       COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :)          :: matrix_c
    1800             :       INTEGER                                            :: n
    1801             :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: eigenvalues
    1802             : 
    1803        3539 :       IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
    1804             :       ! IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 1e-14_dp) CPABORT("Expected anti-symmetric matrix")
    1805             : 
    1806        3539 :       n = SIZE(matrix, 1)
    1807       21234 :       ALLOCATE (matrix_c(n, n), eigenvalues(n))
    1808             : 
    1809      235869 :       matrix_c(:, :) = CMPLX(0.0_dp, -matrix, kind=dp)
    1810        3539 :       CALL diag_complex(matrix_c, evecs, eigenvalues)
    1811       27874 :       evals = CMPLX(0.0_dp, eigenvalues, kind=dp)
    1812             : 
    1813        3539 :       DEALLOCATE (matrix_c, eigenvalues)
    1814        3539 :    END SUBROUTINE diag_antisym
    1815             : 
    1816             : END MODULE mathlib

Generated by: LCOV version 1.15