qshep.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 qshep_oct_m
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::            &
    qshep_t,           &
    qshep_init,        &
    qshep_interpolate, &
    qshep_end
 
  type qshepr_t
    private
    integer(int64) :: npoints, nq, nw, nr, dim
    integer(int64), allocatable :: lcell(:, :, :), lnext(:)
    real(real64), allocatable :: rsq(:), a(:, :)
    real(real64) :: xmin, ymin, dx, dy, rmax, xyzmin(3), xyzdel(3)
 
    real(real64), allocatable :: x(:), y(:), z(:)
  end type qshepr_t
 
  type qshep_t
    private
    integer(int64) :: kind
    type(qshepr_t) :: re, im
  end type qshep_t
 
  interface qshep_init
    module procedure dqshep_init, zqshep_init
  end interface
 
  interface qshep_interpolate
    module procedure dqshep_interpolate, zqshep_interpolate
  end interface qshep_interpolate
 
  interface
    real(real64) function qs2val(px, py, n, x, y, f, nr, lcell, lnext, xmin, &
      ymin, dx, dy, rmax, rsq, a)
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64) :: px
      real(real64) :: py
      integer(int64) :: n
      real(real64) :: x(n)
      real(real64) :: y(n)
      real(real64) :: f(n)
      integer(int64) :: nr
      integer(int64) :: lcell(nr, nr)
      integer(int64) :: lnext(n)
      real(real64) :: xmin
      real(real64) :: ymin
      real(real64) :: dx
      real(real64) :: dy
      real(real64) :: rmax
      real(real64) :: rsq(n)
      real(real64) :: a(5, n)
    end function qs2val
 
    subroutine qs2grd(px, py, n, x, y, f, nr, lcell, lnext, xmin, &
      ymin, dx, dy, rmax, rsq, a, q, qx, qy, ier)
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64) :: px
      real(real64) :: py
      integer(int64) :: n
      real(real64) :: x(n)
      real(real64) :: y(n)
      real(real64) :: f(n)
      integer(int64) :: nr
      integer(int64) :: lcell(nr, nr)
      integer(int64) :: lnext(n)
      real(real64) :: xmin
      real(real64) :: ymin
      real(real64) :: dx
      real(real64) :: dy
      real(real64) :: rmax
      real(real64) :: rsq(n)
      real(real64) :: a(5, n)
      real(real64) :: q
      real(real64) :: qx
      real(real64) :: qy
      integer(int64) :: ier
    end subroutine qs2grd
 
    real(real64) function qs3val(px, py, pz, n, x, y, z, f, nr, lcell, lnext, xyzmin, &
      xyzdel, rmax, rsq, a)
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64) :: px
      real(real64) :: py
      real(real64) :: pz
      integer(int64) :: n
      real(real64) :: x(n)
      real(real64) :: y(n)
      real(real64) :: z(n)
      real(real64) :: f(n)
      integer(int64) :: nr
      integer(int64) :: lcell(nr, nr, nr)
      integer(int64) :: lnext(n)
      real(real64) :: xyzmin(3)
      real(real64) :: xyzdel(3)
      real(real64) :: rmax
      real(real64) :: rsq(n)
      real(real64) :: a(9, n)
    end function qs3val
 
    subroutine qs3grd(px, py, pz, n, x, y, z, f, nr, lcell, lnext, xyzmin, &
      xyzdel, rmax, rsq, a, q, qx, qy, qz, ier)
      use, intrinsic :: iso_fortran_env
      implicit none
      real(real64) :: px
      real(real64) :: py
      real(real64) :: pz
      integer(int64) :: n
      real(real64) :: x(n)
      real(real64) :: y(n)
      real(real64) :: z(n)
      real(real64) :: f(n)
      integer(int64) :: nr
      integer(int64) :: lcell(nr, nr, nr)
      integer(int64) :: lnext(n)
      real(real64) :: xyzmin(3)
      real(real64) :: xyzdel(3)
      real(real64) :: rmax
      real(real64) :: rsq(n)
      real(real64) :: a(9, n)
      real(real64) :: q
      real(real64) :: qx
      real(real64) :: qy
      real(real64) :: qz
      integer(int64) :: ier
    end subroutine qs3grd
  end interface
 
