program test_gaussquad

  use gaussquad
  implicit none

  integer, parameter :: MAXN = 20
  real(kind=WP) :: t(MAXN), w(MAXN), work(MAXN), a(2*MAXN), b(2*MAXN), muzero, &
                   p(MAXN,0:2*MAXN), alpha, beta
  integer :: icode
  
  print "('Testing plain Gauss rules with up to ', i0, ' points:', /)", MAXN

  alpha = HALF
  beta  = TWO
  do icode = LEGENDRE, HERMITE
    call test_rule(icode, 'N')
  end do

  print "(/, 'Testing left Radau rules with up to ', i0, ' points:', /)", MAXN
  do icode = LEGENDRE, LAGUERRE
    call test_rule(icode, 'L')
  end do

  print "(/, 'Testing right Radau rules with up to ', i0, ' points:', /)", MAXN
  do icode = LEGENDRE, JACOBI
    call test_rule(icode, 'R')
  end do

  print "(/, 'Testing Lobatto rules with up to ', i0, ' points:', /)", MAXN
  do icode = LEGENDRE, JACOBI
    call test_rule(icode, 'B')
  end do

contains

  subroutine test_rule(icode, endpt)
    integer,          intent(in) :: icode
    character(len=1), intent(in) :: endpt

    integer          :: info, n, dop
    real(kind=WP)    :: d, maxd


    call recur_coeffs(icode, 2*MAXN, alpha, beta, a, b, muzero, info)
    maxd = ZERO
    select case(endpt)
    case('N')
      do n = 1, MAXN
        dop = 2*n - 1
        call gauss_rule(icode, n, t, w, work, alpha, beta, endpt, info)
        call orthonormal_polynomials(2*n, n, t, a, b, muzero, p, MAXN)
        d = discrepancy(n, dop, w, p(1:n,0:2*n))
        maxd = max(maxd, d)
      end do
    case('B')
      do n = 2, MAXN
        dop = 2*n - 3
        call gauss_rule(icode, n, t, w, work, alpha, beta, endpt, info)
        call orthonormal_polynomials(2*n, n, t, a, b, muzero, p, MAXN)
        d = discrepancy(n, dop, w, p(1:n,0:2*n))
        maxd = max(maxd, d)
      end do
    case('L','R')
      do n = 1, MAXN
        dop = 2*n - 2
        call gauss_rule(icode, n, t, w, work, alpha, beta, endpt, info)
        call orthonormal_polynomials(2*n, n, t, a, b, muzero, p, MAXN)
        d = discrepancy(n, dop, w, p(1:n,0:2*n))
        maxd = max(maxd, d)
      end do
    end select
    print "(t5, a, t28, es12.4)", RULE_NAME(icode), maxd

  end subroutine test_rule

  function discrepancy(n, dop, w, p) result(d)
    integer,       intent(in) :: n, dop
    real(kind=WP), intent(in) :: w(n), p(n,0:2*n)
    real(kind=WP)             :: d
    !
    ! Returns the maximum error in the inner product <p_j,p_k>
    ! over 0 <= j + k <= dop.
    !
    integer       :: k, j, l
    real(kind=WP) :: s

    d = ZERO
    do k = 0, dop
      do j = 0, min(k-1, dop-k)
        s = ZERO
        do l = 1, n
          s = s + w(l) * p(l,j) * p(l,k)
        end do
!        print *, j, k, "  s = ", s
        d = max(d, abs(s))
      end do
      if ( 2*k <= dop ) then
        s = ZERO
        do l = 1, n
          s = s + w(l) * p(l,k)**2
        end do
!        print *, k, k, "  s = ", s
        d = max(d, abs(s-ONE))
      end if
    end do

  end function discrepancy

end program test_gaussquad
