Description

Multidimensional (1D-6D) B-Spline interpolation of data on a regular grid. Basic subroutine interface.

Notes

This module is based on the bspline and spline routines from [1]. The original Fortran 77 routines were converted to free-form source. Some of them are relatively unchanged from the originals, but some have been extensively refactored. In addition, new routines for 1d, 4d, 5d, and 6d interpolation were also created (these are simply extensions of the same algorithm into higher dimensions).

References

DBSPLIN and DTENSBS from the NIST Core Math Library. Original code is public domain.
Carl de Boor, “A Practical Guide to Splines”, Springer-Verlag, New York, 1978.
Carl de Boor, Efficient Computer Manipulation of Tensor Products, ACM Transactions on Mathematical Software, Vol. 5 (1979), p. 173-182.
D.E. Amos, “Computation with Splines and B-Splines”, SAND78-1968, Sandia Laboratories, March, 1979.
Carl de Boor, Package for calculating with B-splines, SIAM Journal on Numerical Analysis 14, 3 (June 1977), p. 441-472.

Uses

- iso_fortran_env
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.

Used by

Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.
Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.

Variables

Type	Visibility	Attributes		Name		Initial
integer,	private,	parameter	::	wp	=	real64	Real precision

Functions

private function dbvalu(t, a, n, k, ideriv, x, inbv, work, iflag)

Evaluates the b-representation (t,a,n,k) of a b-spline at x for the function value on ideriv=0 or any of its derivatives on ideriv=1,2,…,k-1. right limiting values (right derivatives) are returned except at the right end point x=t(n+1) where left limiting values are computed. the spline is defined on t(k) <= x <= t(n+1). dbvalu returns a fatal error message when x is outside of this interval.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(:)	::	t	knot vector of length n+k
real(kind=wp),	intent(in),	dimension(n)	::	a	b-spline coefficient vector of length n
integer,	intent(in)		::	n	number of b-spline coefficients. (sum of knot multiplicities-k)
integer,	intent(in)		::	k	order of the b-spline, k >= 1
integer,	intent(in)		::	ideriv	order of the derivative, 0 <= ideriv <= k-1. ideriv = 0 returns the b-spline value
real(kind=wp),	intent(in)		::	x	argument, t(k) <= x <= t(n+1)
integer,	intent(inout)		::	inbv	an initialization parameter which must be set to 1 the first time dbvalu is called. inbv contains information for efficient process- ing after the initial call and inbv must not be changed by the user. distinct splines require distinct inbv parameters.
real(kind=wp),		dimension(:)	::	work	work vector of length 3*k
integer,	intent(out)		::	iflag	if 0: no errors if 401: k does not satisfy k>=1 if 402: n does not satisfy n>=k if 403: ideriv does not satisfy 0<=ideriv<k if 404: x is not greater than or equal to t(k) if 405: x is not less than or equal to t(n+1) if 406: a left limiting value cannot be obtained at t(k)

Return Value real(kind=wp)

Subroutines

private subroutine check_inputs(routine, iflag, nx, ny, nz, nq, nr, ns, kx, ky, kz, kq, kr, ks, x, y, z, q, r, s, tx, ty, tz, tq, tr, ts, status_ok)

Check the validity of the inputs to the “ink” routines. Prints warning message if there is an error, and also sets iflag and status_ok.

Arguments

