speclib.f90

Go to the documentation of this file.

! Copyright (c) 2008-2020, UCHICAGO ARGONNE, LLC.
!
! The UChicago Argonne, LLC as Operator of Argonne National
! Laboratory holds copyright in the Software. The copyright holder
! reserves all rights except those expressly granted to licensees,
! and U.S. Government license rights.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! 1. Redistributions of source code must retain the above copyright
! notice, this list of conditions and the disclaimer below.
!
! 2. Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the disclaimer (as noted below)
! in the documentation and/or other materials provided with the
! distribution.
!
! 3. Neither the name of ANL nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
! UCHICAGO ARGONNE, LLC, THE U.S. DEPARTMENT OF
! ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
! TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Additional BSD Notice
! ---------------------
! 1. This notice is required to be provided under our contract with
! the U.S. Department of Energy (DOE). This work was produced at
! Argonne National Laboratory under Contract
! No. DE-AC02-06CH11357 with the DOE.
!
! 2. Neither the United States Government nor UCHICAGO ARGONNE,
! LLC nor any of their employees, makes any warranty,
! express or implied, or assumes any liability or responsibility for the
! accuracy, completeness, or usefulness of any information, apparatus,
! product, or process disclosed, or represents that its use would not
! infringe privately-owned rights.
!
! 3. Also, reference herein to any specific commercial products, process,
! or services by trade name, trademark, manufacturer or otherwise does
! not necessarily constitute or imply its endorsement, recommendation,
! or favoring by the United States Government or UCHICAGO ARGONNE LLC.
! The views and opinions of authors expressed
! herein do not necessarily state or reflect those of the United States
! Government or UCHICAGO ARGONNE, LLC, and shall
! not be used for advertising or product endorsement purposes.
!
!==============================================================================
!
!
!------------------------------------------------------------------------------
!
!     ABBRIVIATIONS:
!
!     M   - Set of mesh points
!     Z   - Set of collocation/quadrature points
!     W   - Set of quadrature weights
!     H   - Lagrangian interpolant
!     D   - Derivative operator
!     I   - Interpolation operator
!     GL  - Gauss Legendre
!     GLL - Gauss-Lobatto Legendre
!     GJ  - Gauss Jacobi
!     GLJ - Gauss-Lobatto Jacobi
!
!
!     MAIN ROUTINES:
!
!     Points and weights:
!
!     ZWGL      Compute Gauss Legendre points and weights
!     ZWGLL     Compute Gauss-Lobatto Legendre points and weights
!     ZWGJ      Compute Gauss Jacobi points and weights (general)
!     ZWGLJ     Compute Gauss-Lobatto Jacobi points and weights (general)
!
!     Lagrangian interpolants:
!
!     HGL       Compute Gauss Legendre Lagrangian interpolant
!     HGLL      Compute Gauss-Lobatto Legendre Lagrangian interpolant
!     HGJ       Compute Gauss Jacobi Lagrangian interpolant (general)
!     HGLJ      Compute Gauss-Lobatto Jacobi Lagrangian interpolant (general)
!
!     Derivative operators:
!
!     DGLL      Compute Gauss-Lobatto Legendre derivative matrix
!     DGLLGL    Compute derivative matrix for a staggered mesh (GLL->GL)
!     DGJ       Compute Gauss Jacobi derivative matrix (general)
!     DGLJ      Compute Gauss-Lobatto Jacobi derivative matrix (general)
!     DGLJGJ    Compute derivative matrix for a staggered mesh (GLJ->GJ) (general)
!
!     Interpolation operators:
!
!     IGLM      Compute interpolation operator GL  -> M
!     IGLLM     Compute interpolation operator GLL -> M
!     IGJM      Compute interpolation operator GJ  -> M  (general)
!     IGLJM     Compute interpolation operator GLJ -> M  (general)
!
!     Other:
!
!     PNLEG     Compute Legendre polynomial of degree N
!     legendre_poly Compute Legendre polynomial of degree 0-N
!     PNDLEG    Compute derivative of Legendre polynomial of degree N
!
!     Comments:
!
!     Note that many of the above routines exist in both single and
!     double precision. If the name of the single precision routine is
!     SUB, the double precision version is called SUBD. In most cases
!     all the "low-level" arithmetic is done in double precision, even
!     for the single precsion versions.
!
!     Useful references:
!
! [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
!     Providence, Rhode Island, 1939.
! [2] Abramowitz & Stegun: Handbook of Mathematical Functions,
!     Dover, New York, 1972.
! [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid
!     Dynamics, Springer-Verlag, 1988.
!
!
!==============================================================================
module speclib
  use num_types, only : rp, xp
  use math, only: abscmp
  use utils, only : neko_error
  use, intrinsic :: iso_fortran_env, only : stderr => error_unit
  implicit none
 
contains
 
  subroutine zwgl(Z, W, NP)
    integer, intent(in) :: NP
    real(kind=rp), intent(inout) :: z(np), w(np)
    real(kind=rp) alpha, beta
    alpha = 0.0_rp
    beta = 0.0_rp
    call zwgj(z, w, np, alpha, beta)
  end subroutine zwgl
 
 
  subroutine zwgll(Z, W, NP)
    integer, intent(in) :: NP
    real(kind=rp), intent(inout) :: z(np), w(np)
    real(kind=rp) alpha, beta
    alpha = 0.0_rp
    beta = 0.0_rp
    call zwglj(z, w, np, alpha, beta)
  end subroutine zwgll
 
  subroutine zwgj(Z, W, NP, ALPHA, BETA)
    integer, intent(in) :: NP
    real(kind=rp), intent(inout) :: z(np), w(np)
    real(kind=rp), intent(in) :: alpha, beta
 
    integer, parameter :: NMAX = 84
    integer, parameter :: NZD = nmax
 
    real(kind=xp) zd(nzd), wd(nzd), alphad, betad
    integer :: I, NPMAX
 
    npmax = nzd
    if (np .gt. npmax) then
       write (stderr, *) 'Too large polynomial degree in ZWGJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here NP=', np
       call neko_error
    end if
 
    alphad = real(alpha, kind=xp)
    betad = real(beta, kind=xp)
    call zwgjd(zd, wd, np, alphad, betad)
    do i = 1, np
       z(i) = real(zd(i), kind=rp)
       w(i) = real(wd(i), kind=rp)
    end do
  end subroutine zwgj
 
  subroutine zwgjd(Z, W, NP, ALPHA, BETA)
    integer, intent(in) :: NP
    real(kind=xp), intent(inout) :: z(np), w(np)
    real(kind=xp), intent(in) :: alpha, beta
 
    real(kind=xp) :: dn, apb
    real(kind=xp) :: fac1, fac2, fac3, fnorm
    real(kind=xp) :: rcoef, p, pd, pm1, pdm1, pm2, pdm2
    real(kind=xp) :: dnp1, dnp2
    integer :: N, NP1, NP2, I
 
    n = np - 1
    dn = real(n, kind=xp)
 
    apb = alpha + beta
 
    if (np .le. 0) then
       write (stderr, *) 'ZWGJD: Minimum number of Gauss points is 1', np
       call neko_error
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('ZWGJD: Alpha and Beta must be greater than -1')
    end if
 
    if (np .eq. 1) then
       z(1) = (beta - alpha) / (apb + 2.0_xp)
       w(1) = gammaf(alpha + 1.0_xp) * gammaf(beta + 1.0_xp) / &
            gammaf(apb + 2.0_xp) * 2.0_xp**(apb + 1.0_xp)
       return
    end if
 
    call jacg(z, np, alpha, beta)
 
    np1 = n + 1
    np2 = n + 2
    dnp1 = real(np1, kind=xp)
    dnp2 = real(np2, kind=xp)
    fac1 = dnp1 + alpha + beta + 1.0_xp
    fac2 = fac1 + dnp1
    fac3 = fac2 + 1.0_xp
    fnorm = pnormj(np1, alpha, beta)
    rcoef = (fnorm*fac2*fac3) / (2.0_xp*fac1*dnp2)
    do i = 1, np
       call jacobf(p, pd, pm1, pdm1, pm2, pdm2, np2, alpha, beta, z(i))
       w(i) = -rcoef/(p*pdm1)
    end do
  end subroutine zwgjd
 
  subroutine zwglj(Z, W, NP, ALPHA, BETA)
    integer, intent(in) :: NP
    real(kind=rp), intent(inout) :: z(np), w(np)
    real(kind=rp), intent(in) :: alpha, beta
 
    integer, parameter :: NMAX = 84
    integer, parameter :: NZD = nmax
 
    real(kind=xp) zd(nzd), wd(nzd), alphad, betad
    integer :: I, NPMAX
 
    npmax = nzd
    if (np .gt. npmax) then
       write (stderr, *) 'Too large polynomial degree in ZWGLJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here NP=', np
       call neko_error
    end if
    alphad = real(alpha, kind=xp)
    betad = real(beta, kind=xp)
    call zwgljd(zd, wd, np, alphad, betad)
    do i = 1, np
       z(i) = real(zd(i), kind=rp)
       w(i) = real(wd(i), kind=rp)
    end do
  end subroutine zwglj
 
 
  subroutine zwgljd(Z, W, NP, ALPHA, BETA)
 
    integer, intent(in) :: NP
    real(kind=xp), intent(inout) :: z(np), w(np)
    real(kind=xp), intent(in) :: alpha, beta
 
    real(kind=xp) :: alpg, betg
    real(kind=xp) :: p, pd, pm1, pdm1, pm2, pdm2
    integer :: N, NM1, I
 
    n = np - 1
    nm1 = n - 1
 
    if (np .le. 1) then
       write (stderr, *) 'ZWGLJD: Minimum number of Gauss-Lobatto points is 2'
       write (stderr, *) 'ZWGLJD: alpha, beta:', alpha, beta, np
       call neko_error
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('ZWGLJD: Alpha and Beta must be greater than -1')
    end if
 
    if (nm1 .gt. 0) then
       alpg = alpha + 1.0_xp
       betg = beta + 1.0_xp
       call zwgjd(z(2), w(2), nm1, alpg, betg)
    end if
 
    z(1) = -1.0_xp
    z(np) = 1.0_xp
    do i = 2, np - 1
       w(i) = w(i) / (1.0_xp-z(i)**2)
    end do
    call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(1))
    w(1) = endw1(n, alpha, beta) / (2.0_xp*pd)
    call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(np))
    w(np) = endw2(n, alpha, beta) / (2.0_xp*pd)
 
  end subroutine zwgljd
 
  real(kind=xp) function endw1(N, ALPHA, BETA)
 
    real(kind=xp), intent(in) :: alpha, beta
    integer, intent(in) :: n
 
    real(kind=xp) :: apb
    real(kind=xp) :: f1, f2, f3, fint1, fint2
    real(kind=xp) :: a1, a2, a3, di, abn, abnn
    integer :: i
 
    if (n .eq. 0) then
       endw1 = 0.0_xp
       return
    end if
 
    apb = alpha + beta
    f1 = gammaf(alpha + 2.0_xp)*gammaf(beta + 1.0_xp) / gammaf(apb + 3.0_xp)
    f1 = f1*(apb + 2.0_xp)*2.0_xp**(apb + 2.0_xp)/2.0_xp
    if (n .eq. 1) then
       endw1 = f1
       return
    end if
 
    fint1 = gammaf(alpha + 2.0_xp)*gammaf(beta + 1.0_xp) / gammaf(apb + 3.0_xp)
    fint1 = fint1*2.0_xp**(apb + 2.0_xp)
    fint2 = gammaf(alpha + 2.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 4.0_xp)
    fint2 = fint2*2.0_xp**(apb + 3.0_xp)
    f2 = (-2.0_xp*(beta + 2.0_xp)*fint1 + (apb + 4.0_xp)*fint2) * &
         (apb + 3.0_xp) / 4.0_xp
    if (n .eq. 2) then
       endw1 = f2
       return
    end if
 
    do i = 3, n
       di = real(i - 1, kind=xp)
       abn = alpha + beta + di
       abnn = abn + di
       a1 = -(2.0_xp*(di + alpha) * (di + beta)) / (abn*abnn*(abnn + 1.0_xp))
       a2 = (2.0_xp*(alpha - beta)) / (abnn*(abnn + 2.0_xp))
       a3 = (2.0_xp*(abn + 1.0_xp)) / ((abnn + 2.0_xp) * (abnn + 1.0_xp))
       f3 = -(a2*f2 + a1*f1) / a3
       f1 = f2
       f2 = f3
    end do
    endw1 = f3
  end function endw1
 
  real(kind=xp) function endw2(N, ALPHA, BETA)
 
    real(kind=xp), intent(in) :: alpha, beta
    integer, intent(in) :: n
 
    real(kind=xp) :: apb
    real(kind=xp) :: f1, f2, f3, fint1, fint2
    real(kind=xp) :: a1, a2, a3, di, abn, abnn
    integer :: i
 
 
    if (n .eq. 0) then
       endw2 = 0.0_xp
       return
    end if
 
    apb = alpha + beta
    f1 = gammaf(alpha + 1.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 3.0_xp)
    f1 = f1*(apb + 2.0_xp)*2.0_xp**(apb + 2.0_xp)/2.0_xp
    if (n .eq. 1) then
       endw2 = f1
       return
    end if
 
    fint1 = gammaf(alpha + 1.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 3.0_xp)
    fint1 = fint1*2.0_xp**(apb + 2.0_xp)
    fint2 = gammaf(alpha + 2.0_xp)*gammaf(beta + 2.0_xp) / gammaf(apb + 4.0_xp)
    fint2 = fint2*2.0_xp**(apb + 3.0_xp)
    f2 = (2.0_xp*(alpha + 2.0_xp)*fint1 - (apb + 4.0_xp)*fint2) * &
         (apb + 3.0_xp) / 4.0_xp
    if (n .eq. 2) then
       endw2 = f2
       return
    end if
 
    do i = 3, n
       di = ((i-1))
       abn = alpha + beta + di
       abnn = abn + di
       a1 = -(2.0_xp*(di + alpha) * (di + beta)) / (abn*abnn*(abnn + 1.0_xp))
       a2 = (2.0_xp*(alpha - beta)) / (abnn*(abnn + 2.0_xp))
       a3 = (2.0_xp*(abn + 1.0_xp)) / ((abnn + 2.0_xp) * (abnn + 1.0_xp))
       f3 = -(a2*f2 + a1*f1)/a3
       f1 = f2
       f2 = f3
    end do
    endw2 = f3
  end function endw2
 
  real(kind=xp) function gammaf(X)
    real(kind=xp), intent(in) :: x
    real(kind=xp), parameter :: pi = 4.0_xp*atan(1.0_xp)
 
    gammaf = 1.0_xp
    if (abscmp(x, -0.5_xp)) gammaf = -2.0_xp*sqrt(pi)
    if (abscmp(x, 0.5_xp)) gammaf = sqrt(pi)
    if (abscmp(x, 1.0_xp)) gammaf = 1.0_xp
    if (abscmp(x, 2.0_xp)) gammaf = 1.0_xp
    if (abscmp(x, 1.5_xp)) gammaf = sqrt(pi) / 2.0_xp
    if (abscmp(x, 2.5_xp)) gammaf = 1.5_xp * sqrt(pi) / 2.0_xp
    if (abscmp(x, 3.5_xp)) gammaf = 0.5_xp * (2.5_xp * (1.5_xp * sqrt(pi)))
    if (abscmp(x, 3.0_xp)) gammaf = 2.0_xp
    if (abscmp(x, 4.0_xp)) gammaf = 6.0_xp
    if (abscmp(x, 5.0_xp)) gammaf = 24.0_xp
    if (abscmp(x, 6.0_xp)) gammaf = 120.0_xp
  end function gammaf
 
  real(kind=xp) function pnormj(N, ALPHA, BETA)
    real(kind=xp), intent(in) :: alpha, beta
    integer, intent(in) :: n
 
    real(kind=xp) :: dn, dindx
    real(kind=xp) :: const, prod, frac
    integer :: i
 
    dn = real(n, kind=xp)
    const = alpha + beta + 1.0_xp
    if (n .le. 1) then
       prod = gammaf(dn + alpha)*gammaf(dn + beta)
       prod = prod / (gammaf(dn)*gammaf(dn + alpha + beta))
       pnormj = prod * 2.0_xp**const / (2.0_xp*dn + const)
       return
    end if
 
    prod = gammaf(alpha + 1.0_xp)*gammaf(beta + 1.0_xp)
    prod = prod/(2.0_xp*(1.0_xp + const)*gammaf(const + 1.0_xp))
    prod = prod*(1.0_xp + alpha) * (2.0_xp + alpha)
    prod = prod*(1.0_xp + beta) * (2.0_xp + beta)
    do i = 3, n
       dindx = real(i, kind=xp)
       frac = (dindx + alpha) * (dindx + beta) / (dindx*(dindx + alpha + beta))
       prod = prod*frac
    end do
    pnormj = prod*2.0_xp**const / (2.0_xp*dn + const)
  end function pnormj
 
  subroutine jacg(XJAC, NP, ALPHA, BETA)
    integer, intent(in) :: NP
    real(kind=xp), intent(inout) :: xjac(np)
    real(kind=xp), intent(in) :: alpha, beta
 
    integer, parameter :: KSTOP = 10
    real(kind=rp), parameter :: eps = 1.0e-12_rp
    real(kind=xp), parameter :: pi = 4.0_xp*atan(1.0_xp)
 
    real(kind=xp) :: dth, x, x1, x2, xlast, delx, xmin
    real(kind=xp) :: p, pd, pm1, pdm1, pm2, pdm2
    real(kind=xp) :: recsum, swap
    integer :: I, J, K, N, JM, JMIN
 
    n = np - 1
    dth = pi / (2.0_xp*real(n, kind=xp) + 2.0_xp)
    do j = 1, np
       if (j .eq. 1) then
          x = cos((2.0_xp*(real(j, kind=xp) - 1.0_xp) + 1.0_xp)*dth)
       else
          x1 = cos((2.0_xp*(real(j, kind=xp) - 1.0_xp) + 1.0_xp)*dth)
          x2 = xlast
          x = (x1 + x2) / 2.0_xp
       end if
 
       do k = 1, kstop
          call jacobf(p, pd, pm1, pdm1, pm2, pdm2, np, alpha, beta, x)
          recsum = 0.0_xp
          jm = j - 1
          do i = 1, jm
             recsum = recsum + 1.0_xp / (x-xjac(np - i + 1))
          end do
          delx = -p / (pd - recsum*p)
          x = x + delx
          if (abs(delx) .lt. eps) exit
       end do
 
       xjac(np-j + 1) = x
       xlast = x
    end do
 
    do i = 1, np
       xmin = 2.
       do j = i, np
          if (xjac(j) .lt. xmin) then
             xmin = xjac(j)
             jmin = j
          end if
       end do
       if (jmin .ne. i) then
          swap = xjac(i)
          xjac(i) = xjac(jmin)
          xjac(jmin) = swap
       end if
    end do
  end subroutine jacg
 
  subroutine jacobf(POLY, PDER, POLYM1, PDERM1, POLYM2, PDERM2, N, ALP, BET, X)
 
    real(kind=xp), intent(inout) :: poly, pder, polym1, pderm1, polym2, pderm2
    real(kind=xp), intent(in) :: alp, bet, x
    integer, intent(in) :: N
 
    real(kind=xp) :: apb, polyl, pderl, polyn, pdern
    real(kind=xp) :: psave, pdsave
    real(kind=xp) :: a1, a2, a3, a4, b3
    real(kind=xp) :: dk
    integer :: K
 
    apb = alp + bet
    poly = 1.0_xp
    pder = 0.0_xp
    if (n .eq. 0) return
 
    polyl = poly
    pderl = pder
    poly = (alp - bet + (apb + 2.0_xp)*x) / 2.0_xp
    pder = (apb + 2.0_xp) / 2.0_xp
    if (n .eq. 1) return
 
    do k = 2, n
       dk = real(k, kind=xp)
       a1 = 2.0_xp*dk*(dk + apb) * (2.0_xp*dk + apb - 2.0_xp)
       a2 = (2.0_xp*dk + apb - 1.0_xp) * (alp**2 - bet**2)
       b3 = (2.0_xp*dk + apb - 2.0_xp)
       a3 = b3*(b3 + 1.0_xp) * (b3 + 2.0_xp)
       a4 = 2.0_xp*(dk + alp - 1.0_xp) * (dk + bet - 1.0_xp) * (2.0_xp*dk + apb)
       polyn = ((a2 + a3*x)*poly - a4*polyl) / a1
       pdern = ((a2 + a3*x)*pder - a4*pderl + a3*poly) / a1
       psave = polyl
       pdsave = pderl
       polyl = poly
       poly = polyn
       pderl = pder
       pder = pdern
    end do
    polym1 = polyl
    pderm1 = pderl
    polym2 = psave
    pderm2 = pdsave
  end subroutine jacobf
 
  real(kind=xp) function hgj(II, Z, ZGJ, NP, ALPHA, BETA)
    integer, intent(in) :: np, ii
    real(kind=xp), intent(in) :: z, zgj(np), alpha, beta
 
    integer, parameter :: nmax = 84
    integer, parameter :: nzd = nmax
 
    real(kind=xp) zd, zgjd(nzd)
    integer :: i, npmax
 
    npmax = nzd
    if (np .gt. npmax) then
       write (stderr, *) 'Too large polynomial degree in HGJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here NP=', np
       call neko_error
    end if
 
    zd = z
    do i = 1, np
       zgjd(i) = zgj(i)
    end do
    hgj = hgjd(ii, zd, zgjd, np, alpha, beta)
  end function hgj
 
  real(kind=xp) function hgjd(II, Z, ZGJ, NP, ALPHA, BETA)
    integer, intent(in) :: np, ii
    real(kind=xp), intent(in) :: z, zgj(np), alpha, beta
 
    real(kind=xp) :: eps, zi, dz
    real(kind=xp) :: pz, pdz, pzi, pdzi, pm1, pdm1, pm2, pdm2
 
    eps = 1.0e-5_xp
    zi = zgj(ii)
    dz = z - zi
    if (abs(dz) .lt. eps) then
       hgjd = 1.0_xp
       return
    end if
    call jacobf(pzi, pdzi, pm1, pdm1, pm2, pdm2, np, alpha, beta, zi)
    call jacobf(pz, pdz, pm1, pdm1, pm2, pdm2, np, alpha, beta, z)
    hgjd = pz / (pdzi*(z-zi))
  end function hgjd
 
  real(kind=xp) function hglj(II, Z, ZGLJ, NP, ALPHA, BETA)
    integer, intent(in) :: np, ii
    real(kind=xp), intent(in) :: z, zglj(np), alpha, beta
 
    integer, parameter :: nmax = 84
    integer, parameter :: nzd = nmax
 
    real(kind=xp) zd, zgljd(nzd)
    integer :: i, npmax
 
    npmax = nzd
    if (np .gt. npmax) then
       write (stderr, *) 'Too large polynomial degree in HGLJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here NP=', np
       call neko_error
    end if
    zd = z
    do i = 1, np
       zgljd(i) = zglj(i)
    end do
    hglj = hgljd(ii, zd, zgljd, np, alpha, beta)
  end function hglj
 
  real(kind=xp) function hgljd(I, Z, ZGLJ, NP, ALPHA, BETA)
    integer, intent(in) :: np, i
    real(kind=xp), intent(in) :: z, zglj(np), alpha, beta
 
    real(kind=xp) :: eps, zi, dz, dn
    real(kind=xp) :: p, pd, pi, pdi, pm1, pdm1, pm2, pdm2
    real(kind=xp) :: eigval, const
    integer :: n
 
    eps = 1.0e-5_xp
    zi = zglj(i)
    dz = z-zi
    if (abs(dz) .lt. eps) then
       hgljd = 1.0_xp
       return
    end if
 
    n = np - 1
    dn = real(n, kind=xp)
    eigval = -dn*(dn + alpha + beta + 1.0_xp)
    call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, n, alpha, beta, zi)
    const = eigval*pi + alpha*(1.0_xp + zi)*pdi - beta*(1.0_xp - zi)*pdi
    call jacobf(p, pd, pm1, pdm1, pm2, pdm2, n, alpha, beta, z)
    hgljd = (1.0_xp - z**2)*pd / (const*(z - zi))
  end function hgljd
 
  subroutine dgj(D, DT, Z, NZ, NZD, ALPHA, BETA)
    integer, intent(in) :: NZ, NZD
    real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
    real(kind=xp), intent(in) :: z(nz), alpha, beta
 
    integer, parameter :: NMAX = 84
    integer, parameter :: NZDD = nmax
 
    real(kind=xp) :: dd(nzdd, nzdd), dtd(nzdd, nzdd), zd(nzdd)
    integer :: I, J
 
    if (nz .le. 0) then
       call neko_error('DGJ: Minimum number of Gauss points is 1')
    else if (nz .gt. nmax) then
       write (stderr, *) 'Too large polynomial degree in DGJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here Nz=', nz
       call neko_error
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGJ: Alpha and Beta must be greater than -1')
    end if
 
    do i = 1, nz
       zd(i) = z(i)
    end do
    call dgjd(dd, dtd, zd, nz, nzdd, alpha, beta)
    do i = 1, nz
       do j = 1, nz
          d(i, j) = dd(i, j)
          dt(i, j) = dtd(i, j)
       end do
    end do
  end subroutine dgj
 
  subroutine dgjd(D, DT, Z, NZ, NZD, ALPHA, BETA)
    integer, intent(in) :: NZ, NZD
    real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
    real(kind=xp), intent(in) :: z(nz), alpha, beta
 
    real(kind=xp) :: dn
    real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
    integer :: I, J, N
 
    n = nz - 1
    dn = real(n, kind=xp)
 
 
    if (nz .le. 1) then
       call neko_error('DGJD: Minimum number of Gauss-Lobatto points is 2')
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGJD: Alpha and Beta must be greater than -1')
    end if
 
    do i = 1, nz
       do j = 1, nz
          call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, nz, alpha, beta, z(i))
          call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, nz, alpha, beta, z(j))
          if (i .ne. j) then
             d(i, j) = pdi / (pdj*(z(i) - z(j)))
          else
             d(i, j) = ((alpha + beta + 2.0_xp)*z(i) + alpha - beta) / &
                  (2.0_xp*(1.0_xp - z(i)**2))
          end if
          dt(j, i) = d(i, j)
       end do
    end do
  end subroutine dgjd
 
  subroutine dglj(D, DT, Z, NZ, NZD, ALPHA, BETA)
    integer, parameter :: NMAX = 84
    integer, parameter :: NZDD = nmax
    integer, intent(in) :: NZ, NZD
    real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
    real(kind=xp), intent(in) :: z(nz), alpha, beta
 
    real(kind=xp) :: dd(nzdd, nzdd), dtd(nzdd, nzdd), zd(nzdd)
    integer :: I, J
 
    if (nz .le. 1) then
       call neko_error('DGLJ: Minimum number of Gauss-Lobatto points is 2')
    else if (nz .gt. nmax) then
       write (stderr, *) 'Too large polynomial degree in DGLJ'
       write (stderr, *) 'Maximum polynomial degree is', nmax
       write (stderr, *) 'Here NZ=', nz
       call neko_error
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGLJ: Alpha and Beta must be greater than -1')
    end if
 
    do i = 1, nz
       zd(i) = z(i)
    end do
    call dgljd(dd, dtd, zd, nz, nzdd, alpha, beta)
    do i = 1, nz
       do j = 1, nz
          d(i, j) = dd(i, j)
          dt(i, j) = dtd(i, j)
       end do
    end do
  end subroutine dglj
 
 
  subroutine dgljd(D, DT, Z, NZ, NZD, ALPHA, BETA)
    integer, intent(in) :: NZ, NZD
    real(kind=xp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
    real(kind=xp), intent(in) :: z(nz), alpha, beta
 
    real(kind=xp) :: dn, eigval
    real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
    real(kind=xp) :: ci, cj
    integer :: I, J, N
 
    n = nz - 1
    dn = real(n, kind=xp)
 
    eigval = -dn*(dn + alpha + beta + 1.0_xp)
 
    if (nz .le. 1) then
       call neko_error('DGLJD: Minimum number of Gauss-Lobatto points is 2')
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGLJD: Alpha and Beta must be greater than -1')
    end if
 
    do i = 1, nz
       do j = 1, nz
          call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(i))
          call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, n, alpha, beta, z(j))
          ci = eigval*pi - (beta*(1.0_xp - z(i)) - alpha*(1.0_xp + z(i)))*pdi
          cj = eigval*pj - (beta*(1.0_xp - z(j)) - alpha*(1.0_xp + z(j)))*pdj
 
          ! Todo: This should have some elses in there
          if (i .ne. j) then
             d(i, j) = ci / (cj*(z(i) - z(j)))
          else if (i .eq. 1) then
             d(i, j) = (eigval + alpha) / (2.0_xp*(beta + 2.0_xp))
          else if (i .eq. nz) then
             d(i, j) = -(eigval + beta) / (2.0_xp*(alpha + 2.0_xp))
          else
             d(i, j) = (alpha*(1.0_xp + z(i)) - beta*(1.0_xp - z(i))) / &
                  (2.0_xp*(1.0_xp - z(i)**2))
          end if
          dt(j, i) = d(i, j)
       end do
    end do
  end subroutine dgljd
 
  subroutine dgll(D, DT, Z, NZ, NZD)
 
    integer, intent(in) :: NZ, NZD
    real(kind=rp), intent(inout) :: d(nzd, nzd), dt(nzd, nzd)
    real(kind=rp), intent(in) :: z(nz)
 
    integer, parameter :: NMAX = 84
 
    real(kind=xp) :: d0, fn
    integer :: I, J, N
 
    n = nz - 1
    if (nz .gt. nmax) then
       write (stderr, *) 'Subroutine DGLL'
       write (stderr, *) 'Maximum polynomial degree =', nmax
       write (stderr, *) 'Polynomial degree         =', nz
       call neko_error
    else if (nz .eq. 1) then
       d(1, 1) = 0.0_rp
       return
    end if
 
    fn = real(n, kind=xp)
    d0 = fn*(fn + 1.0_xp)/4.0_xp
    do i = 1, nz
       do j = 1, nz
          if (i .ne. j) then
             d(i, j) = pnleg(real(z(i), xp), n)/ &
                  (pnleg(real(z(j), xp), n) * (z(i) - z(j)))
          else if (i .eq. 1) then
             d(i, j) = -d0
          else if (i .eq. nz) then
             d(i, j) = d0
          else
             d(i, j) = 0.0_rp
          end if
          dt(j, i) = d(i, j)
       end do
    end do
  end subroutine dgll
 
  real(kind=xp) function hgll(I, Z, ZGLL, NZ)
    integer, intent(in) :: i, nz
    real(kind=xp), intent(in) :: zgll(nz), z
 
    real(kind=xp) :: eps, dz
    real(kind=xp) :: alfan
    integer :: n
 
    eps = 1.0e-5_xp
    dz = z - zgll(i)
    if (abs(dz) .lt. eps) then
       hgll = 1.0_xp
       return
    end if
 
    n = nz - 1
    alfan = real(n, kind=xp) * (real(n, kind=xp) + 1.0_xp)
    hgll = -(1.0_xp - z*z)*pndleg(z, n) / (alfan*pnleg(zgll(i), n) * &
         (z - zgll(i)))
  end function hgll
 
  real(kind=xp) function hgl (I, Z, ZGL, NZ)
    integer, intent(in) :: i, nz
    real(kind=xp), intent(in) :: zgl(nz), z
    real(kind=xp) :: eps, dz
 
    integer :: n
 
    eps = 1.0e-5_xp
    dz = z - zgl(i)
    if (abs(dz) .lt. eps) then
       hgl = 1.0_xp
       return
    end if
 
    n = nz - 1
    hgl = pnleg(z, nz) / (pndleg(zgl(i), nz) * (z - zgl(i)))
  end function hgl
 
  real(kind=xp) function pnleg(Z, N)
 
