math.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module math_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use loct_math_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                     &
    is_close,                   &
    diagonal_matrix,            &
    ylmr_cmplx,                 &
    ylmr_real,                  &
    weights,                    &
    factorial,                  &
    hermite,                    &
    dcross_product,             &
    zcross_product,             &
    zdcross_product,            &
    dzcross_product,            &
    ddelta,                     &
    member,                     &
    make_idx_set,               &
    interpolation_coefficients, &
    interpolate,                &
    even,                       &
    odd,                        &
    cartesian2hyperspherical,   &
    hyperspherical2cartesian,   &
    hypersphere_grad_matrix,    &
    pad,                        &
    pad_pow2,                   &
    log2,                       &
    exponential,                &
    phi1,                       &
    phi2,                       &
    square_root,                &
    is_prime,                   &
    is_diagonal,                &
    generate_rotation_matrix,   &
    numder_ridders,             &
    generalized_laguerre_polynomial, &
    upper_triangular_to_hermitian,   &
    lower_triangular_to_hermitian,   &
    dsymmetrize_matrix,              &
    dzero_small_elements_matrix,     &
    convert_to_base,                 &
    convert_from_base,               &
    clamp
 
 
  interface is_close
    module procedure dis_close_scalar, zis_close_scalar
  end interface
 
  interface diagonal_matrix
    module procedure idiagonal_matrix, ddiagonal_matrix, zdiagonal_matrix
  end interface diagonal_matrix
 
  interface upper_triangular_to_hermitian
    module procedure dupper_triangular_to_hermitian, zupper_triangular_to_hermitian
  end interface upper_triangular_to_hermitian
 
  interface lower_triangular_to_hermitian
    module procedure dlower_triangular_to_hermitian, zlower_triangular_to_hermitian
  end interface lower_triangular_to_hermitian
 
 
  ! ------------------------------------------------------------------------------
  interface interpolate
    module procedure dinterpolate_0, dinterpolate_1, dinterpolate_2
    module procedure zinterpolate_0, zinterpolate_1, zinterpolate_2
  end interface interpolate
  ! ------------------------------------------------------------------------------
 
  interface log2
    module procedure dlog2, ilog2, llog2
  end interface log2
 
  interface pad
    module procedure pad4, pad8, pad48, pad88
  end interface pad
 