Type	Intent	Optional	Attributes		Name
character(len=*),	intent(in)			::	routine
integer,	intent(inout)			::	iflag
integer,	intent(in),	optional		::	nx
integer,	intent(in),	optional		::	ny
integer,	intent(in),	optional		::	nz
integer,	intent(in),	optional		::	nq
integer,	intent(in),	optional		::	nr
integer,	intent(in),	optional		::	ns
integer,	intent(in),	optional		::	kx
integer,	intent(in),	optional		::	ky
integer,	intent(in),	optional		::	kz
integer,	intent(in),	optional		::	kq
integer,	intent(in),	optional		::	kr
integer,	intent(in),	optional		::	ks
real(kind=wp),	intent(in),	optional,	dimension(:)	::	x
real(kind=wp),	intent(in),	optional,	dimension(:)	::	y
real(kind=wp),	intent(in),	optional,	dimension(:)	::	z
real(kind=wp),	intent(in),	optional,	dimension(:)	::	q
real(kind=wp),	intent(in),	optional,	dimension(:)	::	r
real(kind=wp),	intent(in),	optional,	dimension(:)	::	s
real(kind=wp),	intent(in),	optional,	dimension(:)	::	tx
real(kind=wp),	intent(in),	optional,	dimension(:)	::	ty
real(kind=wp),	intent(in),	optional,	dimension(:)	::	tz
real(kind=wp),	intent(in),	optional,	dimension(:)	::	tq
real(kind=wp),	intent(in),	optional,	dimension(:)	::	tr
real(kind=wp),	intent(in),	optional,	dimension(:)	::	ts
logical,	intent(out)			::	status_ok

public subroutine db1ink(x, nx, fcn, kx, tx, bcoef, iflag)

Determines the parameters of a function that interpolates the one-dimensional gridded data $[x(i),\mathrm{fcn}(i)] ~\mathrm{for}~ i=1,..,n_x$ The interpolating function and its derivatives may subsequently be evaluated by the function db1val.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(nx)	::	x	Array of x abcissae. Must be strictly increasing.
integer,	intent(in)		::	nx	Number of x abcissae
real(kind=wp),	intent(in),	dimension(nx)	::	fcn	Array of function values to interpolate. fcn(i) should contain the function value at the point x(i)
integer,	intent(in)		::	kx	The order of spline pieces in x (>= 2, < nx). (order = polynomial degree + 1)
real(kind=wp),	intent(inout),	dimension(nx+kx)	::	tx	The knots in the x direction for the spline interpolant. If iflag=0 these are chosen by db1ink. If iflag=1 these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(out),	dimension(nx)	::	bcoef	Array of coefficients of the b-spline interpolant.
integer,	intent(inout)		::	iflag	on input: 0 = knot sequence chosen by db1ink. 1 = knot sequence chosen by user. on output: 1 = successful execution. 2 = iflag out of range. 3 = nx out of range. 4 = kx out of range. 5 = x not strictly increasing. 6 = tx not non-decreasing.

public subroutine db1val(xval, idx, tx, nx, kx, bcoef, f, iflag, inbvx)

Evaluates the tensor product piecewise polynomial interpolant constructed by the routine db1ink or one of its derivatives at the point xval.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	xval	x coordinate of evaluation point.
integer,	intent(in)		::	idx	x derivative of piecewise polynomial to evaluate.
real(kind=wp),	intent(in),	dimension(nx+kx)	::	tx	sequence of knots defining the piecewise polynomial in the x direction. (same as in last call to db1ink)
integer,	intent(in)		::	nx	the number of interpolation points in x. (same as in last call to db1ink)
integer,	intent(in)		::	kx	order of polynomial pieces in x. (same as in last call to db1ink)
real(kind=wp),	intent(in),	dimension(nx)	::	bcoef	the b-spline coefficients computed by db1ink.
real(kind=wp),	intent(out)		::	f	interpolated value
integer,	intent(out)		::	iflag	status flag: 0 : no errors, /=0 : error
integer,	intent(inout)		::	inbvx	initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.

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

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

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(nx)	::	x	Array of x abcissae. Must be strictly increasing.
integer,	intent(in)		::	nx	Number of x abcissae
real(kind=wp),	intent(in),	dimension(ny)	::	y	Array of y abcissae. Must be strictly increasing.
integer,	intent(in)		::	ny	Number of y abcissae
real(kind=wp),	intent(in),	dimension(nx,ny)	::	fcn	Array of function values to interpolate. fcn(i,j) should contain the function value at the point (x(i),y(j))
integer,	intent(in)		::	kx	The order of spline pieces in x (>= 2, < nx). (order = polynomial degree + 1)
integer,	intent(in)		::	ky	The order of spline pieces in y (>= 2, < ny). (order = polynomial degree + 1)
real(kind=wp),	intent(inout),	dimension(nx+kx)	::	tx	The knots in the x direction for the spline interpolant. If iflag=0 these are chosen by db2ink. If iflag=1 these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(inout),	dimension(ny+ky)	::	ty	The knots in the y direction for the spline interpolant. If iflag=0 these are chosen by db2ink. If iflag=1 these are specified by the user. Must be non-decreasing.
real(kind=wp),	intent(out),	dimension(nx,ny)	::	bcoef	Array of coefficients of the b-spline interpolant.
integer,	intent(inout)		::	iflag	on input: 0 = knot sequence chosen by db2ink. 1 = knot sequence chosen by user. on output: 1 = successful execution. 2 = iflag 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.