!---------------------------------------------------------------------
!
!     This next statement is to overcome the underflow bug in the i860.
!     It can be removed at a later date.  11 Aug 1990   pff.
!
 
    real(kind=xp), intent(in) :: z
    integer, intent(in) :: n
 
    real(kind=xp) :: p1, p2, p3, fk
    integer :: k
 
    p1 = 1.0_xp
    if (n .eq. 0) then
       pnleg = p1
       return
    end if
 
    p2 = z
    p3 = p2
    do k = 1, n-1
       fk = real(k, kind=xp)
       p3 = ((2.0_xp*fk + 1.0_xp)*z*p2 - fk*p1) / (fk + 1.0_xp)
       p1 = p2
       p2 = p3
    end do
    pnleg = p3
  end function pnleg
 
  subroutine legendre_poly(L, x, N)
    integer, intent(in) :: N
    real(kind=rp), intent(inout) :: l(0:n)
    real(kind=rp), intent(in) :: x
 
    real(kind=rp) :: dj
    integer :: j
 
    l(0) = 1.0_rp
    if (n .eq. 0) return
    l(1) = x
 
    do j = 1, n-1
       dj = real(j, kind=rp)
       l(j + 1) = ((2.0_rp*dj + 1.0_rp)*x*l(j) - dj*l(j-1)) / (dj + 1.0_rp)
    end do
  end subroutine legendre_poly
 
  real(kind=xp) function pndleg(Z, N)
    real(kind=xp), intent(in) :: z
    integer, intent(in) :: n
 
    real(kind=xp) :: p1, p2, p3, p1d, p2d, p3d, fk
    integer :: k
 
    if (n .eq. 0) then
       pndleg = 0.0_xp
       return
    end if
 
    p1 = 1.0_xp
    p2 = z
    p1d = 0.0_xp
    p2d = 1.0_xp
    p3d = 1.0_xp
    do k = 1, n-1
       fk = real(k, kind=xp)
       p3 = ((2.0_xp*fk + 1.0_xp)*z*p2 - fk*p1) / (fk + 1.0_xp)
       p3d = ((2.0_xp*fk + 1.0_xp)*p2 + (2.0_xp*fk + 1.0_xp)*z*p2d - fk*p1d) / &
            (fk + 1.0_xp)
       p1 = p2
       p2 = p3
       p1d = p2d
       p2d = p3d
    end do
    pndleg = p3d
  end function pndleg
 
  subroutine dgllgl(D, DT, ZM1, ZM2, IM12, NZM1, NZM2, ND1, ND2)
    integer, intent(in) :: NZM1, NZM2, ND1, ND2
    real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
    real(kind=xp), intent(in) :: zm1(nd1), zm2(nd2), im12(nd2, nd1)
 
    real(kind=xp) eps, zp, zq
    integer :: IP, JQ, NM1
 
    if (nzm1 .eq. 1) then
       d(1, 1) = 0.0_xp
       dt(1, 1) = 0.0_xp
       return
    end if
    eps = 1.0e-6_xp
    nm1 = nzm1 - 1
    do ip = 1, nzm2
       do jq = 1, nzm1
          zp = zm2(ip)
          zq = zm1(jq)
          if ((abs(zp) .lt. eps) .and. (abs(zq) .lt. eps)) then
             d(ip, jq) = 0.0_xp
          else
             d(ip, jq) = (pnleg(zp, nm1) / pnleg(zq, nm1) - im12(ip, jq)) / &
                  (zp - zq)
          end if
          dt(jq, ip) = d(ip, jq)
       end do
    end do
  end subroutine dgllgl
 
  subroutine dgljgj(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
    integer, intent(in) :: NPGL, NPG, ND1, ND2
    real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
    real(kind=xp), intent(in) :: zgl(nd1), zg(nd2), iglg(nd2, nd1), alpha, beta
 
    integer, parameter :: NMAX = 84
    integer, parameter :: NDD = nmax
 
    real(kind=xp) dd(ndd, ndd), dtd(ndd, ndd)
    real(kind=xp) zgd(ndd), zgld(ndd), iglgd(ndd, ndd)
    integer :: I, J
 
    if (npgl .le. 1) then
       call neko_error('DGLJGJ: Minimum number of Gauss-Lobatto points is 2')
    else if (npgl .gt. nmax) then
       write(stderr, *) 'Polynomial degree too high in DGLJGJ'
       write(stderr, *) 'Maximum polynomial degree is', nmax
       write(stderr, *) 'Here NPGL=', npgl
       call neko_error
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGLJGJ: Alpha and Beta must be greater than -1')
    end if
 
    do i = 1, npg
       zgd(i) = zg(i)
       do j = 1, npgl
          iglgd(i, j) = iglg(i, j)
       end do
    end do
    do i = 1, npgl
       zgld(i) = zgl(i)
    end do
    call dgljgjd(dd, dtd, zgld, zgd, iglgd, npgl, npg, ndd, ndd, alpha, beta)
    do i = 1, npg
       do j = 1, npgl
          d(i, j) = dd(i, j)
          dt(j, i) = dtd(j, i)
       end do
    end do
  end subroutine dgljgj
 
  subroutine dgljgjd(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
    integer, intent(in) :: NPGL, NPG, ND1, ND2
    real(kind=xp), intent(inout) :: d(nd2, nd1), dt(nd1, nd2)
    real(kind=xp), intent(in) :: zgl(nd1), zg(nd2), iglg(nd2, nd1), alpha, beta
 
    real(kind=xp) :: eps, eigval, dn
    real(kind=xp) :: pdi, pdj, pi, pj, pm1, pdm1, pm2, pdm2
    real(kind=xp) :: dz, faci, facj, const
    integer :: I, J, NGL
 
    if (npgl .le. 1) then
       call neko_error('DGLJGJD: Minimum number of Gauss-Lobatto points is 2')
    else if ((alpha .le. -1.0_xp) .or. (beta .le. -1.0_xp)) then
       call neko_error('DGLJGJD: Alpha and Beta must be greater than -1')
    end if
 
    eps = 1.0e-6_xp
 
    ngl = npgl-1
    dn = real(ngl, kind=xp)
    eigval = -dn*(dn + alpha + beta + 1.0_xp)
 
    do i = 1, npg
       do j = 1, npgl
          dz = abs(zg(i)-zgl(j))
          if (dz .lt. eps) then
             d(i, j) = (alpha*(1.0_xp + zg(i)) - beta*(1.0_xp - zg(i))) / &
                  (2.0_xp*(1.0_xp - zg(i)**2))
          else
             call jacobf(pi, pdi, pm1, pdm1, pm2, pdm2, ngl, alpha, beta, zg(i))
             call jacobf(pj, pdj, pm1, pdm1, pm2, pdm2, ngl, alpha, beta, zgl(j))
             faci = alpha*(1.0_xp + zg(i)) - beta*(1.0_xp - zg(i))
             facj = alpha*(1.0_xp + zgl(j)) - beta*(1.0_xp - zgl(j))
             const = eigval*pj + facj*pdj
             d(i, j) = ((eigval*pi + faci*pdi) * (zg(i) - zgl(j)) - &
                  (1.0_xp - zg(i)**2)*pdi) / (const*(zg(i) - zgl(j))**2)
          end if
          dt(j, i) = d(i, j)
       end do
    end do
  end subroutine dgljgjd
 
  subroutine iglm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
    integer, intent(in) :: NZ1, NZ2, ND1, ND2
    real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
    real(kind=xp), intent(in) :: z1(nd1), z2(nd2)
    real(kind=xp) :: zi
    integer :: I, J
 
    if (nz1 .eq. 1) then
       i12(1, 1) = 1.0_xp
       it12(1, 1) = 1.0_xp
       return
    end if
 
    do i = 1, nz2
       zi = z2(i)
       do j = 1, nz1
          i12(i, j) = hgl(j, zi, z1, nz1)
          it12(j, i) = i12(i, j)
       end do
    end do
  end subroutine iglm
 
  subroutine igllm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
    integer, intent(in) :: NZ1, NZ2, ND1, ND2
    real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
    real(kind=xp), intent(in) :: z1(nd1), z2(nd2)
    real(kind=xp) :: zi
    integer :: I, J
 
    if (nz1 .eq. 1) then
       i12(1, 1) = 1.0_xp
       it12(1, 1) = 1.0_xp
       return
    end if
 
    do i = 1, nz2
       zi = z2(i)
       do j = 1, nz1
          i12(i, j) = hgll(j, zi, z1, nz1)
          it12(j, i) = i12(i, j)
       end do
    end do
  end subroutine igllm
 
  subroutine igjm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
    integer, intent(in) :: NZ1, NZ2, ND1, ND2
    real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
    real(kind=xp), intent(in) :: z1(nd1), z2(nd2), alpha, beta
    real(kind=xp) :: zi
    integer :: I, J
 
    if (nz1 .eq. 1) then
       i12(1, 1) = 1.0_xp
       it12(1, 1) = 1.0_xp
       return
    end if
 
    do i = 1, nz2
       zi = z2(i)
       do j = 1, nz1
          i12(i, j) = hgj(j, zi, z1, nz1, alpha, beta)
          it12(j, i) = i12(i, j)
       end do
    end do
  end subroutine igjm
 
  subroutine igljm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
    integer, intent(in) :: NZ1, NZ2, ND1, ND2
    real(kind=xp), intent(inout) :: i12(nd2, nd1), it12(nd1, nd2)
    real(kind=xp), intent(in) :: z1(nd1), z2(nd2), alpha, beta
    real(kind=xp) :: zi
    integer :: I, J
 
    if (nz1 .eq. 1) then
       i12(1, 1) = 1.0_xp
       it12(1, 1) = 1.0_xp
       return
    end if
 
    do i = 1, nz2
       zi = z2(i)
       do j = 1, nz1
          i12(i, j) = hglj(j, zi, z1, nz1, alpha, beta)
          it12(j, i) = i12(i, j)
       end do
    end do
  end subroutine igljm
end module speclib