program qshep2d_prb
!
!***********************************************************************
!
!! QSHEP2D_PRB is a test for QSHEP2D.
!
  implicit none
!
  integer, parameter :: n = 36
  integer, parameter :: nq = 13
  integer, parameter :: nr = 3
  integer, parameter :: nw = 19
!
  real a(5,n)
  real dx
  real dy
  real eps
  real eq
  real eqx
  real eqy
  real f(n)
  real fq
  real fx
  real fy
  integer i
  integer ier
  integer j
  integer k
  integer lcell(3,3)
  integer lnext(n)
  real p(10)
  real px
  real py
  real q
  real q1
  real qs2val
  real qx
  real qy
  real rmax
  real rq
  real rsq(n)
  real x(n)
  real xmin
  real xx
  real y(n)
  real ymin
  real yy
!
!  Quadratic test function and partial derivatives.
!
  fq(xx,yy) = (         ( xx + 2.0E+00 * yy ) / 3.0E+00 )**2
  fx(xx,yy) = 2.0E+00 * ( xx + 2.0E+00 * yy ) / 9.0E+00
  fy(xx,yy) = 4.0E+00 * ( xx + 2.0E+00 * yy ) / 9.0E+00
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QSHEP2D_PRB'
  write ( *, '(a)' ) '  Tests for QSHEP2D.'
!
!  Generate a 6 by 6 grid of nodes in the unit square with
!  the natural ordering.
!
  k = 0
  do j = 5, 0, -1
    do i = 0, 5
      k = k + 1
      x(k) = real ( i ) / 5.0E+00
      y(k) = real ( j ) / 5.0E+00
    end do
  end do
!
!  Compute the data values.
!
  do k = 1, n
    f(k) = fq ( x(k), y(k) )
  end do
!
!  Call QSHEP2 to define the interpolant Q to the data.
!
  call qshep2 ( n, x, y, f, nq, nw, nr, lcell, lnext, xmin, ymin, &
    dx, dy, rmax, rsq, a, ier )

  if ( ier /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QSHEP2D_PRB - Error!'
    write ( *, '(a,i6)' ) '  Error in QSHEP2D, IER = ', ier
    stop
  end if
!
!  Generate a 10 by 10 uniform grid of interpolation points
!  (p(i),p(j)) in the unit square.  
!
  do i = 1, 10
    p(i) = real ( i - 1 ) / 9.0E+00
  end do
!
!  Compute the machine precision EPS.
!
  eps = epsilon ( eps )
!
!  Compute the interpolation errors.
!
  eq = 0.0E+00
  eqx = 0.0E+00
  eqy = 0.0E+00

  do j = 1, 10

    py = p(j)

    do i = 1, 10

      px = p(i)

      q1 = qs2val ( px, py, n, x, y, f, nr, lcell, lnext, xmin, &
        ymin, dx, dy, rmax, rsq, a )

      call qs2grd ( px, py, n, x, y, f, nr, lcell, lnext, xmin, &
        ymin, dx, dy, rmax, rsq, a, q, qx, qy, ier )

      if ( ier /= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'QSHEP2D_PRB - Error!'
        write ( *, '(a,i6)' ) '  Error in QS2GRD, IER = ', ier
        stop
      end if

      if ( abs ( q1 - q ) > 3.0E+00 * abs ( q ) * eps ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'QSHEP2D_PRB - Error!'
        write ( *, '(a)' ) '  The interpolated values Q1 (QS2VAL)'
        write ( *, '(a)' ) '  and Q (QS2GRD) differ.'
        write ( *, '(a,g14.6)' ) '  Q1 = ', q1
        write ( *, '(a,g14.6)' ) '  Q  = ', q
        stop
      end if

      eq = max ( eq, abs ( fq ( px, py ) - q ) )
      eqx = max ( eqx, abs ( fx ( px, py ) - qx ) )
      eqy = max ( eqy, abs ( fy ( px, py ) - qy ) )

    end do

  end do
!
!  Print the maximum errors and the ratio EQ / EPS.
!
  rq = eq / eps

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Maximum absolute errors in the interpolant Q and partial'
  write ( *, '(a)' ) 'derivatives QX and QY relative to machine precision EPS.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Function   Max error   Max error/EPS'
  write ( *, '(a)' ) ' '
  write ( *, '(a4,2x,e9.3,2x,f4.2)' ) 'Q   ', eq, rq
  write ( *, '(a4,2x,e9.3)'      ) 'dQdX', eqx
  write ( *, '(a4,2x,e9.3)'      ) 'dQdY', eqy

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QSHEP2D_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end