public subroutine db2val(xval, yval, idx, idy, tx, ty, nx, ny, kx, ky, bcoef, f, iflag, inbvx, inbvy, iloy)

Evaluates the tensor product piecewise polynomial interpolant constructed by the routine db2ink or one of its derivatives at the point (xval,yval).

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in)		::	xval	x coordinate of evaluation point.
real(kind=wp),	intent(in)		::	yval	y coordinate of evaluation point.
integer,	intent(in)		::	idx	x derivative of piecewise polynomial to evaluate.
integer,	intent(in)		::	idy	y derivative of piecewise polynomial to evaluate.
real(kind=wp),	intent(in),	dimension(nx+kx)	::	tx	sequence of knots defining the piecewise polynomial in the x direction. (same as in last call to db2ink)
real(kind=wp),	intent(in),	dimension(ny+ky)	::	ty	sequence of knots defining the piecewise polynomial in the y direction. (same as in last call to db2ink)
integer,	intent(in)		::	nx	the number of interpolation points in x. (same as in last call to db2ink)
integer,	intent(in)		::	ny	the number of interpolation points in y. (same as in last call to db2ink)
integer,	intent(in)		::	kx	order of polynomial pieces in x. (same as in last call to db2ink)
integer,	intent(in)		::	ky	order of polynomial pieces in y. (same as in last call to db2ink)
real(kind=wp),	intent(in),	dimension(nx,ny)	::	bcoef	the b-spline coefficients computed by db2ink.
real(kind=wp),	intent(out)		::	f	interpolated value
integer,	intent(out)		::	iflag	status flag: 0 : no errors, /=0 : error
integer,	intent(inout)		::	inbvx	initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer,	intent(inout)		::	inbvy	initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.
integer,	intent(inout)		::	iloy	initialization parameter which must be set to 1 the first time this routine is called, and must not be changed by the user.

private subroutine dbintk(x, y, t, n, k, bcoef, q, work, iflag)

dbintk produces the b-spline coefficients, bcoef, of the b-spline of order k with knots t(i), i=1,…,n+k, which takes on the value y(i) at x(i), i=1,…,n. the spline or any of its derivatives can be evaluated by calls to dbvalu.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(n)	::	x	vector of length n containing data point abscissa in strictly increasing order.
real(kind=wp),	intent(in),	dimension(n)	::	y	corresponding vector of length n containing data point ordinates.
real(kind=wp),	intent(in),	dimension(*)	::	t	knot vector of length n+k since t(1),..,t(k) <= x(1) and t(n+1),..,t(n+k) Read more…
integer,	intent(in)		::	n	number of data points, n >= k
integer,	intent(in)		::	k	order of the spline, k >= 1
real(kind=wp),	intent(out),	dimension(n)	::	bcoef	a vector of length n containing the b-spline coefficients
real(kind=wp),	intent(out),	dimension(*)	::	q	a work vector of length (2k-1)n, containing the triangular factorization of the coefficient matrix of the linear system being solved. the coefficients for the interpolant of an additional data set (x(i),yy(i)), i=1,…,n with the same abscissa can be obtained by loading yy into bcoef and then executing call dbnslv(q,2k-1,n,k-1,k-1,bcoef)
real(kind=wp),	intent(out),	dimension(*)	::	work	work vector of length 2*k
integer,	intent(out)		::	iflag	if 0: no errors. if 100: k does not satisfy k>=1. if 101: n does not satisfy n>=k. if 102: x(i) does not satisfy x(i)<x(i+1) for some i. if 103: some abscissa was not in the support of the. corresponding basis function and the system is singular. if 104: the system of solver detects a singular system. although the theoretical conditions for a solution were satisfied.

