db2ink Subroutine

public pure subroutine db2ink(x, nx, y, ny, fcn, kx, ky, iknot, tx, ty, bcoef, iflag)

Determines the parameters of a function that interpolates the two-dimensional gridded data The interpolating function and its derivatives may subsequently be evaluated by the function db2val.

The interpolating function is a piecewise polynomial function represented as a tensor product of one-dimensional b-splines. the form of this function is

where the functions and are one-dimensional b-spline basis functions. the coefficients are chosen so that

Note that for each fixed value of , is a piecewise polynomial function of alone, and for each fixed value of , is a piecewise polynomial function of alone. in one dimension a piecewise polynomial may be created by partitioning a given interval into subintervals and defining a distinct polynomial piece on each one. the points where adjacent subintervals meet are called knots. each of the functions and above is a piecewise polynomial.

Users of db2ink choose the order (degree+1) of the polynomial pieces used to define the piecewise polynomial in each of the and directions (kx and ky). users also may define their own knot sequence in and separately (tx and ty). if iflag=0, however, db2ink will choose sequences of knots that result in a piecewise polynomial interpolant with kx-2 continuous partial derivatives in and ky-2 continuous partial derivatives in . (kx knots are taken near each endpoint in the direction, not-a-knot end conditions are used, and the remaining knots are placed at data points if kx is even or at midpoints between data points if kx is odd. the direction is treated similarly.)

After a call to db2ink, all information necessary to define the interpolating function are contained in the parameters nx, ny, kx, ky, tx, ty, and bcoef. These quantities should not be altered until after the last call of the evaluation routine db2val.

History

  • Boisvert, Ronald, NBS : 25 may 1982 : Author of original routine.
  • JEC : 000330 modified array declarations.
  • Jacob Williams, 2/24/2015 : extensive refactoring of CMLIB routine.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(:) :: x

(nx) array of abcissae. Must be strictly increasing.

integer(kind=ip), intent(in) :: nx

Number of abcissae

real(kind=wp), intent(in), dimension(:) :: y

(ny) array of abcissae. Must be strictly increasing.

integer(kind=ip), intent(in) :: ny

Number of abcissae

real(kind=wp), intent(in), dimension(:,:) :: fcn

(nx,ny) matrix of function values to interpolate. fcn(i,j) should contain the function value at the point (x(i),y(j))

integer(kind=ip), intent(in) :: kx

The order of spline pieces in ( ) (order = polynomial degree + 1)

integer(kind=ip), intent(in) :: ky

The order of spline pieces in ( ) (order = polynomial degree + 1)

integer(kind=ip), intent(in) :: iknot

knot sequence flag:

  • 0 = knot sequence chosen by db1ink.
  • 1 = knot sequence chosen by user.
real(kind=wp), intent(inout), dimension(:) :: tx

The (nx+kx) knots in the direction for the spline interpolant.

  • If iknot=0 these are chosen by db2ink.
  • If iknot=1 these are specified by the user.

Must be non-decreasing.

real(kind=wp), intent(inout), dimension(:) :: ty

The (ny+ky) knots in the direction for the spline interpolant.

  • If iknot=0 these are chosen by db2ink.
  • If iknot=1 these are specified by the user.

Must be non-decreasing.

real(kind=wp), intent(out), dimension(:,:) :: bcoef

(nx,ny) matrix of coefficients of the b-spline interpolant.

integer(kind=ip), intent(out) :: iflag
  • 0 = successful execution.
  • 2 = iknot out of range.
  • 3 = nx out of range.
  • 4 = kx out of range.
  • 5 = x not strictly increasing.
  • 6 = tx not non-decreasing.
  • 7 = ny out of range.
  • 8 = ky out of range.
  • 9 = y not strictly increasing.
  • 10 = ty not non-decreasing.
  • 700 = size(x) size(fcn,1)
  • 701 = size(y) size(fcn,2)
  • 706 = size(x) nx
  • 707 = size(y) ny
  • 712 = size(tx) nx+kx
  • 713 = size(ty) ny+ky
  • 800 = size(x) size(bcoef,1)
  • 801 = size(y) size(bcoef,2)

Calls

proc~~db2ink~~CallsGraph proc~db2ink bspline_sub_module::db2ink proc~check_inputs bspline_sub_module::check_inputs proc~db2ink->proc~check_inputs proc~dbknot bspline_sub_module::dbknot proc~db2ink->proc~dbknot proc~dbtpcf bspline_sub_module::dbtpcf proc~db2ink->proc~dbtpcf proc~dbintk bspline_sub_module::dbintk proc~dbtpcf->proc~dbintk proc~dbnslv bspline_sub_module::dbnslv proc~dbtpcf->proc~dbnslv proc~dbintk->proc~dbnslv proc~dbnfac bspline_sub_module::dbnfac proc~dbintk->proc~dbnfac proc~dbspvn bspline_sub_module::dbspvn proc~dbintk->proc~dbspvn

Called by

