program qshep3d_prb
!
!***********************************************************************
!
!! QSHEP3D_PRB tests QSHEP3D.
!
!      algorithm 661, collected algorithms from acm.
!      this work published in transactions on mathematical software,
!      vol. 14, no. 2, p.151.
!
!         qs3test
!
!   this program tests the scattered data interpolation
! package qshep3d by printing the maximum errors associated
! with interpolated values and gradients on a 5 by 5 by 5
! uniform grid in the unit cube.  the data set consists
! of 64 nodes with data values taken from a quadratic func-
! tion for which the method is exact.  the ratio of maximum
! interpolation error relative to the machine precision is
! also printed.  this should be o(1).  the interpolated
! values from qs3val and qs3grd are compared for agreement.
!
  implicit none
!
  integer, parameter :: n = 64
  integer, parameter :: nq = 17
  integer, parameter :: nr = 3
  integer, parameter :: nw = 32
!
  real a(9,n)
  real eps
  real eq
  real eqx
  real eqy
  real eqz
  real f(n)
  real fq
  real fx
  real fy
  real fz
  integer i
  integer ier
  integer j
  integer k
  integer l
  integer lcell(3,3,3)
  integer lnext(n)
  real p(5)
  real px
  real py
  real pz
  real q
  real q1
  real qs3val
  real qx
  real qy
  real qz
  real rmax
  real rq
  real rsq(n)
  real x(n)
  real xx
  real xyzdel(3)
  real xyzmin(3)
  real y(n)
  real yl
  real yy
  real z(n)
  real zl
  real zz
!
!  Quadratic test function and its partial derivatives.
!
  fq(xx,yy,zz) = ( ( xx + 2.0 * yy + 3.0 * zz ) / 6.0 )**2
  fx(xx,yy,zz) = ( xx + 2.0 * yy + 3.0 * zz ) / 18.0
  fy(xx,yy,zz) = ( xx + 2.0 * yy + 3.0 * zz ) / 9.0
  fz(xx,yy,zz) = ( xx + 2.0 * yy + 3.0 * zz ) / 6.0
!
  call timestamp ( )

  write ( *, * ) ' '
  write ( *, * ) 'QSHEP3D_PRB'
  write ( *, * ) '  A set of tests for QSHEP3D.'
!
!  Generate a 4 by 4 by 4 grid of nodes in the unit cube.
!
  l = 0
  do k = 1,4
    zl = real ( k - 1 ) / 3.0E+00
    do j = 1, 4
      yl = real ( j - 1 ) / 3.0E+00
      do i = 1, 4
        l = l + 1
        x(l) = real ( i - 1 ) / 3.0E+00
        y(l) = yl
        z(l) = zl
      end do
    end do
  end do
!
!  Compute the data values.
!
  do l = 1, n
    f(l) = fq ( x(l), y(l), z(l) )
  end do
!
!  Compute parameters defining the interpolant Q.
!
  call qshep3 ( n, x, y, z, f, nq, nw, nr, lcell, lnext, xyzmin, &
    xyzdel, rmax, rsq, a, ier )

  if ( ier /= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'QSHEP3D_PRB'
    write ( *, * ) '  Error return from QSHEP3, IER = ', ier
    stop
  end if
!
!  Generate a 5 by 5 by 5 uniform grid of interpolation
!  points (p(i),p(j),p(k)) in the unit cube.  The eight
!  corners coincide with nodes.
!
  do i = 1, 5
    p(i) = real ( i - 1 ) / 4.0
  end do
!
!  Compute the machine precision eps.
!
  eps = epsilon ( eps )
!
!  Compute interpolation errors and test for agreement in the
!  Q values returned by qs3val and qs3grd.
!
  eq = 0.0
  eqx = 0.0
  eqy = 0.0
  eqz = 0.0
  do k = 1, 5
    pz = p(k)
    do j = 1, 5
      py = p(j)
      do i = 1, 5
        px = p(i)

        q1 = qs3val ( px, py, pz, n, x, y, z, f, nr, lcell, lnext, &
          xyzmin, xyzdel, rmax, rsq, a )

        call qs3grd ( px, py, pz, n, x, y, z, f, nr, lcell, lnext, &
          xyzmin, xyzdel, rmax, rsq, a, q, qx, qy, qz, ier )

        if ( ier /= 0 ) then
          write ( *, * ) ' '
          write ( *, * ) 'QSHEP3D_PRB'
          write ( *, * ) '  Error return from QS3GRD, IER = ', ier
          stop
        end if

        if ( abs ( q1 - q ) > 3.0 * abs ( q ) * eps ) then
          write ( *, * ) ' '
          write ( *, * ) 'QSHEP3D_PRB - Error.'
          write ( *, * ) '  The interpolated values Q1 (QS3VAL) and'
          write ( *, * ) '  Q (QS3GRD) differ.'
          write ( *, * ) '  Q1 = ', q1
          write ( *, * ) '  Q  = ', q
          stop
        end if

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

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

  write ( *, * ) ' '
  write ( *, * ) 'Maximum absolute errors in the interpolant Q '
  write ( *, * ) 'and partial derivatives (Qx,Qy,Qz) relative '
  write ( *, * ) 'to machine precision EPS'
  write ( *, * ) ' '
  write ( *, * ) '  Function   Max error  Max error/EPS'
  write ( *, * ) ' '
  write ( *, '(''      Q       '',e9.3,''       '',f4.2)' ) eq, rq
  write ( *, '(''      Qx      '',e9.3)' ) eqx
  write ( *, '(''      Qy      '',e9.3)' ) eqy
  write ( *, '(''      Qz      '',e9.3)' ) eqz

  write ( *, * ) ' '
  write ( *, * ) 'QSHEP3D_PRB'
  write ( *, * ) '  Normal end of execution.'

  stop
end