private subroutine dbknot(x, n, k, t)

dbknot chooses a knot sequence for interpolation of order k at the data points x(i), i=1,..,n. the n+k knots are placed in the array t. k knots are placed at each endpoint and not-a-knot end conditions are used. the remaining knots are placed at data points if n is even and between data points if n is odd. the rightmost knot is shifted slightly to the right to insure proper interpolation at x(n) (see page 350 of the reference).

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(n)	::	x
integer,	intent(in)		::	n
integer,	intent(in)		::	k
real(kind=wp),	intent(out),	dimension(:)	::	t

private subroutine dbnfac(w, nroww, nrow, nbandl, nbandu, iflag)

Returns in w the LU-factorization (without pivoting) of the banded matrix a of order nrow with (nbandl + 1 + nbandu) bands or diagonals in the work array w .

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(inout),	dimension(nroww,nrow)	::	w	work array. See header for details.
integer,	intent(in)		::	nroww	row dimension of the work array w. must be >= nbandl + 1 + nbandu.
integer,	intent(in)		::	nrow	matrix order
integer,	intent(in)		::	nbandl	number of bands of a below the main diagonal
integer,	intent(in)		::	nbandu	number of bands of a above the main diagonal
integer,	intent(out)		::	iflag	indicating success(=1) or failure (=2)

private subroutine dbnslv(w, nroww, nrow, nbandl, nbandu, b)

Companion routine to dbnfac. it returns the solution x of the linear system a*x = b in place of b, given the lu-factorization for a in the work array w from dbnfac.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(nroww,nrow)	::	w	describes the lu-factorization of a banded matrix a of order nrow as constructed in dbnfac.
integer,	intent(in)		::	nroww	describes the lu-factorization of a banded matrix a of order nrow as constructed in dbnfac.
integer,	intent(in)		::	nrow	describes the lu-factorization of a banded matrix a of order nrow as constructed in dbnfac.
integer,	intent(in)		::	nbandl	describes the lu-factorization of a banded matrix a of order nrow as constructed in dbnfac.
integer,	intent(in)		::	nbandu	describes the lu-factorization of a banded matrix a of order nrow as constructed in dbnfac.
real(kind=wp),	intent(inout),	dimension(nrow)	::	b	in: right side of the system to be solved out: the solution x, of order nrow

private subroutine dbspvn(t, jhigh, k, index, x, ileft, vnikx, work, iwork, iflag)

Calculates the value of all (possibly) nonzero basis functions at x of order max(jhigh,(j+1)*(index-1)), where t(k) <= x <= t(n+1) and j=iwork is set inside the routine on the first call when index=1. ileft is such that t(ileft) <= x < t(ileft+1). a call to dintrv(t,n+1,x,ilo,ileft,mflag) produces the proper ileft. dbspvn calculates using the basic algorithm needed in dbspvd. if only basis functions are desired, setting jhigh=k and index=1 can be faster than calling dbspvd, but extra coding is required for derivatives (index=2) and dbspvd is set up for this purpose.

Arguments

Type	Intent		Name
real(kind=wp),	intent(in)	::	t(*)	knot vector of length n+k, where n = number of b-spline basis functions n = sum of knot multiplicities-k dimension t(ileft+jhigh)
integer,	intent(in)	::	jhigh	order of b-spline, 1 <= jhigh <= k
integer,	intent(in)	::	k	highest possible order
integer,	intent(in)	::	index	index = 1 gives basis functions of order jhigh = 2 denotes previous entry with work, iwork values saved for subsequent calls to dbspvn.
real(kind=wp),	intent(in)	::	x	argument of basis functions, t(k) <= x <= t(n+1)
integer,	intent(in)	::	ileft	largest integer such that t(ileft) <= x < t(ileft+1)
real(kind=wp),	intent(out)	::	vnikx(k)	vector of length k for spline values.
real(kind=wp),	intent(out)	::	work(*)	a work vector of length 2*k
integer,	intent(out)	::	iwork	a work parameter. both work and iwork contain information necessary to continue for index = 2. when index = 1 exclusively, these are scratch variables and can be used for other purposes.
integer,	intent(out)	::	iflag	if 0: no errors if 201: k does not satisfy k>=1 if 202: jhigh does not satisfy 1<=jhigh<=k if 203: index is not 1 or 2 if 204: x does not satisfy t(ileft)<=x<=t(ileft+1)

