db4ink – bspline-fortran

public pure subroutine db4ink(x, nx, y, ny, z, nz, q, nq, fcn, kx, ky, kz, kq, iknot, tx, ty, tz, tq, bcoef, iflag)

Determines the parameters of a function that interpolates the four-dimensional gridded data $[x(i),y(j),z(k),q(l),\mathrm{fcn}(i,j,k,l)] ~\mathrm{for}~ i=1,..,n_x ~\mathrm{and}~ j=1,..,n_y, ~\mathrm{and}~ k=1,..,n_z, ~\mathrm{and}~ l=1,..,n_q$ The interpolating function and its derivatives may subsequently be evaluated by the function db4val.

See db3ink header for more details.

History

Jacob Williams, 2/24/2015 : Created this routine.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(:)	::	x	`(nx)` array of $x$ abcissae. must be strictly increasing.
integer(kind=ip),	intent(in)		::	nx	number of $x$ abcissae ( $\ge 3$ )
real(kind=wp),	intent(in),	dimension(:)	::	y	`(ny)` array of $y$ abcissae. must be strictly increasing.
integer(kind=ip),	intent(in)		::	ny	number of $y$ abcissae ( $\ge 3$ )
real(kind=wp),	intent(in),	dimension(:)	::	z	`(nz)` array of $z$ abcissae. must be strictly increasing.
integer(kind=ip),	intent(in)		::	nz	number of $z$ abcissae ( $\ge 3$ )
real(kind=wp),	intent(in),	dimension(:)	::	q	`(nq)` array of $q$ abcissae. must be strictly increasing.
integer(kind=ip),	intent(in)		::	nq	number of $q$ abcissae ( $\ge 3$ )
real(kind=wp),	intent(in),	dimension(:,:,:,:)	::	fcn	`(nx,ny,nz,nq)` matrix of function values to interpolate. `fcn(i,j,k,q)` should contain the function value at the point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
integer(kind=ip),	intent(in)		::	kx	the order of spline pieces in $x$ ( $2 \le k_x < n_x$ ). (order = polynomial degree + 1)
integer(kind=ip),	intent(in)		::	ky	the order of spline pieces in $y$ ( $2 \le k_y < n_y$ ). (order = polynomial degree + 1)
integer(kind=ip),	intent(in)		::	kz	the order of spline pieces in $z$ ( $2 \le k_z < n_z$ ). (order = polynomial degree + 1)
integer(kind=ip),	intent(in)		::	kq	the order of spline pieces in $q$ ( $2 \le k_q < n_q$ ). (order = polynomial degree + 1)
integer(kind=ip),	intent(in)		::	iknot	knot sequence flag: 0 = knot sequence chosen by db4ink. 1 = knot sequence chosen by user.
real(kind=wp),	intent(inout),	dimension(:)	::	tx	The `(nx+kx)` knots in the x direction for the spline interpolant. If `iknot=0` these are chosen by db4ink. 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 y direction for the spline interpolant. If `iknot=0` these are chosen by db4ink. If `iknot=1` these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(inout),	dimension(:)	::	tz	The `(nz+kz)` knots in the z direction for the spline interpolant. If `iknot=0` these are chosen by db4ink. If `iknot=1` these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(inout),	dimension(:)	::	tq	The `(nq+kq)` knots in the q direction for the spline interpolant. If `iknot=0` these are chosen by db4ink. If `iknot=1` these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(out),	dimension(:,:,:,:)	::	bcoef	`(nx,ny,nz,nq)` 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. 11 = `nz` out of range. 12 = `kz` out of range. 13 = `z` not strictly increasing. 14 = `tz` not non-decreasing. 15 = `nq` out of range. 16 = `kq` out of range. 17 = `q` not strictly increasing. 18 = `tq` not non-decreasing. 700 = `size(x)` $\ne$ `size(fcn,1)` 701 = `size(y)` $\ne$ `size(fcn,2)` 702 = `size(z)` $\ne$ `size(fcn,3)` 703 = `size(q)` $\ne$ `size(fcn,4)` 706 = `size(x)` $\ne$ `nx` 707 = `size(y)` $\ne$ `ny` 708 = `size(z)` $\ne$ `nz` 709 = `size(q)` $\ne$ `nq` 712 = `size(tx`) $\ne$ `nx+kx` 713 = `size(ty`) $\ne$ `ny+ky` 714 = `size(tz`) $\ne$ `nz+kz` 715 = `size(tq`) $\ne$ `nq+kq` 800 = `size(x)` $\ne$ `size(bcoef,1)` 801 = `size(y)` $\ne$ `size(bcoef,2)` 802 = `size(z)` $\ne$ `size(bcoef,3)` 803 = `size(q)` $\ne$ `size(bcoef,4)`

Calls

Help

Called by

Help