contains
 
  elemental logical function dis_close_scalar(x, y, rtol, atol)
    real(real64), intent(in) :: x, y
    real(real64), optional, intent(in) :: rtol
    real(real64), optional, intent(in) :: atol
    real(real64) :: atol_, rtol_
 
    atol_ = optional_default(atol, 1.e-8_real64)
    rtol_ = optional_default(rtol, 1.e-5_real64)
 
    dis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
  end function dis_close_scalar
 
  elemental logical function zis_close_scalar(x, y, rtol, atol)
    complex(real64), intent(in) :: x, y
    real(real64), optional, intent(in) :: rtol
    real(real64), optional, intent(in) :: atol
    real(real64) :: atol_, rtol_
 
    atol_ = optional_default(atol, 1.e-8_real64)
    rtol_ = optional_default(rtol, 1.e-5_real64)
 
    zis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
  end function zis_close_scalar
 
 
  ! ---------------------------------------------------------
  pure function idiagonal_matrix(dim, diag) result(matrix)
    integer, intent(in) :: dim
    integer, intent(in) :: diag
    integer :: matrix(dim, dim)
 
    integer :: ii
 
    matrix = 0
    do ii = 1, dim
      matrix(ii, ii) = diag
    end do
 
  end function idiagonal_matrix
 
  ! ---------------------------------------------------------
  recursive function hermite(n, x) result (h)
    integer, intent(in) :: n
    real(real64),   intent(in) :: x
 
    real(real64) :: h
 
    ! no push_sub, called too frequently
 
    if (n <= 0) then
      h = m_one
    elseif (n == 1) then
      h = m_two*x
    else
      h = m_two*x*hermite(n-1,x) - m_two*(n-1)*hermite(n-2,x)
    end if
 
  end function hermite
 
 
  ! ---------------------------------------------------------
  recursive function factorial (n) RESULT (fac)
    integer, intent(in) :: n
 
    integer :: fac
 
    ! no push_sub for recursive function
 
    if (n <= 1) then
      fac = 1
    else
      fac = n * factorial(n-1)
    end if
  end function factorial
 
  ! ---------------------------------------------------------
  subroutine ylmr_cmplx(xx, li, mi, ylm)
    real(real64),   intent(in) :: xx(3)
    integer, intent(in) :: li, mi
    complex(real64),   intent(out) :: ylm
 
    integer :: i
    real(real64) :: dx, dy, dz, r, plm, cosm, sinm, cosmm1, sinmm1, cosphi, sinphi
    real(real64),   parameter :: tiny = 1.e-15_real64
 
    !   no push_sub: this routine is called too many times
 
    ! if l=0, no calculations are required
    if (li == 0) then
      ylm = cmplx(sqrt(m_one/(m_four*m_pi)), m_zero, real64)
      return
    end if
 
    r = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
 
    ! if r=0, direction is undefined => make ylm=0 except for l=0
    if (r <= tiny) then
      ylm = m_z0
      return
    end if
    dx = xx(1)/r
    dy = xx(2)/r
    dz = xx(3)/r
 
    ! get the associated Legendre polynomial (including the normalization factor)
    ! We need to avoid getting outside [-1:1] due to rounding errors
    if(abs(dz) > m_one) dz = sign(m_one, dz)
    plm = loct_legendre_sphplm(li, abs(mi), dz)
 
    ! compute sin(|m|*phi) and cos(|m|*phi)
    r = sqrt(dx*dx + dy*dy)
    if(r > tiny) then
      cosphi = dx/r
      sinphi = dy/r
    else ! In this case the values are ill defined so we choose \phi=0
      cosphi = m_zero
      sinphi = m_zero
    end if
 
    cosm = m_one
    sinm = m_zero
    do i = 1, abs(mi)
      cosmm1 = cosm
      sinmm1 = sinm
      cosm = cosmm1*cosphi - sinmm1*sinphi
      sinm = cosmm1*sinphi + sinmm1*cosphi
    end do
 
    !And now ylm
    ylm = plm*cmplx(cosm, sinm, real64)
 
    if (mi < 0) then
      ylm = conjg(ylm)
      do i = 1, abs(mi)
        ylm = -ylm
      end do
    end if
 
  end subroutine ylmr_cmplx
 
 
  ! ---------------------------------------------------------
  subroutine ylmr_real(xx, li, mi, ylm)
    real(real64),    intent(in)  :: xx(3)
    integer,         intent(in)  :: li, mi
    real(real64),    intent(out) :: ylm
 
    integer, parameter :: lmaxd = 20
    real(real64),   parameter :: tiny = 1.e-15_real64
    integer :: i, ilm0, l, m, mabs
    integer, save :: lmax = -1
 
    real(real64) :: cmi, cosm, cosmm1, cosphi, dphase, dplg, fac, &
      fourpi, plgndr, phase, pll, pmm, pmmp1, sinm, &
      sinmm1, sinphi, rsize, rx, ry, rz, xysize
    real(real64), save :: c(0:(lmaxd+1)*(lmaxd+1))
 
    ! no push_sub: called too frequently
 
    ! evaluate normalization constants once and for all
    if (li > lmax) then
      fourpi = m_four*m_pi
      do l = 0, li
        ilm0 = l*l + l
        do m = 0, l
          fac = (2*l+1)/fourpi
          do i = l - m + 1, l + m
            fac = fac/i
          end do
          c(ilm0 + m) = sqrt(fac)
          ! next line because ylm are real combinations of m and -m
          if (m /= 0) c(ilm0 + m) = c(ilm0 + m)*sqrt(m_two)
          c(ilm0 - m) = c(ilm0 + m)
        end do
      end do
      lmax = li
    end if
 
    ! if l=0, no calculations are required
    if (li == 0) then
      ylm = c(0)
      return
    end if
 
    ! if r=0, direction is undefined => make ylm=0 except for l=0
    rsize = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
    if (rsize < tiny) then
      ylm = m_zero
      return
    end if
 
    rx = xx(1)/rsize
    ry = xx(2)/rsize
    rz = xx(3)/rsize
 
    ! explicit formulas for l=1 and l=2
    if (li == 1) then
      select case (mi)
      case (-1)
        ylm = (-c(1))*ry
      case (0)
        ylm = c(2)*rz
      case (1)
        ylm = (-c(3))*rx
      end select
      return
    end if
 
    if (li == 2) then
      select case (mi)
      case (-2)
        ylm = c(4)*6.0_real64*rx*ry
      case (-1)
        ylm = (-c(5))*m_three*ry*rz
      case (0)
        ylm = c(6)*m_half*(m_three*rz*rz - m_one)
      case (1)
        ylm = (-c(7))*m_three*rx*rz
      case (2)
        ylm = c(8)*m_three*(rx*rx - ry*ry)
      end select
      return
    end if
 
    ! general algorithm based on routine plgndr of numerical recipes
    mabs = abs(mi)
    xysize = sqrt(max(rx*rx + ry*ry, tiny))
    cosphi = rx/xysize
    sinphi = ry/xysize
    cosm = m_one
    sinm = m_zero
    do m = 1, mabs
      cosmm1 = cosm
      sinmm1 = sinm
      cosm = cosmm1*cosphi - sinmm1*sinphi
      sinm = cosmm1*sinphi + sinmm1*cosphi
    end do
 
    if (mi < 0) then
      phase = sinm
      dphase = mabs*cosm
    else
      phase = cosm
      dphase = (-mabs)*sinm
    end if
 
    pmm = m_one
    fac = m_one
 
    if (mabs > m_zero) then
      do i = 1, mabs
        pmm = (-pmm)*fac*xysize
        fac = fac + m_two
      end do
    end if
 
    if (li == mabs) then
      plgndr = pmm
      dplg = (-li)*rz*pmm/(xysize**2)
    else
      pmmp1 = rz*(2*mabs + 1)*pmm
      if (li == mabs + 1) then
        plgndr = pmmp1
        dplg = -((li*rz*pmmp1 - (mabs + li)*pmm)/(xysize**2))
      else
        do l = mabs + 2, li
          pll = (rz*(2*l - 1)*pmmp1 - (l + mabs - 1)*pmm)/(l - mabs)
          pmm = pmmp1
          pmmp1 = pll
        end do
        plgndr = pll
        dplg = -((li*rz*pll - (l + mabs - 1)*pmm)/(xysize**2))
      end if
    end if
 
    ilm0 = li*li + li
    cmi = c(ilm0 + mi)
    ylm = cmi*plgndr*phase
 
  end subroutine ylmr_real
 
 
  ! ---------------------------------------------------------
  subroutine weights(N, M, cc, side)
    integer,           intent(in)  :: n, m
    real(real64),      intent(out) :: cc(0:,0:,0:)
    integer, optional, intent(in)  :: side
 
    integer :: i, j, k, mn, side_
    real(real64) :: c1, c2, c3, c4, c5, xi
    real(real64), allocatable :: x(:)
 
    push_sub(weights)
 
    safe_allocate(x(0:m))
 
    if (present(side)) then
      side_ = side
    else
      side_ = 0
    end if
 
    select case (side_)
    case (-1)
      ! grid-points for left-side finite-difference formulas on an equi-spaced grid
      mn = m
      x(:) = (/(-i,i=0,mn)/)
    case (+1)
      ! grid-points for right-side finite-difference formulas on an equi-spaced grid
      mn = m
      x(:) = (/(i,i=0,mn)/)
    case default
      ! grid-points for centered finite-difference formulas on an equi-spaced grid
      mn = m/2
      x(:) = (/0,(-i,i,i=1,mn)/)
    end select
 
 
 
    xi = m_zero  ! point at which the approx. are to be accurate
 
    cc = m_zero
    cc(0,0,0) = m_one
 
    c1 = m_one
    c4 = x(0) - xi
 
    do j = 1, m
      mn = min(j,n)
      c2 = m_one
      c5 = c4
      c4 = x(j) - xi
 
      do k = 0, j - 1
        c3 = x(j) - x(k)
        c2 = c2*c3
 
        if (j <= n) cc(k, j - 1, j) = m_zero
        cc(k, j, 0) = c4*cc(k, j - 1, 0)/c3
 
        do i = 1, mn
          cc(k, j, i) = (c4*cc(k, j - 1, i) - i*cc(k, j - 1, i - 1))/c3
        end do
 
        cc(j, j, 0) = -c1*c5*cc(j - 1, j - 1, 0) / c2
      end do
 
      do i = 1, mn
        cc(j, j, i) = c1*(i*cc(j - 1, j - 1, i - 1) - c5*cc(j - 1, j - 1, i))/c2
      end do
 
      c1 = c2
    end do
 
    safe_deallocate_a(x)
    pop_sub(weights)
  end subroutine weights
 
  ! ---------------------------------------------------------
  real(real64) pure function ddelta(i, j)
    integer, intent(in) :: i
    integer, intent(in) :: j
 
    ! no push_sub in pure function
 
    if (i == j) then
      ddelta = m_one
    else
      ddelta = m_zero
    end if
 
  end function ddelta
 
 
  ! ---------------------------------------------------------
  subroutine make_idx_set(n, out, length, in)
    integer,              intent(in)  :: n
    integer, allocatable, intent(out) :: out(:)
    integer,              intent(out) :: length
    integer, optional,    intent(in)  :: in(:)
 
    integer :: i
 
    push_sub(make_idx_set)
 
    if (present(in)) then
      length = size(in, 1)
      safe_allocate(out(1:length))
      out(:) = in(:)
    else
      length = n
      safe_allocate(out(1:length))
      do i = 1, length
        out(i) = i
      end do
    end if
 
    pop_sub(make_idx_set)
  end subroutine make_idx_set
 
 
  ! ---------------------------------------------------------
  logical function member(n, a)
    integer, intent(in) :: n
    integer, intent(in) :: a(:)
 
    integer :: i
 
    ! no push_sub, called too frequently
 
    member = .false.
 
    do i = 1, ubound(a, 1)
      if (a(i) == n) then
        member = .true.
        exit
      end if
    end do
 
  end function member
 
 
  ! ---------------------------------------------------------
  subroutine interpolation_coefficients(nn, xa, xx, cc)
    integer, intent(in)  :: nn
    real(real64),   intent(in)  :: xa(:)
    real(real64),   intent(in)  :: xx
    real(real64),   intent(out) :: cc(:)
 
    integer :: ii, kk
 
    ! no push_sub, called too frequently
 
    do ii = 1, nn
      cc(ii) = m_one
      do kk = 1, nn
        if (kk == ii) cycle
        cc(ii) = cc(ii)*(xx - xa(kk))/(xa(ii) - xa(kk))
      end do
    end do
 
  end subroutine interpolation_coefficients
 
 
  ! ---------------------------------------------------------
  logical pure function even(n)
    integer, intent(in) :: n
 
    even = (mod(n, 2) == 0)
 
  end function even
 
 
  ! ---------------------------------------------------------
  logical pure function odd(n)
    integer, intent(in) :: n
 
    odd = .not. even(n)
 
  end function odd
 
  ! ---------------------------------------------------------
  subroutine cartesian2hyperspherical(x, u)
    real(real64), intent(in)  :: x(:)
    real(real64), intent(out) :: u(:)
 
    integer :: n, k, j
    real(real64) :: sumx2
    real(real64), allocatable :: xx(:)
 
    push_sub(cartesian2hyperspherical)
 
    n = size(x)
    assert(n>1)
    assert(size(u) == n-1)
 
    ! These lines make the code less machine-dependent.
    safe_allocate(xx(1:n))
    do j = 1, n
      if (abs(x(j))<1.0e-8_real64) then
        xx(j) = m_zero
      else
        xx(j) = x(j)
      end if
    end do
 
    u = m_zero
    do k = 1, n-1
      sumx2 = m_zero
      do j = k+1, n
        sumx2 = sumx2 + xx(j)**2
      end do
      if (abs(sumx2) <= m_epsilon .and. abs(xx(k)) <= m_epsilon) exit
      if (k < n - 1) then
        u(k) = atan2(sqrt(sumx2), xx(k))
      else
        u(n-1) = atan2(xx(n), xx(n-1))
      end if
    end do
 
 
    pop_sub(cartesian2hyperspherical)
  end subroutine cartesian2hyperspherical
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine hyperspherical2cartesian(u, x)
    real(real64), intent(in)  :: u(:)
    real(real64), intent(out) :: x(:)
 
    integer :: n, j, k
 
    push_sub(hyperspherical2cartesian)
 
    n = size(x)
    assert(n>1)
    assert(size(u) == n-1)
 
    if (n == 2) then
      x(1) = cos(u(1))
      x(2) = sin(u(1))
    elseif (n == 3) then
      x(1) = cos(u(1))
      x(2) = sin(u(1))*cos(u(2))
      x(3) = sin(u(1))*sin(u(2))
    else
      x(1) = cos(u(1))
      x(2) = sin(u(1))*cos(u(2))
      x(3) = sin(u(1))*sin(u(2))*cos(u(3))
      do j = 4, n - 1
        x(j) = m_one
        do k = 1, j - 1
          x(j) = x(j) * sin(u(k))
        end do
        x(j) = x(j) * cos(u(j))
      end do
      x(n) = m_one
      do k = 1, n - 2
        x(n) = x(n) * sin(u(k))
      end do
      x(n) = x(n) * sin(u(n-1))
    end if
 
    pop_sub(hyperspherical2cartesian)
  end subroutine hyperspherical2cartesian
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine  hypersphere_grad_matrix(grad_matrix, r, x)
    real(real64), intent(out)  :: grad_matrix(:,:)
    real(real64), intent(in)   :: r
    real(real64), intent(in)   :: x(:)
 
    integer :: n, l, m
 
    push_sub(hypersphere_grad_matrix)
 
    n = size(x)+1  ! the dimension of the matrix is (n-1)x(n)
 
    ! --- l=1 ---
    grad_matrix = m_one
    grad_matrix(1,1) = -r*sin(x(1))
    do m = 2, n-1
      grad_matrix(m,1) = m_zero
    end do
 
    ! --- l=2..(n-1) ---
    do l=2, n-1
      do m=1, n-1
        if (m == l) then
          grad_matrix(m,l) = -r*grad_matrix(m,l)*product(sin(x(1:l)))
        elseif (m < l) then
          grad_matrix(m,l) = grad_matrix(m,l)*r*cos(x(m))*cos(x(l))*product(sin(x(1:m-1)))*product(sin(x(m+1:l-1)))
        else
          grad_matrix(m,l) = m_zero
        end if
      end do
    end do
 
    ! --- l=n ---
    do m=1, n-1
      grad_matrix(m,n) = r*cos(x(m))*grad_matrix(m,n)*product(sin(x(1:m-1)))*product(sin(x(m+1:n-1)))
    end do
 
    pop_sub(hypersphere_grad_matrix)
  end subroutine  hypersphere_grad_matrix
 
  ! ---------------------------------------------------------
  integer(int64) pure function pad88(size, blk)
    integer(int64), intent(in) :: size
    integer(int64), intent(in) :: blk
 
    integer(int64) :: mm
 
    mm = mod(size, blk)
    if (mm == 0) then
      pad88 = size
    else
      pad88 = size + blk - mm
    end if
  end function pad88
 
  ! ---------------------------------------------------------
  integer(int64) pure function pad48(size, blk)
    integer,     intent(in) :: size
    integer(int64), intent(in) :: blk
 
    pad48 = pad88(int(size, int64), blk)
  end function pad48
 
  ! ---------------------------------------------------------
  integer(int64) pure function pad8(size, blk)
    integer(int64), intent(in) :: size
    integer, intent(in) :: blk
 
    pad8 = pad88(size, int(blk, int64))
  end function pad8
 
  ! ---------------------------------------------------------
  integer pure function pad4(size, blk)
    integer, intent(in) :: size
    integer, intent(in) :: blk
 
    pad4 = int(pad88(int(size, int64), int(blk, int64)), int32)
  end function pad4
 
  ! ---------------------------------------------------------
 
  integer pure function pad_pow2(size)
    integer, intent(in) :: size
 
    integer :: mm, mm2
 
    mm = size
    pad_pow2 = 1
 
    ! loop below never terminates otherwise! just pick 1 as value.
    if (size == 0) return
 
    ! first we divide by two all the times we can, so we catch exactly
    ! the case when size is already a power of two
    do
      mm2 = mm/2
      if (mm - 2*mm2 /= 0) exit
      pad_pow2 = pad_pow2*2
      mm = mm2
    end do
 
    ! the rest is handled by log2
    if (mm /= 1) then
      pad_pow2 = pad_pow2*2**(ceiling(log(real(mm, real64) )/log(2.0_8)))
    end if
 
  end function pad_pow2
 
  ! -------------------------------------------------------
 
  real(real64) pure function dlog2(xx)
    real(real64), intent(in) :: xx
 
    dlog2 = log(xx)/log(m_two)
  end function dlog2
 
  ! -------------------------------------------------------
 
  integer pure function ilog2(xx)
    integer, intent(in) :: xx
 
    ilog2 = nint(log2(real(xx, real64) ))
  end function ilog2
 
  ! -------------------------------------------------------
 
  integer(int64) pure function llog2(xx)
    integer(int64), intent(in) :: xx
 
    llog2 = nint(log2(real(xx, real64) ), kind=int64)
  end function llog2
 
  ! -------------------------------------------------------
  complex(real64) pure function exponential(z)
    complex(real64), intent(in) :: z
 
    exponential = exp(z)
  end function exponential
 
  ! -------------------------------------------------------
  complex(real64) pure function phi1(z)
    complex(real64), intent(in) :: z
 
    real(real64), parameter :: cut = 0.002_real64
 
    if (abs(z) < cut) then
      phi1 = m_one + m_half*z*(m_one + z/m_three*(m_one + m_fourth*z*(m_one + z/m_five)))
    else
      phi1 = (exp(z) - m_z1)/z
    end if
  end function phi1
 
  ! -------------------------------------------------------
  complex(real64) pure function phi2(z)
    complex(real64), intent(in) :: z
 
    real(real64), parameter :: cut = 0.002_real64
 
    if (abs(z) < cut) then
      phi2 = m_half + z/6.0_real64*(m_one + m_fourth*z*(m_one + z/m_five*(m_one + z/6.0_real64)))
    else
      phi2 = (exp(z) - z - m_z1)/z**2
    end if
  end function phi2
 
  ! -------------------------------------------------------
  real(real64) pure function square_root(x)
    real(real64), intent(in) :: x
 
    square_root = sqrt(x)
  end function square_root
 
  ! -------------------------------------------------------
 
  logical function is_prime(n)
    integer, intent(in) :: n
 
    integer :: i, root
 
    push_sub(is_prime)
 
    if (n < 1) then
      message(1) = "Internal error: is_prime does not take negative numbers."
      call messages_fatal(1)
    end if
    if (n == 1) then
      is_prime = .false.
      pop_sub(is_prime)
      return
    end if
 
    root = nint(sqrt(real(n, real64) ))
    do i = 2, root
      if (mod(n,i) == 0) then
        is_prime = .false.
        pop_sub(is_prime)
        return
      end if
    end do
 
    is_prime = .true.
    pop_sub(is_prime)
  end function is_prime
 
  ! ---------------------------------------------------------
  subroutine generate_rotation_matrix(R, ff, tt)
    real(real64),   intent(out)  :: r(:,:)
    real(real64),   intent(in)   :: ff(:)
    real(real64),   intent(in)   :: tt(:)
 
    integer            :: dim, i, j
    real(real64)       :: th, uv, uu, vv, ft
    real(real64), allocatable :: axis(:), u(:), v(:), f(:), t(:), p(:)
 
    push_sub(generate_rotation_matrix)
 
    dim = size(ff,1)
 
    assert((dim == 2) .or. (dim == 3))
    assert(size(tt,1) == dim)
    assert((size(r,1) == dim) .and. (size(r,2) == dim))
 
    safe_allocate(u(1:dim))
    safe_allocate(v(1:dim))
    safe_allocate(f(1:dim))
    safe_allocate(t(1:dim))
 
 
    !normalize
    f = ff / norm2(ff)
    t = tt / norm2(tt)
 
    ft = dot_product(f,t)
    ! Note : if ff and tt are antiparallel, we get ft = -1 and R = Identity.
    ! This is wrong.
 
    if (abs(ft) < m_one) then
      select case (dim)
      case (2)
        th = acos(ft)
        r(1,1) = cos(th)
        r(1,2) = -sin(th)
 
        r(2,1) = sin(th)
        r(2,2) = cos(th)
 
      case (3)
        if (.false.) then
          !Old implementation
          safe_allocate(axis(1:dim))
          th = acos(ft)
 
          u = f / dot_product(f,f)
          v = t /dot_product(t,t)
 
          axis = dcross_product(u,v)
          axis = axis / norm2(axis)
 
          r(1,1) = cos(th) + axis(1)**2 * (1 - cos(th))
          r(1,2) = axis(1)*axis(2)*(1-cos(th)) + axis(3)*sin(th)
          r(1,3) = axis(1)*axis(3)*(1-cos(th)) - axis(2)*sin(th)
 
          r(2,1) = axis(2)*axis(1)*(1-cos(th)) - axis(3)*sin(th)
          r(2,2) = cos(th) + axis(2)**2 * (1 - cos(th))
          r(2,3) = axis(2)*axis(3)*(1-cos(th)) + axis(1)*sin(th)
 
          r(3,1) = axis(3)*axis(1)*(1-cos(th)) + axis(2)*sin(th)
          r(3,2) = axis(3)*axis(2)*(1-cos(th)) - axis(1)*sin(th)
          r(3,3) = cos(th) + axis(3)**2 * (1 - cos(th))
 
          safe_deallocate_a(axis)
        end if
 
        if (.true.) then
          !Naive implementation
          th = acos(ft)
          u = dcross_product(f,t)
          u = u / norm2(u)
 
          r(1,1) = u(1)**2 + (1-u(1)**2)*cos(th)
          r(1,2) = u(1)*u(2)*(1-cos(th)) - u(3)*sin(th)
          r(1,3) = u(1)*u(3) + u(2)*sin(th)
 
          r(2,1) = u(1)*u(2)*(1-cos(th)) + u(3)*sin(th)
          r(2,2) = u(2)**2 + (1-u(2)**2)*cos(th)
          r(2,3) = u(2)*u(3)*(1-cos(th)) - u(1)*sin(th)
 
          r(3,1) = u(1)*u(3)*(1-cos(th)) - u(2)*sin(th)
          r(3,2) = u(2)*u(3)*(1-cos(th)) + u(1)*sin(th)
          r(3,3) = u(3)**2 + (1-u(3)**2)*cos(th)
        end if
 
        if (.false.) then
          !Fast
          safe_allocate(p(1:dim))
 
          if (abs(f(1)) <= abs(f(2)) .and. abs(f(1)) < abs(f(3))) then
            p = (/m_one, m_zero, m_zero/)
          else if (abs(f(2)) < abs(f(1)) .and. abs(f(2)) <= abs(f(3))) then
            p = (/m_zero, m_one, m_zero/)
          else if (abs(f(3)) <= abs(f(1)) .and. abs(f(3)) < abs(f(2))) then
            p = (/m_zero, m_zero, m_one/)
          end if
 
          u = p - f
          v = p - t
 
          uu = dot_product(u,u)
          vv = dot_product(v,v)
          uv = dot_product(u,v)
 
          do i = 1,3
            do j = 1,3
 
              r(i,j) = ddelta(i,j) - m_two * u(i)*u(j)/uu - m_two * v(i)*v(j)/vv &
                + 4_real64*uv * v(i)*u(j) /(uu*vv)
 
            end do
          end do
 
          safe_deallocate_a(p)
        end if
 
 
      end select
 
    else
 
      r = m_zero
      do i = 1,dim
        r(i,i) = m_one
      end do
 
    end if
 
 
    safe_deallocate_a(u)
    safe_deallocate_a(v)
 
    pop_sub(generate_rotation_matrix)
  end subroutine generate_rotation_matrix
 
  ! ---------------------------------------------------------
  subroutine numder_ridders(x, h, res, err, f)
    implicit none
    real(real64), intent(in)  :: x, h
    real(real64), intent(out) :: res, err
    interface
      subroutine f(x, fx)
        use, intrinsic :: iso_fortran_env
        implicit none
        real(real64), intent(in)    :: x
        real(real64), intent(inout) :: fx
      end subroutine f
    end interface
 
    real(real64), parameter :: con = 1.4_real64, big = 1.0e30_real64, safe = m_two
    integer, parameter :: ntab = 20
    integer :: i, j
    real(real64) :: errt, fac, hh, fx1, fx2
    real(real64), allocatable :: a(:, :)
 
    push_sub(numder_ridders)
 
    if (abs(h) <= m_epsilon) then
      message(1) = "h must be nonzero in numder_ridders"
      call messages_fatal(1)
    end if
 
    safe_allocate(a(1:ntab, 1:ntab))
 
    hh = h
    call f(x+hh, fx1)
    call f(x-hh, fx2)
    a(1,1) = (fx1-fx2) / (m_two*hh)
    err = big
    do i = 2, ntab
      hh = hh / con
      call f(x+hh, fx1)
      call f(x-hh, fx2)
      a(1,i) = (fx1-fx2) / (m_two*hh)
      fac = con**2
      do j = 2, i
        a(j,i) = (a(j-1,i)*fac-a(j-1,i-1)) / (fac-m_one)
        fac = con**2*fac
        errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
        if (errt .le. err) then
          err = errt
          res = a(j,i)
        end if
      end do
      if (abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) then
        safe_deallocate_a(a)
        pop_sub(numder_ridders)
        return
      end if
    end do
 
    safe_deallocate_a(a)
    pop_sub(numder_ridders)
  end subroutine numder_ridders
 
  ! ---------------------------------------------------------
  pure function dzcross_product(a, b) result(c)
    real(real64), intent(in) :: a(:)
    complex(real64), intent(in) :: b(:)
 
    complex(real64) :: c(1:3)
 
    c(1) = a(2)*b(3) - a(3)*b(2)
    c(2) = a(3)*b(1) - a(1)*b(3)
    c(3) = a(1)*b(2) - a(2)*b(1)
 
  end function dzcross_product
 
  ! ---------------------------------------------------------
  pure function zdcross_product(a, b) result(c)
    complex(real64), intent(in) :: a(:)
    real(real64), intent(in) :: b(:)
 
    complex(real64) :: c(1:3)
 
    c(1) = a(2)*b(3) - a(3)*b(2)
    c(2) = a(3)*b(1) - a(1)*b(3)
    c(3) = a(1)*b(2) - a(2)*b(1)
 
  end function zdcross_product
 
 