private subroutine dbtpcf(x, n, fcn, ldf, nf, t, k, bcoef, work, iflag)

dbtpcf computes b-spline interpolation coefficients for nf sets of data stored in the columns of the array fcn. the b-spline coefficients are stored in the rows of bcoef however. each interpolation is based on the n abcissa stored in the array x, and the n+k knots stored in the array t. the order of each interpolation is k. the work array must be of length at least 2k(n+1).

Arguments

Type	Intent		Name
real(kind=wp)		::	x(n)
integer,	intent(in)	::	n
real(kind=wp)		::	fcn(ldf,nf)
integer,	intent(in)	::	ldf
integer,	intent(in)	::	nf
real(kind=wp)		::	t(*)
integer,	intent(in)	::	k
real(kind=wp)		::	bcoef(nf,n)
real(kind=wp)		::	work(*)
integer,	intent(out)	::	iflag	if 0: no errors if 301: n should be >0

private subroutine dintrv(XT, lxt, x, ilo, ileft, mflag)

Computes the largest integer ileft in 1 <= ileft <= lxt such that XT(ileft) <= x where XT(*) is a subdivision of the x interval. precisely,

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(lxt)	::	XT	a knot or break point vector of length lxt
integer,	intent(in)		::	lxt	length of the XT vector
real(kind=wp),	intent(in)		::	x	argument
integer,	intent(inout)		::	ilo	an initialization parameter which must be set to 1 the first time the spline array XT is processed by dintrv. ilo contains information for efficient processing after the initial call and ilo must not be changed by the user. distinct splines require distinct ilo parameters.
integer,	intent(out)		::	ileft	largest integer satisfying XT(ileft) <= x
integer,	intent(out)		::	mflag	signals when x lies out of bounds

bspline Module

Contents

Variables

Functions

Subroutines

Description

Notes

See also

References

Uses

Used by

Variables

Functions

private function dbvalu(t, a, n, k, ideriv, x, inbv, work, iflag)

Arguments

Return Value real(kind=wp)

Subroutines

private subroutine check_inputs(routine, iflag, nx, ny, nz, nq, nr, ns, kx, ky, kz, kq, kr, ks, x, y, z, q, r, s, tx, ty, tz, tq, tr, ts, status_ok)

Arguments

public subroutine db1ink(x, nx, fcn, kx, tx, bcoef, iflag)

Arguments

public subroutine db1val(xval, idx, tx, nx, kx, bcoef, f, iflag, inbvx)

Arguments

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

Arguments

public subroutine db2val(xval, yval, idx, idy, tx, ty, nx, ny, kx, ky, bcoef, f, iflag, inbvx, inbvy, iloy)

Arguments

private subroutine dbintk(x, y, t, n, k, bcoef, q, work, iflag)

Arguments

private subroutine dbknot(x, n, k, t)

Arguments

private subroutine dbnfac(w, nroww, nrow, nbandl, nbandu, iflag)

Arguments

private subroutine dbnslv(w, nroww, nrow, nbandl, nbandu, b)

Arguments

private subroutine dbspvn(t, jhigh, k, index, x, ileft, vnikx, work, iwork, iflag)

Arguments

private subroutine dbtpcf(x, n, fcn, ldf, nf, t, k, bcoef, work, iflag)

Arguments

private subroutine dintrv(XT, lxt, x, ilo, ileft, mflag)

Arguments