Source Code

    pure subroutine db4ink(x,nx,y,ny,z,nz,q,nq,&
                           fcn,&
                           kx,ky,kz,kq,&
                           iknot,&
                           tx,ty,tz,tq,&
                           bcoef,iflag)

    implicit none

    integer(ip),intent(in)                      :: nx    !! number of \(x\) abcissae ( \( \ge 3 \) )
    integer(ip),intent(in)                      :: ny    !! number of \(y\) abcissae ( \( \ge 3 \) )
    integer(ip),intent(in)                      :: nz    !! number of \(z\) abcissae ( \( \ge 3 \) )
    integer(ip),intent(in)                      :: nq    !! number of \(q\) abcissae ( \( \ge 3 \) )
    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)
    integer(ip),intent(in)                      :: kz    !! the order of spline pieces in \(z\)
                                                         !! ( \( 2 \le k_z < n_z \) ).
                                                         !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: kq    !! the order of spline pieces in \(q\)
                                                         !! ( \( 2 \le k_q < n_q \) ).
                                                         !! (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)            :: z     !! `(nz)` array of \(z\) abcissae. must be strictly increasing.
    real(wp),dimension(:),intent(in)            :: q     !! `(nq)` array of \(q\) abcissae. must be strictly increasing.
    real(wp),dimension(:,:,:,:),intent(in)      :: fcn   !! `(nx,ny,nz,nq)` matrix of function values to interpolate.
                                                         !! `fcn(i,j,k,q)` should contain the function value at the
                                                         !!  point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
    integer(ip),intent(in)                      :: iknot !! knot sequence flag:
                                                         !!
                                                         !! * 0 = knot sequence chosen by [[db4ink]].
                                                         !! * 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 [[db4ink]].
                                                         !! * 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 [[db4ink]].
                                                         !! * If `iknot=1` these are specified by the user.
                                                         !!
                                                         !! Must be non-decreasing.
    real(wp),dimension(:),intent(inout)         :: tz    !! The `(nz+kz)` knots in the z direction for the spline
                                                         !! interpolant.
                                                         !!
                                                         !! * If `iknot=0` these are chosen by [[db4ink]].
                                                         !! * If `iknot=1` these are specified by the user.
                                                         !!
                                                         !! Must be non-decreasing.
    real(wp),dimension(:),intent(inout)         :: tq    !! The `(nq+kq)` knots in the q direction for the spline
                                                         !! interpolant.
                                                         !!
                                                         !! * If `iknot=0` these are chosen by [[db4ink]].
                                                         !! * If `iknot=1` these are specified by the user.
                                                         !!
                                                         !! Must be non-decreasing.
    real(wp),dimension(:,:,:,:),intent(out)     :: bcoef !! `(nx,ny,nz,nq)` 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.
                                                         !! * 11 = `nz` out of range.
                                                         !! * 12 = `kz` out of range.
                                                         !! * 13 = `z` not strictly increasing.
                                                         !! * 14 = `tz` not non-decreasing.
                                                         !! * 15 = `nq` out of range.
                                                         !! * 16 = `kq` out of range.
                                                         !! * 17 = `q` not strictly increasing.
                                                         !! * 18 = `tq` not non-decreasing.
                                                         !! * 700 = `size(x)`  \( \ne \) `size(fcn,1)`
                                                         !! * 701 = `size(y)`  \( \ne \) `size(fcn,2)`
                                                         !! * 702 = `size(z)`  \( \ne \) `size(fcn,3)`
                                                         !! * 703 = `size(q)`  \( \ne \) `size(fcn,4)`
                                                         !! * 706 = `size(x)`  \( \ne \) `nx`
                                                         !! * 707 = `size(y)`  \( \ne \) `ny`
                                                         !! * 708 = `size(z)`  \( \ne \) `nz`
                                                         !! * 709 = `size(q)`  \( \ne \) `nq`
                                                         !! * 712 = `size(tx`) \( \ne \) `nx+kx`
                                                         !! * 713 = `size(ty`) \( \ne \) `ny+ky`
                                                         !! * 714 = `size(tz`) \( \ne \) `nz+kz`
                                                         !! * 715 = `size(tq`) \( \ne \) `nq+kq`
                                                         !! * 800 = `size(x)`  \( \ne \) `size(bcoef,1)`
                                                         !! * 801 = `size(y)`  \( \ne \) `size(bcoef,2)`
                                                         !! * 802 = `size(z)`  \( \ne \) `size(bcoef,3)`
                                                         !! * 803 = `size(q)`  \( \ne \) `size(bcoef,4)`

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

    ! check validity of input

    call check_inputs(  iknot,&
                        iflag,&
                        nx=nx,ny=ny,nz=nz,nq=nq,&
                        kx=kx,ky=ky,kz=kz,kq=kq,&
                        x=x,y=y,z=z,q=q,&
                        tx=tx,ty=ty,tz=tz,tq=tq,&
                        f4=fcn,&
                        bcoef4=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)
            call dbknot(z,nz,kz,tz)
            call dbknot(q,nq,kq,tq)
        end if

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

        ! construct b-spline coefficients
                         call dbtpcf(x,nx,fcn,  nx,ny*nz*nq,tx,kx,temp, work,iflag)
        if (iflag==0_ip) call dbtpcf(y,ny,temp, ny,nx*nz*nq,ty,ky,bcoef,work,iflag)
        if (iflag==0_ip) call dbtpcf(z,nz,bcoef,nz,nx*ny*nq,tz,kz,temp, work,iflag)
        if (iflag==0_ip) call dbtpcf(q,nq,temp, nq,nx*ny*nz,tq,kq,bcoef,work,iflag)

        deallocate(temp)
        deallocate(work)

     end if

    end subroutine db4ink

db4ink Subroutine

Contents

Source Code

public pure subroutine db4ink(x, nx, y, ny, z, nz, q, nq, fcn, kx, ky, kz, kq, iknot, tx, ty, tz, tq, bcoef, iflag)

History

Arguments

Calls

Called by

Source Code