proc~~db2ink~~CalledByGraph proc~db2ink bspline_sub_module::db2ink proc~initialize_2d_auto_knots bspline_oo_module::bspline_2d%initialize_2d_auto_knots proc~initialize_2d_auto_knots->proc~db2ink proc~initialize_2d_specify_knots bspline_oo_module::bspline_2d%initialize_2d_specify_knots proc~initialize_2d_specify_knots->proc~db2ink proc~bspline_2d_constructor_auto_knots bspline_oo_module::bspline_2d_constructor_auto_knots proc~bspline_2d_constructor_auto_knots->proc~initialize_2d_auto_knots proc~bspline_2d_constructor_specify_knots bspline_oo_module::bspline_2d_constructor_specify_knots proc~bspline_2d_constructor_specify_knots->proc~initialize_2d_specify_knots interface~bspline_2d bspline_oo_module::bspline_2d interface~bspline_2d->proc~bspline_2d_constructor_auto_knots interface~bspline_2d->proc~bspline_2d_constructor_specify_knots

Source Code

    pure subroutine db2ink(x,nx,y,ny,fcn,kx,ky,iknot,tx,ty,bcoef,iflag)

    implicit none

    integer(ip),intent(in)                  :: nx     !! Number of \(x\) abcissae
    integer(ip),intent(in)                  :: ny     !! Number of \(y\) abcissae
    integer(ip),intent(in)                  :: kx     !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                  :: ky     !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)        :: x      !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in)        :: y      !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:),intent(in)      :: fcn    !! `(nx,ny)` matrix of function values to interpolate.
                                                      !! `fcn(i,j)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)                  :: iknot  !! knot sequence flag:
                                                      !!
                                                      !! * 0 = knot sequence chosen by [[db1ink]].
                                                      !! * 1 = knot sequence chosen by user.
    real(wp),dimension(:),intent(inout)     :: tx     !! The `(nx+kx)` knots in the \(x\) direction for the spline
                                                      !! interpolant.
                                                      !!
                                                      !! * If `iknot=0` these are chosen by [[db2ink]].
                                                      !! * If `iknot=1` these are specified by the user.
                                                      !!
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(inout)     :: ty     !! The `(ny+ky)` knots in the \(y\) direction for the spline
                                                      !! interpolant.
                                                      !!
                                                      !! * If `iknot=0` these are chosen by [[db2ink]].
                                                      !! * If `iknot=1` these are specified by the user.
                                                      !!
                                                      !! Must be non-decreasing.
    real(wp),dimension(:,:),intent(out)     :: bcoef  !! `(nx,ny)` matrix of coefficients of the b-spline interpolant.
    integer(ip),intent(out)                 :: iflag  !! *  0 = successful execution.
                                                      !! *  2 = `iknot` out of range.
                                                      !! *  3 = `nx` out of range.
                                                      !! *  4 = `kx` out of range.
                                                      !! *  5 = `x` not strictly increasing.
                                                      !! *  6 = `tx` not non-decreasing.
                                                      !! *  7 = `ny` out of range.
                                                      !! *  8 = `ky` out of range.
                                                      !! *  9 = `y` not strictly increasing.
                                                      !! * 10 = `ty` not non-decreasing.
                                                      !! * 700 = `size(x)`  \( \ne \) `size(fcn,1)`
                                                      !! * 701 = `size(y)`  \( \ne \) `size(fcn,2)`
                                                      !! * 706 = `size(x)`  \( \ne \) `nx`
                                                      !! * 707 = `size(y)`  \( \ne \) `ny`
                                                      !! * 712 = `size(tx)` \( \ne \) `nx+kx`
                                                      !! * 713 = `size(ty)` \( \ne \) `ny+ky`
                                                      !! * 800 = `size(x)`  \( \ne \) `size(bcoef,1)`
                                                      !! * 801 = `size(y)`  \( \ne \) `size(bcoef,2)`

    logical :: status_ok
    real(wp),dimension(:),allocatable :: temp !! work array of length `nx*ny`
    real(wp),dimension(:),allocatable :: work !! work array of length `max(2*kx*(nx+1),2*ky*(ny+1))`

    !check validity of inputs

    call check_inputs(  iknot,&
                        iflag,&
                        nx=nx,ny=ny,&
                        kx=kx,ky=ky,&
                        x=x,y=y,&
                        tx=tx,ty=ty,&
                        f2=fcn,&
                        bcoef2=bcoef,&
                        status_ok=status_ok)

    if (status_ok) then

        !choose knots
        if (iknot == 0_ip) then
            call dbknot(x,nx,kx,tx)
            call dbknot(y,ny,ky,ty)
        end if

        allocate(temp(nx*ny))
        allocate(work(max(2_ip*kx*(nx+1_ip),2_ip*ky*(ny+1_ip))))

        !construct b-spline coefficients
                         call dbtpcf(x,nx,fcn, nx,ny,tx,kx,temp, work,iflag)
        if (iflag==0_ip) call dbtpcf(y,ny,temp,ny,nx,ty,ky,bcoef,work,iflag)

        deallocate(temp)
        deallocate(work)

    end if

    end subroutine db2ink