do i = 5, 1, -1

      t = 0.0D+00

      do j = i+1, 5
        t = t + b(j,i) * a(j,k)
      end do

      a(i,k) = ( b(6,i) - t ) / b(i,i)

    end do
!
!  Scale the coefficients to adjust for the column scaling.
!
    a(1:3,k) = a(1:3,k) / avsq
    a(4,k) = a(4,k) / av
    a(5,k) = a(5,k) / av
!
!  Unmark K and the elements of NPTS.
!
    lnext(k) = - lnext(k)

    do i = 1, lnp
      np = npts(i)
      lnext(np) = - lnext(np)
    end do

  end do
!
!  No errors encountered.
!
  xmin = xmn
  ymin = ymn
  dx = ddx
  dy = ddy
  rmax = sqrt ( rsmx )
  ier = 0

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

!***********************************************************************
!
!! QS2VAL evaluates the interpolant function at a point.
!
!  Discussion:
!
!    QS2VAL returns the value Q(PX,PY) where Q is the weighted sum of 
!    quadratic nodal functions defined by QSHEP2.  If the spatial 
!    derivatives of Q are also desired, call QS2GRD instead.
!
!    Input parameters are not altered by this function.  The
!    parameters other than PX and PY should be input unaltered
!    from their values on output from QSHEP2.  This function
!    should not be called if a nonzero error flag was returned
!    by QSHEP2.
!
!  Modified:
!
!    10 July 1999
!
!  Author:
!
!    Robert Renka,
!    University of North Texas
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 149-150.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) PX, PY, the (X,Y) coordinates of the point P at
!    which Q is to be evaluated.
!
!    Input, integer N, the number of nodes and data values to be 
!    interpolated.  N must be at least 6.
!
!    Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes at which
!    data has been supplied.
!
!    Input, real ( kind = 8 ) F(N), the data values at the nodes.
!
!    Input, integer NR, the number of rows and columns in the cell grid.
!    Refer to subroutine STORE2.  NR must be at least 1.
!
!    Input, integer LCELL(NR,NR), the array of nodal indices associated
!    with cells.  Refer to STORE2.
!
!    Input, integer LNEXT(N), the next-node indices.  Refer to STORE2.
!
!    Input, real ( kind = 8 ) XMIN, YMIN, DX, DY, the minimum nodal X, Y
!    coordinates, and the X, Y dimensions of a cell.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) RMAX, the square root of the largest element
!    in RSQ, the maximum radius of influence.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) RSQ(N), the squared radii which enter into the
!    weights defining the interpolant Q.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) A(5,N), the coefficients for the nodal functions 
!    defining the interpolant Q.  Computed by QSHEP2.
!
!    Output, real ( kind = 8 ) QS2VAL, the interpolated function value
!    at (PX,PY).
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) a(5,n)
  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  real ( kind = 8 ) f(n)
  integer i
  integer imax
  integer imin
  integer j
  integer jmax
  integer jmin
  real ( kind = 8 ) ds
  integer k
  integer kp
  integer lcell(nr,nr)
  integer lnext(n)
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) qs2val
  real ( kind = 8 ) rd
  real ( kind = 8 ) rds
  real ( kind = 8 ) rmax
  real ( kind = 8 ) rs
  real ( kind = 8 ) rsq(n)
  real ( kind = 8 ) sw
  real ( kind = 8 ) swq
  real ( kind = 8 ) w
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xp
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymin
  real ( kind = 8 ) yp

  xp = px
  yp = py
  qs2val = 0.0D+00

  if ( n < 6  ) then
    return  
  else if ( nr < 1  ) then
    return
  else if ( dx <= 0.0D+00 ) then
    return
  else if ( dy <= 0.0D+00 ) then
    return
  else if ( rmax < 0.0D+00 ) then
    return
  end if
!
!  Set imin, imax, jmin, and jmax to cell indices defining
!  the range of the search for nodes whose radii include
!  p.  The cells which must be searched are those intersected
!  by (or contained in) a circle of radius rmax
!  centered at P.
!
  imin = int ( ( xp - xmin - rmax ) / dx ) + 1
  imin = max ( imin, 1 )

  imax = int ( ( xp - xmin + rmax ) / dx ) + 1
  imax = min ( imax, nr )

  jmin = int ( ( yp - ymin - rmax ) / dy ) + 1
  jmin = max ( jmin, 1 )

  jmax = int ( ( yp - ymin + rmax ) / dy ) + 1
  jmax = min ( jmax, nr )
!
!  Test for no cells within the circle of radius RMAX.
!
  if ( imax < imin .or. jmax < jmin ) then
    qs2val = 0.0D+00
    return
  end if
!
!  Accumulate weight values in SW and weighted nodal function
!  values in swq.  the weights are w(k) = ((r-d)+/(r*d))**2
!  for r**2 = rsq(k) and d = distance between p and node K.
!
  sw = 0.0D+00
  swq = 0.0D+00

  do j = jmin, jmax

    do i = imin, imax

      k = lcell(i,j)

      if ( k /= 0 ) then

        do

          delx = xp - x(k)
          dely = yp - y(k)
          ds = delx * delx + dely * dely
          rs = rsq(k)

          if ( ds < rs ) then

            if ( ds == 0.0D+00 ) then
              qs2val = f(k)
              return
            end if

            rds = rs * ds
            rd = sqrt ( rds )
            w = ( rs + ds - rd - rd ) / rds
            sw = sw + w

            swq = swq + w * ( f(k) + a(1,k) * delx * delx &
              + a(2,k) * delx * dely + a(3,k) * dely * dely &
              + a(4,k) * delx + a(5,k) * dely )

          end if

          kp = k
          k = lnext(kp)

          if ( k == kp ) then
            exit
          end if

        end do

      end if

    end do

  end do
!
!  SW = 0 if and only if P is not within the radius R(K) for any node K.
!
  if ( sw == 0.0D+00 ) then
    qs2val = 0.0D+00
  else
    qs2val = swq / sw
  end if

  return
end
subroutine rotate ( n, c, s, x, y )

!***********************************************************************
!
!! ROTATE applies a Givens rotation.
!
!  Discussion:
!
!    The rotation has the form:
!
!      (   C  S )
!      ( - S  C )
!
!    and is essentially applied to a 2 by N matrix:
!
!      ( X(1) X(2) ... X(N) )
!      ( Y(1) Y(2) ... Y(N) )
!
!  Modified:
!
!    28 June 1999
!
!  Parameters:
!
!    Input, integer N, the dimension of the vectors.
!
!    Input, real ( kind = 8 ) C, S, the cosine and sine entries of the Givens
!    rotation matrix.  These may be determined by subroutine GIVENS.
!
!    Input/output, real ( kind = 8 ) X(N), Y(N), the rotated vectors. 
!
  implicit none

  integer n

  real ( kind = 8 ) c
  integer i
  real ( kind = 8 ) s
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xi
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) yi

  if ( n <= 0 ) then
    return
  else if ( c == 1.0D+00 .and. s == 0.0D+00 ) then
    return
  end if

  do i = 1, n
    xi = x(i)
    yi = y(i)
    x(i) =   c * xi + s * yi
    y(i) = - s * xi + c * yi
  end do

  return
end
subroutine setup2 ( xk, yk, fk, xi, yi, fi, s1, s2, r, row )

!***********************************************************************
!
!! SETUP2 sets up a row of the least squares regression matrix.
!
!  Discussion:
!
!    SETUP2 sets up the I-th row of an augmented regression matrix for 
!    a weighted least-squares fit of a quadratic function Q(X,Y) to a set 
!    of data values F, where Q(XK,YK) = FK.  
!
!    The first 3 columns are quadratic terms, and are scaled by 1/S2.
!    The fourth and fifth columns represent linear terms, and are scaled 
!    by 1/S1.  
!
!    If D = 0, or R <= D, the weight is
!      0,
!    else if D < R, the weight is 
!      (R-D)/(R*D), 
!    where D is the distance between nodes I and K, and R is a maximum
!    influence distance.
!
!  Modified:
!
!    05 July 1999
!
!  Author:
!
!    Robert Renka,
!    University of North Texas
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 149-150.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XK, YK, FK, the coordinates and value of the data
!    at data node K.
!
!    Input, real ( kind = 8 ) XI, YI, FI, the coorindates and value of the data
!    at data node I.
!
!    Input, real ( kind = 8 ) S1, S2, reciprocals of the scale factors.
!
!    Input, real ( kind = 8 ) R, the maximum radius of influence about node K.
!
!    Output, real ( kind = 8 ) ROW(6), a row of the augmented regression matrix.
!
  implicit none

  real ( kind = 8 ) d
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  real ( kind = 8 ) fi
  real ( kind = 8 ) fk
  integer i
  real ( kind = 8 ) r
  real ( kind = 8 ) row(6)
  real ( kind = 8 ) s1
  real ( kind = 8 ) s2
  real ( kind = 8 ) w
  real ( kind = 8 ) xi
  real ( kind = 8 ) xk
  real ( kind = 8 ) yi
  real ( kind = 8 ) yk

  dx = xi - xk
  dy = yi - yk

  d = sqrt ( dx * dx + dy * dy )

  if ( d <= 0.0D+00 .or. r <= d ) then

    row(1:6) = 0.0D+00

  else

    w = ( r - d ) / r / d

    row(1) = dx * dx * w / s2
    row(2) = dx * dy * w / s2
    row(3) = dy * dy * w / s2
    row(4) = dx * w / s1
    row(5) = dy * w / s1
    row(6) = ( fi - fk ) * w

  end if

  return
end
subroutine store2 ( n, x, y, nr, lcell, lnext, xmin, ymin, dx, dy, ier )

!***********************************************************************
!
!! STORE2 creates a cell data structure for the scattered data.
!
!  Discussion:
!
!    STORE2 is given a set of N arbitrarily distributed nodes in the 
!    plane and creates a data structure for a cell-based method of 
!    solving closest-point problems.  The smallest rectangle containing 
!    all the nodes is partitioned into an NR by NR uniform grid of cells, 
!    and nodes are associated with cells.      
!
!    In particular, the data structure stores the indices of the nodes 
!    contained in each cell.  For a uniform random distribution of nodes, 
!    the nearest node to an arbitrary point can be determined in constant
!    expected time.
!
!  Modified:
!
!    05 July 1999
!
!  Author:
!
!    Robert Renka
!    University of North Texas
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 149-150.
!
!  Parameters:
!
!    Input, integer N, the number of data nodes.  N must be at least 2.
!
!    Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the data nodes.
!
!    Input, integer NR, the number of rows and columns in the grid.  The
!    cell density, or average number of data nodes per cell, is
!      D = N / ( NR * NR ).
!    A recommended value, based on empirical evidence, is 
!      D = 3. 
!    Hence, the corresponding value of NR is recommended to be about
!      NR = SQRT ( N / 3 ).  
!    NR must be at least 1.
!
!    Output, integer LCELL(NR,NR), an array set up so that LCELL(I,J)
!    contains the index (for X and Y) of the first data node (that is, the
!    data node with smallest index) in the (I,J) cell.  LCELL(I,J) will be 0 if 
!    no data nodes are contained in the (I,J) cell.  The upper right corner of 
!    the (I,J) cell has coordinates 
!      ( XMIN + I * DX, YMIN + J * DY ).
!
!    Output, integer LNEXT(N), an array of next-node indices.  LNEXT(K)
!    contains the index of the next node in the cell which contains node K, 
!    or LNEXT(K) = K if K is the last node in the cell.
!    The data nodes contained in a cell are ordered by their indices.
!    If, for example, cell (I,J) contains nodes 2, 3, and 5 and no others, 
!    then:
!
!      LCELL(I,J) = 2, (index of the first data node)
!
!      LNEXT(2) = 3, 
!      LNEXT(3) = 5,
!      LNEXT(5) = 5.
!
!    Output, real ( kind = 8 ) XMIN, YMIN, the X, Y coordinates of the lower
!    left corner of the rectangle defined by the data nodes.  The upper right 
!    corner is ( XMAX, YMAX ), where
!      XMAX = XMIN + NR * DX,
!      YMAX = YMIN + NR * DY.
!
!    Output, real ( kind = 8 ) DX, DY, the X and Y dimensions of the 
!    individual cells.
!      DX = ( XMAX - XMIN ) / NR
!      DY = ( YMAX - YMIN ) / NR,
!    where XMIN, XMAX, YMIN and YMAX are the extrema of X and Y.
!
!    Output, integer IER, an error indicator.
!    0, if no errors were encountered.
!    1, if N < 2 or NR < 1.
!    2, if DX = 0 or DY = 0.
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  integer i
  integer ier
  integer j
  integer k
  integer l
  integer lcell(nr,nr)
  integer lnext(n)
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin

  ier = 0

  if ( n < 2 ) then
    ier = 1
    return
  end if

  if ( nr < 1 ) then
    ier = 1
    return
  end if
!
!  Compute the dimensions of the (X,Y) rectangle containing all the data nodes.
!
  xmin = minval ( x(1:n) )
  xmax = maxval ( x(1:n) )
  ymin = minval ( y(1:n) )
  ymax = maxval ( y(1:n) )
!
!  Compute the dimensions of a single cell.
!
  dx = ( xmax - xmin ) / real ( nr, kind = 8 )
  dy = ( ymax - ymin ) / real ( nr, kind = 8 )
!
!  Test for zero area.
!
  if ( dx == 0.0D+00 .or. dy == 0.0D+00 ) then
    ier = 2
    return
  end if
!
!  Initialize LCELL.
!
  lcell(1:nr,1:nr) = 0
!
!  Loop on nodes, storing indices in LCELL and LNEXT.
!
  do k = n, 1, -1

    i = int ( ( x(k) - xmin ) / dx ) + 1
    i = min ( i, nr )

    j = int ( ( y(k) - ymin ) / dy ) + 1
    j = min ( j, nr )

    l = lcell(i,j)

    if ( l /= 0 ) then
      lnext(k) = l
    else
      lnext(k) = k
    end if

    lcell(i,j) = k

  end do

  return
end