lgwt Subroutine

    subroutine lgwt(a, b, num_pts, x, w)

        !*  This function is a fortran90 port of the matlab function,
        !   lgwt.m The source code of lgwt.m was originally found at:
        !   http://www.mathworks.com/matlabcentral/fileexchange/4540

        ! Variables in/out
        integer, intent(in)   :: num_pts
        real(wp), intent(in)  :: a, b
        real(wp), intent(out) :: x(:), w(:)

        ! Local variables
        integer:: ii, jj, N, N1, N2
        ! real(wp), parameter:: eps=sqrt(epsilon(1.0_wp))
        real(wp), parameter                       :: eps=1d-10
        real(wp), dimension(:),     allocatable   :: xu, array1, y, y0, Lpp
        real(wp), dimension(:, :),  allocatable   :: L, Lp
        real(wp), parameter                       :: pi = 4.0_wp*datan(1.0_wp)

        N = num_pts - 1
        N1 = N + 1
        N2 = N + 2

        ! Allocate and initialize arrays
        allocate(xu(N1), array1(N1), y(N1), y0(N1))
        allocate(L(N1, N2), Lp(N1, N2), Lpp(N1))

        L = 0.0_wp
        Lp = 0.0_wp

        call linspace(-1.0_wp, 1.0_wp, xu)

        array1 = [ (ii, ii = 0, N) ]
        ! Initial guess of the roots of the Legendre polynomial of order N
        y = cos((2.0_wp * dble(array1) + 1.0_wp) * &
            pi / (2.0_wp * dble(N) + 2.0_wp)) + &
            (0.27_wp/dble(N1)) * sin(pi*xu*dble(N)/dble(N2))

        y0 = 2.0_wp

        outer: do
            L(:, 1) = 1.0_wp
            Lp(:, 1) = 0.0_wp

            L(:, 2) = y
            Lp(:, 2) = 1.0_wp

            inner: do jj = 2, N1
                L(:, jj+1) = ( &
                    (2.0_wp*dble(jj)-1.0_wp) * y * L(:, jj) - &
                    dble(jj-1) * L(:, jj-1)) &
                    / dble(jj)
            enddo inner

            Lpp = dble(N2) * (L(:,N1) - y * L(:,N2)) &
                / (1.0_wp - y**2.0_wp)

            y0 = y
            y = y0 - L(:,N2)/Lpp

            if ( norm2(y-y0) < eps ) then
                exit outer
            endif
        enddo outer

        x = ( a*(1.0_wp-y) + b*(1.0_wp+y) ) / 2.0_wp
        w = ( b-a ) / ((1.0_wp - y**2.0_wp)*Lpp**2.0_wp) * &
            (dble(N2) / dble(N1))**2.0_wp

        return
    end subroutine lgwt