contains
 
  ! ------------------------------------------------------------------------------
  subroutine dqshep_init(interp, npoints, f, x, y, z)
    type(qshep_t), intent(out) :: interp
    integer(int64),    intent(in)  :: npoints
    real(real64),       intent(in)  :: f(:)
    real(real64)                    :: x(:), y(:)
    real(real64),          optional :: z(:)
 
    push_sub(dqshep_init)
 
    interp%kind = 0
    call qshepr_init(interp%re, npoints, f, x, y, z = z)
 
    pop_sub(dqshep_init)
  end subroutine dqshep_init
 
 
  ! ------------------------------------------------------------------------------
  subroutine zqshep_init(interp, npoints, f, x, y, z)
    type(qshep_t), intent(out) :: interp
    integer(int64),    intent(in)  :: npoints
    complex(real64),    intent(in)  :: f(:)
    real(real64)                    :: x(:), y(:)
    real(real64),       optional    :: z(:)
 
    push_sub(zqshep_init)
 
    interp%kind = 1
    call qshepr_init(interp%re, npoints, real(f),  x, y, z = z)
    call qshepr_init(interp%im, npoints, aimag(f), x, y, z = z)
 
    pop_sub(zqshep_init)
  end subroutine zqshep_init
 
 
  ! ------------------------------------------------------------------------------
  subroutine qshepr_init(interp, npoints, f, x, y, z)
    type(qshepr_t), intent(out) :: interp
    integer(int64),     intent(in) :: npoints
    real(real64),        intent(in) :: f(:)
    real(real64),        target     :: x(:), y(:)
    real(real64),  target, optional :: z(:)
 
    integer(int64) :: ier
 
    push_sub(qshepr_init)
 
    interp%npoints = npoints
 
    if (present(z)) then
      interp%dim = 3
    else
      interp%dim = 2
    end if
 
    select case (interp%dim)
    case (2)
      interp%nq = 13
      interp%nw = 19
      interp%nr = nint(sqrt(interp%npoints/3.0_8))
      safe_allocate(interp%lcell(1:interp%nr, 1:interp%nr, 1))
      safe_allocate(interp%a(1:5, 1:npoints))
    case (3)
      interp%nq = 17 ! This is the recommended value in qshep3d.f90
      interp%nw = 16 ! The recommended value in qshep3d.f90 is 32, but this speeds up things.
      interp%nr = nint((interp%npoints/3.0_8)**(1.0_8/3.0_8))
      safe_allocate(interp%lcell(1:interp%nr, 1:interp%nr, 1:interp%nr))
      safe_allocate(interp%a(1:9, 1:interp%npoints))
    end select
 
    safe_allocate(interp%lnext(1:npoints))
    safe_allocate(interp%rsq(1:npoints))
 
    select case (interp%dim)
    case (2)
      call qshep2(npoints, x, y, f, interp%nq, interp%nw, interp%nr, interp%lcell(:, :, 1), &
        interp%lnext, interp%xmin, interp%ymin, interp%dx, interp%dy, &
        interp%rmax, interp%rsq, interp%a, ier)
      safe_allocate(interp%x(1:npoints))
      safe_allocate(interp%y(1:npoints))
      interp%x(1:npoints) = x(1:npoints)
      interp%y(1:npoints) = y(1:npoints)
    case (3)
      call qshep3(npoints, x, y, z, f, interp%nq, interp%nw, interp%nr, interp%lcell, &
        interp%lnext, interp%xyzmin, interp%xyzdel, interp%rmax, interp%rsq, &
        interp%a, ier)
      safe_allocate(interp%x(1:npoints))
      safe_allocate(interp%y(1:npoints))
      safe_allocate(interp%z(1:npoints))
      interp%x(1:npoints) = x(1:npoints)
      interp%y(1:npoints) = y(1:npoints)
      interp%z(1:npoints) = z(1:npoints)
    end select
 
    pop_sub(qshepr_init)
  end subroutine qshepr_init
 
 
  ! ------------------------------------------------------------------------------
  real(real64) function qshep_interpolater(interp, f, p, gf) result(v)
    type(qshepr_t),    intent(in)    :: interp
    real(real64),           intent(in)    :: f(:)
    real(real64),           intent(in)    :: p(:)
    real(real64), optional, intent(inout) :: gf(:)
 
    integer(int64) :: ier
 
    push_sub(qshep_interpolater)
 
    select case (interp%dim)
    case (2)
      if (present(gf)) then
        call qs2grd( p(1), p(2), interp%npoints, interp%x, interp%y, &
          f, interp%nr, interp%lcell(:, :, 1), interp%lnext, interp%xmin, &
          interp%ymin, interp%dx, interp%dy, interp%rmax, interp%rsq, interp%a, &
          v, gf(1), gf(2), ier)
      else
        v = qs2val( p(1), p(2), interp%npoints, interp%x, interp%y, &
          f, interp%nr, interp%lcell(:, :, 1), interp%lnext, interp%xmin, &
          interp%ymin, interp%dx, interp%dy, interp%rmax, interp%rsq, interp%a)
      end if
    case (3)
      if (present(gf)) then
        call qs3grd( p(1), p(2), p(3), interp%npoints, interp%x, interp%y, interp%z, &
          f, interp%nr, interp%lcell, interp%lnext, interp%xyzmin, &
          interp%xyzdel, interp%rmax, interp%rsq, interp%a , &
          v, gf(1), gf(2), gf(3), ier)
      else
        v = qs3val(p(1), p(2), p(3), interp%npoints, interp%x, interp%y, interp%z, &
          f, interp%nr, interp%lcell, interp%lnext, interp%xyzmin, &
          interp%xyzdel, interp%rmax, interp%rsq, interp%a)
      end if
    end select
 
    pop_sub(qshep_interpolater)
  end function qshep_interpolater
 
 
  ! ------------------------------------------------------------------------------
  real(real64) function dqshep_interpolate(interp, f, p, gf) result(v)
    type(qshep_t),     intent(in)    :: interp
    real(real64),           intent(in)    :: f(:)
    real(real64),           intent(in)    :: p(:)
    real(real64), optional, intent(inout) :: gf(:)
 
    push_sub(dqshep_interpolate)
 
    v = qshep_interpolater(interp%re, f, p, gf = gf)
 
    pop_sub(dqshep_interpolate)
  end function dqshep_interpolate
 
 
  ! ------------------------------------------------------------------------------
  complex(real64) function zqshep_interpolate(interp, f, p, gf) result(v)
    type(qshep_t),        intent(in)    :: interp
    complex(real64),           intent(in)    :: f(:)
    real(real64),              intent(in)    :: p(:)
    complex(real64), optional, intent(inout) :: gf(:)
 
    integer(int64) :: i
    real(real64), allocatable :: rgf(:), igf(:)
 
    push_sub(zqshep_interpolate)
 
    if (present(gf)) then
      safe_allocate(rgf(1:size(gf)))
      safe_allocate(igf(1:size(gf)))
      v = cmplx(qshep_interpolater(interp%re, real(f), p, rgf), &
        qshep_interpolater(interp%im, aimag(f), p, igf), 8)
      do i = 1, size(gf)
        gf(i) = cmplx( rgf(i), igf(i), 8)
      end do
      safe_deallocate_a(rgf)
      safe_deallocate_a(igf)
    else
      v = cmplx(qshep_interpolater(interp%re, real(f), p), &
        qshep_interpolater(interp%im, aimag(f), p), 8)
    end if
 
    pop_sub(zqshep_interpolate)
  end function zqshep_interpolate
 
 
  ! ------------------------------------------------------------------------------
  subroutine qshep_end(interp)
    type(qshep_t), intent(inout) :: interp
 
    push_sub(qshep_end)
 
    call qshepr_end(interp%re)
    if (interp%kind == 1) call qshepr_end(interp%im)
 
    pop_sub(qshep_end)
  end subroutine qshep_end
 
 
  ! ------------------------------------------------------------------------------
  subroutine qshepr_end(interp)
    type(qshepr_t), intent(inout) :: interp
 
    push_sub(qshepr_end)
 
    safe_deallocate_a(interp%lcell)
    safe_deallocate_a(interp%lnext)
    safe_deallocate_a(interp%rsq)
    safe_deallocate_a(interp%a)
    safe_deallocate_a(interp%x)
    safe_deallocate_a(interp%y)
    safe_deallocate_a(interp%z)
 
    pop_sub(qshepr_end)
  end subroutine qshepr_end
 
end module qshep_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: