Description:
     Subroutine pda_bispev evaluates on a grid (x(i),y(j)),i=1,...,mx;
     j=1,...,my a bivariate spline s(x,y) of degrees kx and ky, given in
     the b-spline representation.
  
  Calling Sequence:
     call pda_bispev( tx, nx, ty, ny, c, kx, ky, x, mx, y, my, z, wrk, lwrk,
                      iwrk, kwrk, ier )
  
  Input Parameters:
     tx    : Real array, length nx, which contains the position of the
             knots in the x-direction.
  
     nx    : Integer, giving the total number of knots in the x-direction
  
     ty    : Real array, length ny, which contains the position of the
             knots in the y-direction.
  
     ny    : Integer, giving the total number of knots in the y-direction
  
     c     : Real array, length (nx-kx-1)*(ny-ky-1), which contains the
             b-spline coefficients.
  
     kx,ky : Integer values, giving the degrees of the spline.
  
     x     : Real array of dimension (mx). Before entry x(i) must be set to
             the x co-ordinate of the i-th grid point along the x-axis.
             tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
  
     mx    : Integer. On entry mx must specify the number of grid points along
             the x-axis. Mx >=1.
  
     y     : Real array of dimension (my). Before entry y(j) must be set to
             the y co-ordinate of the j-th grid point along the y-axis.
             ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
  
     my    : Integer. On entry my must specify the number of grid points along
             the y-axis. My >=1.
  
     wrk   : Real array of dimension lwrk. Used as workspace.
  
     lwrk  : Integer, specifying the dimension of wrk.
             lwrk >= mx*(kx+1)+my*(ky+1)
  
     iwrk  : Integer array of dimension kwrk. Used as workspace.
  
     kwrk  : Integer, specifying the dimension of iwrk. Kwrk >= mx+my.
  
  Output Parameters:
     z     : Real array of dimension (mx*my). On successful exit z(my*(i-1)+j)
             contains the value of s(x,y) at the point
             (x(i),y(j)),i=1,...,mx;j=1,...,my.
  
     ier   : Integer error flag:
  
             0 : Normal return.
            10 : Invalid input data (see restrictions).
  
  Restrictions:
     mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
     tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
     ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
  
  Other Subroutines Required:
     pda_fpbisp, pda_fpbspl
  
  References :
     de Boor C : "On calculating with b-splines", j. Approximation theory
                  6 (1972) 50-62.
     Cox M.G.  : "The numerical evaluation of b-splines", j. inst. maths
                  applics 10 (1972) 134-149.
     Dierckx P.: "Curve and surface fitting with splines", monographs on
                   numerical analysis, Oxford University Press, 1993.
  
  Author :
     P. Dierckx
     Dept. Computer Science, k.u.leuven
     celestijnenlaan 200a, b-3001 heverlee, Belgium.
     e-mail : Paul.Dierckx@cs.kuleuven.ac.be
  
  Latest Update : march 1987