!*****************************************************************************80
!
!! LM_POLYNOMIAL evaluates Laguerre polynomials Lm(n,m,x).
!
!  First terms:
!
!    M = 0
!
!    Lm(0,0,X) =   1
!    Lm(1,0,X) =  -X   +  1
!    Lm(2,0,X) =   X^2 -  4 X   +  2
!    Lm(3,0,X) =  -X^3 +  9 X^2 -  18 X   +    6
!    Lm(4,0,X) =   X^4 - 16 X^3 +  72 X^2 -   96 X +     24
!    Lm(5,0,X) =  -X^5 + 25 X^4 - 200 X^3 +  600 X^2 -  600 x   +  120
!    Lm(6,0,X) =   X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720
!
!    M = 1
!
!    Lm(0,1,X) =    0
!    Lm(1,1,X) =   -1,
!    Lm(2,1,X) =    2 X - 4,
!    Lm(3,1,X) =   -3 X^2 + 18 X - 18,
!    Lm(4,1,X) =    4 X^3 - 48 X^2 + 144 X - 96
!
!    M = 2
!
!    Lm(0,2,X) =    0
!    Lm(1,2,X) =    0,
!    Lm(2,2,X) =    2,
!    Lm(3,2,X) =   -6 X + 18,
!    Lm(4,2,X) =   12 X^2 - 96 X + 144
!
!    M = 3
!
!    Lm(0,3,X) =    0
!    Lm(1,3,X) =    0,
!    Lm(2,3,X) =    0,
!    Lm(3,3,X) =   -6,
!    Lm(4,3,X) =   24 X - 96
!
!    M = 4
!
!    Lm(0,4,X) =    0
!    Lm(1,4,X) =    0
!    Lm(2,4,X) =    0
!    Lm(3,4,X) =    0
!    Lm(4,4,X) =   24
!
!  Recursion:
!
!    Lm(0,M,X)   = 1
!    Lm(1,M,X)   = (M+1-X)
!
!    if 2 <= N:
!
!      Lm(N,M,X)   = ( (M+2*N-1-X) * Lm(N-1,M,X)
!                   +   (1-M-N)    * Lm(N-2,M,X) ) / N
!
!  Special values:
!
!    For M = 0, the associated Laguerre polynomials Lm(N,M,X) are equal
!    to the Laguerre polynomials L(N,X).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    08 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    National Bureau of Standards, 1964,
!    ISBN: 0-486-61272-4,
!    LC: QA47.A34.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) MM, the number of evaluation points.
!
!    Input, integer ( kind = 4 ) N, the highest order polynomial to compute.
!    Note that polynomials 0 through N will be computed.
!
!    Input, integer ( kind = 4 ) M, the parameter.  M must be nonnegative.
!
!    Input, real ( kind = rk ) X(MM), the evaluation points.
!
!    Output, real ( kind = rk ) CX(MM,0:N), the associated Laguerre polynomials
!    of degrees 0 through N evaluated at the evaluation points.
!
! Taken from
! https://people.sc.fsu.edu/~jburkardt/f_src/laguerre_polynomial/laguerre_polynomial.html
! and adapted to Octopus by N. Tancogne-Dejean
  subroutine generalized_laguerre_polynomial ( np, nn, mm, xx, cx )
    integer, intent(in)  :: np
    integer, intent(in)  :: nn, mm
    real(real64),   intent(in)  :: xx(np)
    real(real64),   intent(out) :: cx(np)
 
    integer ii
    real(real64), allocatable :: cx_tmp(:,:)
 
    assert(mm >= 0)
 
    safe_allocate(cx_tmp(1:np, 0:nn))
 
    if (nn < 0) then
      cx = m_zero
      return
    end if
 
    if (nn == 0) then
      cx(1:np) = m_one
      return
    end if
 
    cx_tmp(1:np,0) = m_one
    cx_tmp(1:np,1) = real(mm + 1, real64)  - xx(1:np)
 
    do ii = 2, nn
      cx_tmp(1:np, ii) = &
        (( real(mm + 2*ii - 1, real64)  - xx(1:np)) * cx_tmp(1:np,ii-1)   &
        + real(-mm    -ii + 1, real64)              * cx_tmp(1:np,ii-2)) / real(ii, real64)
    end do
 
    cx(1:np) = cx_tmp(1:np, nn)
 
    safe_deallocate_a(cx_tmp)
 
  end subroutine generalized_laguerre_polynomial
 
  subroutine dupper_triangular_to_hermitian(nn, aa)
    integer, intent(in)     :: nn
    real(real64),   intent(inout)  :: aa(:, :)
 
    integer :: ii, jj
 
    do jj = 1, nn
      do ii = 1, jj-1
        aa(jj, ii) = aa(ii, jj)
      end do
    end do
  end subroutine dupper_triangular_to_hermitian
 
  subroutine zupper_triangular_to_hermitian(nn, aa)
    integer, intent(in)     :: nn
    complex(real64),   intent(inout)  :: aa(:, :)
 
    integer :: ii, jj
 
    do jj = 1, nn
      do ii = 1, jj-1
        aa(jj, ii) = conjg(aa(ii, jj))
      end do
      aa(jj, jj) = real(aa(jj, jj), real64)
    end do
  end subroutine zupper_triangular_to_hermitian
 
  subroutine dlower_triangular_to_hermitian(nn, aa)
    integer, intent(in)     :: nn
    real(real64),   intent(inout)  :: aa(:, :)
 
    integer :: ii, jj
 
    do jj = 1, nn
      do ii = 1, jj-1
        aa(ii, jj) = aa(jj, ii)
      end do
    end do
  end subroutine dlower_triangular_to_hermitian
 
  subroutine zlower_triangular_to_hermitian(nn, aa)
    integer, intent(in)     :: nn
    complex(real64),   intent(inout)  :: aa(:, :)
 
    integer :: ii, jj
 
    do jj = 1, nn
      do ii = 1, jj-1
        aa(ii, jj) = conjg(aa(jj, ii))
      end do
      aa(jj, jj) = real(aa(jj, jj), real64)
    end do
  end subroutine zlower_triangular_to_hermitian
 
 
  subroutine dsymmetrize_matrix(nn, aa)
    integer, intent(in)     :: nn
    real(real64),   intent(inout)  :: aa(:, :)
 
    integer :: ii, jj
 
    do ii = 1, nn
      do jj = ii+1, nn
        aa(ii, jj) = m_half * (aa(ii, jj) + aa(jj, ii))
        aa(jj, ii) = aa(ii, jj)
      end do
    end do
  end subroutine dsymmetrize_matrix
 
  subroutine dzero_small_elements_matrix(nn, aa, tol)
    integer,        intent(in)     :: nn
    real(real64),   intent(inout)  :: aa(:, :)
    real(real64),   intent(in)     :: tol
 
    where(abs(aa)<tol)
      aa = m_zero
    end where
  end subroutine dzero_small_elements_matrix
 
  ! ---------------------------------------------------------
  logical pure function is_diagonal(dim, matrix, tol)
    integer,        intent(in)    :: dim
    real(real64),   intent(in)    :: matrix(:, :)
    real(real64),   intent(in)    :: tol
 
    integer :: ii, jj
    real(real64) :: max_diag
 
    max_diag = m_zero
    do ii = 1, dim
      max_diag = max(abs(matrix(ii, ii)), max_diag)
    end do
 
    is_diagonal = .true.
    do ii = 1, dim
      do jj = 1, dim
        if (ii == jj) cycle
        if (abs(matrix(ii, jj)) > tol*max_diag) is_diagonal = .false.
      end do
    end do
 
  end function is_diagonal
 
  subroutine convert_to_base(num, base, converted)
    integer, intent(in) :: num
    integer, intent(in) :: base
    integer, intent(out) :: converted(:)
 
    integer :: i, num_, length
 
    ! no PUSH_SUB, called too often
 
    converted = 0
    if (num > 0) then
      ! length of converted number
      length = floor(log(real(num, real64)) / log(real(base, real64))) + 1
      assert(size(converted) >= length)
 
      num_ = num
      ! do conversion by repeated division
      do i = 1, size(converted)
        converted(i) = mod(num_, base)
        num_ = num_ / base
      end do
    end if
 
  end subroutine convert_to_base
 
  subroutine convert_from_base(converted, base, num)
    integer, intent(in) :: converted(:)
    integer, intent(in) :: base
    integer, intent(out) :: num
 
    integer :: i
 
    ! no PUSH_SUB, called too often
 
    num = m_zero
    do i = 1, size(converted)
      num = num + converted(i) * base**(i-1)
    end do
 
  end subroutine convert_from_base
 
  elemental real(real64) function clamp(x, xmin, xmax)
    real(real64), intent(in) :: x, xmin, xmax
    clamp = min(max(x, xmin), xmax)
  end function clamp
 
#include "undef.F90"
#include "complex.F90"
#include "math_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "math_inc.F90"
 
end module math_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: