JacobiGQ Subroutine

    subroutine JacobiGQ(alpha, beta, n, x, w)
    !! Compute the n'th order Gauss quadrature points, \(x\) and weights,
    !! \(w\), associated with the Jacobi polynomial of type
    !! \((\alpha,\beta)>-1\).
      real(wp), intent(in) :: alpha
      !! The value of \(\alpha\).
      real(wp), intent(in) :: beta
      !! The value of \(\beta\).
      integer(ip), intent(in) :: n
      !! The order of the Jacobi polynomial.
      real(wp), dimension(n+1), intent(out) :: x
      !! On output the location of the Gauss quadrature points.
      real(wp), dimension(n+1), intent(out) :: w
      !! On output the weights.

      real(wp), dimension(1:n+1) :: h1, jdiag1, ireal
      real(wp), dimension(1:n) :: jdiag2, d1, d2, d3
      real(wp), dimension(max(1,2*n)) :: work
      real(wp), dimension(1:n+1,1:n+1) :: vect
      real(wp) :: lngammaab, lngammaa, lngammab
      integer(ip) :: info, i

      if (n==0) then
        x(1) = ( alpha - beta ) / ( alpha + beta + 2.0_wp )
        w(1) = 2.0_wp
        return
      end if

      forall(i=0:n) ireal(i+1) = real(i,wp)
      h1 = 2.0_wp*ireal+alpha+beta
      where (h1>0.0_wp)
        jdiag1 = -(alpha**2-beta**2)/(h1*(h1+2.0_wp))
      elsewhere
        jdiag1 = 0.0_wp
      end where
      d1 = 2.0_wp/(h1(1:n)+2.0_wp)
      d2 = ireal(2:n+1)*(ireal(2:n+1)+alpha+beta)* &
                        (ireal(2:n+1)+alpha)*(ireal(2:n+1)+beta)
      d3 = 1.0_wp/((h1(1:n)+1)*(h1(1:n)+3))
      jdiag2 = d1*sqrt(d2*d3)
      if (wp == dp ) then
        call dsteqr('I', n+1, jdiag1, jdiag2, vect, n+1, work, info )
      else if (wp == qp ) then
        call qsteqr('I', n+1, jdiag1, jdiag2, vect, n+1, work, info )
      end if
      if (info <0) then
        print*,'Parameter ', i, ' in call to dsteqr has illegal value'
        stop
      end if
      if (info >0) then
        print*, i, ' off-diagonal elements have not converged to zero in call to dsteqr'
        stop
      end if
      x = jdiag1
      lngammaab = log_gamma(alpha+beta+1.0_wp)
      lngammaa  = log_gamma(alpha+1.0_wp)
      lngammab  = log_gamma(beta+1.0_wp)
      w = vect(1,:)**2*2.0_wp**(alpha+beta+1)/(alpha+beta+1)* &
                       exp(lngammaa+lngammab-lngammaab)
      return
    end subroutine JacobiGQ