db5ink_default – bspline-fortran

private pure subroutine db5ink_default(x, nx, y, ny, z, nz, q, nq, r, nr, fcn, kx, ky, kz, kq, kr, iknot, tx, ty, tz, tq, tr, bcoef, iflag)

Determines the parameters of a function that interpolates the five-dimensional gridded data:

$[x(i),y(j),z(k),q(l),r(m),\mathrm{fcn}(i,j,k,l,m)]$

for:

$i=1,..,n_x ~\mathrm{and}~ j=1,..,n_y, ~\mathrm{and}~ k=1,..,n_z, ~\mathrm{and}~ l=1,..,n_q, ~\mathrm{and}~ m=1,..,n_r$

The interpolating function and its derivatives may subsequently be evaluated by the function db5val.

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(:)	::	r	`(nr)` array of $r$ abcissae. must be strictly increasing.
integer(kind=ip),	intent(in)		::	nr	number of $r$ abcissae ( $\ge 3$ )
real(kind=wp),	intent(in),	dimension(:,:,:,:,:)	::	fcn	`(nx,ny,nz,nq,nr)` matrix of function values to interpolate. `fcn(i,j,k,q,r)` should contain the function value at the point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
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)		::	kr	the order of spline pieces in $r$ ( $2 \le k_r < n_r$ ). (order = polynomial degree + 1)
integer(kind=ip),	intent(in)		::	iknot	knot sequence flag: 0 = knot sequence chosen by db5ink. 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 db5ink. 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 db5ink. 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 db5ink. 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 db5ink. If `iknot=1` these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(inout),	dimension(:)	::	tr	The `(nr+kr)` knots in the $r$ direction for the spline interpolant. If `iknot=0` these are chosen by db5ink. If `iknot=1` these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(out),	dimension(:,:,:,:,:)	::	bcoef	`(nx,ny,nz,nq,nr)` 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. 19 = `nr` out of range. 20 = `kr` out of range. 21 = `r` not strictly increasing. 22 = `tr` 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)` 704 = `size(r)` $\ne$ `size(fcn,5)` 706 = `size(x)` $\ne$ `nx` 707 = `size(y)` $\ne$ `ny` 708 = `size(z)` $\ne$ `nz` 709 = `size(q)` $\ne$ `nq` 710 = `size(r)` $\ne$ `nr` 712 = `size(tx)` $\ne$ `nx+kx` 713 = `size(ty)` $\ne$ `ny+ky` 714 = `size(tz)` $\ne$ `nz+kz` 715 = `size(tq)` $\ne$ `nq+kq` 716 = `size(tr)` $\ne$ `nr+kr` 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)` 804 = `size(r)` $\ne$ `size(bcoef,5)`

Calls

Help

Called by

Help

Source Code

    pure subroutine db5ink_default(x,nx,y,ny,z,nz,q,nq,r,nr,&
                           fcn,&
                           kx,ky,kz,kq,kr,&
                           iknot,&
                           tx,ty,tz,tq,tr,&
                           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)                         :: nr    !! number of \(r\) 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)
    integer(ip),intent(in)                         :: kr    !! the order of spline pieces in \(r\)
                                                            !! ( \( 2 \le k_r < n_r \) ).
                                                            !! (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)               :: r     !! `(nr)` array of \(r\) abcissae. must be strictly increasing.
    real(wp),dimension(:,:,:,:,:),intent(in)       :: fcn   !! `(nx,ny,nz,nq,nr)` matrix of function values to interpolate.
                                                            !! `fcn(i,j,k,q,r)` should contain the function value at the
                                                            !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
    integer(ip),intent(in)                         :: iknot !! knot sequence flag:
                                                            !!
                                                            !! * 0 = knot sequence chosen by [[db5ink]].
                                                            !! * 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 [[db5ink]].
                                                            !! * 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 [[db5ink]].
                                                            !! * 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 [[db5ink]].
                                                            !! * 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 [[db5ink]].
                                                            !! * If `iknot=1` these are specified by the user.
                                                            !!
                                                            !! Must be non-decreasing.
    real(wp),dimension(:),intent(inout)            :: tr    !! The `(nr+kr)` knots in the \(r\) direction for the spline
                                                            !! interpolant.
                                                            !!
                                                            !! * If `iknot=0` these are chosen by [[db5ink]].
                                                            !! * If `iknot=1` these are specified by the user.
                                                            !!
                                                            !! Must be non-decreasing.
    real(wp),dimension(:,:,:,:,:),intent(out)      :: bcoef !! `(nx,ny,nz,nq,nr)` 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.
                                                            !! * 19 = `nr` out of range.
                                                            !! * 20 = `kr` out of range.
                                                            !! * 21 = `r` not strictly increasing.
                                                            !! * 22 = `tr` 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)`
                                                            !! * 704 = `size(r)`  \( \ne \) `size(fcn,5)`
                                                            !! * 706 = `size(x)`  \( \ne \) `nx`
                                                            !! * 707 = `size(y)`  \( \ne \) `ny`
                                                            !! * 708 = `size(z)`  \( \ne \) `nz`
                                                            !! * 709 = `size(q)`  \( \ne \) `nq`
                                                            !! * 710 = `size(r)`  \( \ne \) `nr`
                                                            !! * 712 = `size(tx)` \( \ne \) `nx+kx`
                                                            !! * 713 = `size(ty)` \( \ne \) `ny+ky`
                                                            !! * 714 = `size(tz)` \( \ne \) `nz+kz`
                                                            !! * 715 = `size(tq)` \( \ne \) `nq+kq`
                                                            !! * 716 = `size(tr)` \( \ne \) `nr+kr`
                                                            !! * 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)`
                                                            !! * 804 = `size(r)`  \( \ne \) `size(bcoef,5)`

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

    !  check validity of input
    call check_inputs(  iknot,&
                        iflag,&
                        nx=nx,ny=ny,nz=nz,nq=nq,nr=nr,&
                        kx=kx,ky=ky,kz=kz,kq=kq,kr=kr,&
                        x=x,y=y,z=z,q=q,r=r,&
                        tx=tx,ty=ty,tz=tz,tq=tq,tr=tr,&
                        f5=fcn,&
                        bcoef5=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)
            call dbknot(r,nr,kr,tr)
        end if

        allocate(temp(nx*ny*nz*nq*nr))
        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),2_ip*kr*(nr+1_ip))))

        !  construct b-spline coefficients
        !  (odd number of passes, so start in bcoef to end up there)
                         call dbtpcf(x,nx,fcn,   nx,ny*nz*nq*nr,tx,kx,bcoef,work,iflag)
        if (iflag==0_ip) call dbtpcf(y,ny,bcoef, ny,nx*nz*nq*nr,ty,ky,temp, work,iflag)
        if (iflag==0_ip) call dbtpcf(z,nz,temp,  nz,nx*ny*nq*nr,tz,kz,bcoef,work,iflag)
        if (iflag==0_ip) call dbtpcf(q,nq,bcoef, nq,nx*ny*nz*nr,tq,kq,temp, work,iflag)
        if (iflag==0_ip) call dbtpcf(r,nr,temp,  nr,nx*ny*nz*nq,tr,kr,bcoef,work,iflag)

        deallocate(temp)
        deallocate(work)

    end if

    end subroutine db5ink_default

db5ink_default Subroutine

Contents

Source Code

private pure subroutine db5ink_default(x, nx, y, ny, z, nz, q, nq, r, nr, fcn, kx, ky, kz, kq, kr, iknot, tx, ty, tz, tq, tr, bcoef, iflag)

History

Arguments

Calls

Called by

Source Code