bspline_sub_module – bspline-fortran

Description

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

Notes

This module is based on the B-spline 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.
D.E. Amos, "Quadrature subroutines for splines and B-splines", Report SAND79-1825, Sandia Laboratories, December 1979.

Uses

- iso_fortran_env
- bspline_kinds_module
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.

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.

Variables

Type	Visibility	Attributes		Name		Initial
integer(kind=ip),	public,	parameter	::	bspline_order_linear	=	2_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_quadratic	=	3_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_cubic	=	4_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_quartic	=	5_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_quintic	=	6_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_hexic	=	7_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_heptic	=	8_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]
integer(kind=ip),	public,	parameter	::	bspline_order_octic	=	9_ip	spline order `k` parameter (for input to the `db*ink` routines) [order = polynomial degree + 1]

Interfaces

public interface db1ink

1D initialization routines.

private pure subroutine db1ink_default(x, nx, fcn, kx, iknot, 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.

History

Jacob Williams, 10/30/2015 : Created 1D 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
real(kind=wp),	intent(in),	dimension(:)	::	fcn	`(nx)` array of function values to interpolate. `fcn(i)` should contain the function value at the point `x(i)`
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)		::	iknot	knot sequence flag: Read more…
real(kind=wp),	intent(inout),	dimension(:)	::	tx	The `(nx+kx)` knots in the $x$ direction for the spline interpolant: Read more…
real(kind=wp),	intent(out),	dimension(:)	::	bcoef	`(nx)` array of coefficients of the b-spline interpolant.
integer(kind=ip),	intent(out)		::	iflag	status flag: Read more…

private pure subroutine db1ink_alt(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, kntopt, tx, bcoef, iflag)

Alternate version of db1ink_default, where the boundary conditions can be specified.

History

Jacob Williams, 9/4/2018 : created this routine.

private pure subroutine db1ink_alt_2(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, tleft, tright, tx, bcoef, iflag)

Alternate version of db1ink_alt, where the first and last 3 knots are specified by the user.

History

Jacob Williams, 9/4/2018 : created this routine.

public interface db1val

1D evaluation routines.

private pure subroutine db1val_default(xval, idx, tx, nx, kx, bcoef, f, iflag, inbvx, w0, extrap)

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

To evaluate the interpolant itself, set idx=0, to evaluate the first partial with respect to x, set idx=1, and so on.

db1val returns 0.0 if (xval,yval) is out of range. that is, if

   xval < tx(1) .or. xval > tx(nx+kx)

if the knots tx were chosen by db1ink, then this is equivalent to:

   xval < x(1) .or. xval > x(nx)+epsx

where

   epsx = 0.1*(x(nx)-x(nx-1))

The input quantities tx, nx, kx, and bcoef should be unchanged since the last call of db1ink.

History

Jacob Williams, 10/30/2015 : Created 1D routine.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	xval	$x$ coordinate of evaluation point.
integer(kind=ip),	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(kind=ip),	intent(in)			::	nx	the number of interpolation points in $x$ . (same as in last call to db1ink)
integer(kind=ip),	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(kind=ip),	intent(out)			::	iflag	status flag: Read more…
integer(kind=ip),	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.
real(kind=wp),	intent(inout),		dimension(3_ip*kx)	::	w0	work array
logical,	intent(in),	optional		::	extrap	if extrapolation is allowed (if not present, default is False)

private pure subroutine db1val_alt(xval, idx, tx, nx, n, kx, bcoef, f, iflag, inbvx, w0, extrap)

Alternate version of db1val_default for use with db1ink_alt and db1ink_alt_2.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	xval	$x$ coordinate of evaluation point.
integer(kind=ip),	intent(in)			::	idx	$x$ derivative of piecewise polynomial to evaluate.
real(kind=wp),	intent(in),		dimension(n+kx)	::	tx	sequence of knots defining the piecewise polynomial in the $x$ direction.
integer(kind=ip),	intent(in)			::	nx	the number of interpolation points in $x$ .
integer(kind=ip),	intent(in)			::	n	length of `bcoef`: `nx+2`
integer(kind=ip),	intent(in)			::	kx	order of polynomial pieces in $x$ . (same as in last call to db1ink)
real(kind=wp),	intent(in),		dimension(n)	::	bcoef	the b-spline coefficients computed by db1ink.
real(kind=wp),	intent(out)			::	f	interpolated value
integer(kind=ip),	intent(out)			::	iflag	status flag: Read more…
integer(kind=ip),	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.
real(kind=wp),	intent(inout),		dimension(3_ip*kx)	::	w0	work array
logical,	intent(in),	optional		::	extrap	if extrapolation is allowed (if not present, default is False)

Abstract Interfaces

abstract interface

public function b1fqad_func(x) result(f)

interface for the input function in dbfqad

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	x

Return Value real(kind=wp)

f(x)

Functions

private pure function check_value(x, t, i, extrap) result(iflag)

Checks if the value is withing the range of the knot vectors. This is called by the various db*val routines.

Arguments

Type	Intent	Optional	Attributes		Name
real(kind=wp),	intent(in)			::	x	the value to check
real(kind=wp),	intent(in),		dimension(:)	::	t	the knot vector
integer(kind=ip),	intent(in)			::	i	1=x, 2=y, 3=z, 4=q, 5=r, 6=s
logical,	intent(in),	optional		::	extrap	if extrapolation is allowed (if not present, default is False)

Return Value integer(kind=ip)

returns 0 if value is OK, otherwise returns 600+i

private pure function get_temp_x_for_extrap(x, tmin, tmax, extrap) result(xt)

Returns the value of x to use for computing the interval in t, depending on if extrapolation is allowed or not.

Type	Intent	Optional		Name
real(kind=wp),	intent(in)		::	x	variable value
real(kind=wp),	intent(in)		::	tmin	first knot vector element for b-splines
real(kind=wp),	intent(in)		::	tmax	last knot vector element for b-splines
logical,	intent(in),	optional	::	extrap	if extrapolation is allowed (if not present, default is False)

Type	Intent	Attributes		Name
procedure(b1fqad_func)			::	fun	external function of one argument for the integrand `bf(x)=fun(x)*dbvalu(tx,bcoef,nx,kx,id,x,inbv,work)`
real(kind=wp),	intent(in),	dimension(nx+kx)	::	tx	knot array
real(kind=wp),	intent(in),	dimension(nx)	::	bcoef	b-spline coefficient array
integer(kind=ip),	intent(in)		::	nx	length of coefficient array
integer(kind=ip),	intent(in)		::	kx	order of b-spline, `kx >= 1`
integer(kind=ip),	intent(in)		::	idx	order of the spline derivative, `0 <= idx <= k-1` `idx=0` gives the spline function
real(kind=wp),	intent(in)		::	x1	left point of quadrature interval in `t(k) <= x <= t(n+1)`
real(kind=wp),	intent(in)		::	x2	right point of quadrature interval in `t(k) <= x <= t(n+1)`
real(kind=wp),	intent(in)		::	tol	desired accuracy for the quadrature, suggest `10*dtol < tol <= 0.1` where `dtol` is the maximum of `1.0e-300` and real(wp) unit roundoff for the machine
real(kind=wp),	intent(out)		::	f	integral of `bf(x)` on `(x1,x2)`
integer(kind=ip),	intent(out)		::	iflag	status flag: Read more…
real(kind=wp),	intent(inout),	dimension(3_ip*kx)	::	w0	work array for dbfqad

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

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(*)	::	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(kind=ip),	intent(in)		::	jhigh	order of b-spline, `1 <= jhigh <= k`
integer(kind=ip),	intent(in)		::	k	highest possible order
integer(kind=ip),	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(kind=ip),	intent(in)		::	ileft	largest integer such that `t(ileft) <= x < t(ileft+1)`
real(kind=wp),	intent(out),	dimension(k)	::	vnikx	vector of length `k` for spline values.
real(kind=wp),	intent(inout),	dimension(*)	::	work	a work vector of length `2*k`
integer(kind=ip),	intent(inout)		::	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(kind=ip),	intent(out)		::	iflag	Read more…

Type	Intent	Attributes		Name
procedure(b1fqad_func)			::	fun	name of external function of one argument which multiplies dbvalu.
real(kind=wp),	intent(in),	dimension(:)	::	xt	knot array for dbvalu
real(kind=wp),	intent(in),	dimension(n)	::	bc	b-coefficient array for dbvalu
integer(kind=ip),	intent(in)		::	n	number of b-coefficients for dbvalu
integer(kind=ip),	intent(in)		::	kk	order of the spline, `kk>=1`
integer(kind=ip),	intent(in)		::	id	Order of the spline derivative, `0<=id<=kk-1`
real(kind=wp),	intent(in)		::	a	lower limit of integral
real(kind=wp),	intent(in)		::	b	upper limit of integral (may be less than `a`)
integer(kind=ip),	intent(inout)		::	inbv	initialization parameter for dbvalu
real(kind=wp),	intent(inout)		::	err	IN: is a requested pseudorelative error tolerance. normally pick a value of `abs(err)<1e-3`. `ans` will normally have no more error than `abs(err)` times the integral of the absolute value of `fun(x)*[[dbvalu]]()`. Read more…
real(kind=wp),	intent(out)		::	ans	computed value of integral
integer(kind=ip),	intent(out)		::	iflag	a status code: Read more…
real(kind=wp),	intent(inout),	dimension(:)	::	work	work vector of length `3*k` for dbvalu

bspline_sub_module Module

Contents

Variables

Interfaces

Abstract Interfaces

Functions

Subroutines

Description

Notes

See also

References

Uses

Used by

Variables

Interfaces

public interface db1ink

private pure subroutine db1ink_default(x, nx, fcn, kx, iknot, tx, bcoef, iflag)

History

Arguments

private pure subroutine db1ink_alt(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, kntopt, tx, bcoef, iflag)

History

See also

Arguments

private pure subroutine db1ink_alt_2(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, tleft, tright, tx, bcoef, iflag)

History

See also

Arguments

public interface db1val

private pure subroutine db1val_default(xval, idx, tx, nx, kx, bcoef, f, iflag, inbvx, w0, extrap)

History

Arguments

private pure subroutine db1val_alt(xval, idx, tx, nx, n, kx, bcoef, f, iflag, inbvx, w0, extrap)

Arguments

Abstract Interfaces

abstract interface

public function b1fqad_func(x) result(f)

Arguments

Return Value real(kind=wp)

Functions

private pure function check_value(x, t, i, extrap) result(iflag)

Arguments

Return Value integer(kind=ip)

private pure function get_temp_x_for_extrap(x, tmin, tmax, extrap) result(xt)

Arguments

Return Value real(kind=wp)

public pure function get_status_message(iflag) result(msg)

Arguments

Return Value character(len=:), allocatable

Subroutines

private pure subroutine db1ink_default(x, nx, fcn, kx, iknot, tx, bcoef, iflag)

Arguments

private pure subroutine db1ink_alt(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, kntopt, tx, bcoef, iflag)

Arguments

private pure subroutine db1ink_alt_2(x, nx, fcn, kx, ibcl, ibcr, fbcl, fbcr, tleft, tright, tx, bcoef, iflag)

Arguments

private pure subroutine db1val_default(xval, idx, tx, nx, kx, bcoef, f, iflag, inbvx, w0, extrap)

Arguments

private pure subroutine db1val_alt(xval, idx, tx, nx, n, kx, bcoef, f, iflag, inbvx, w0, extrap)

Arguments

public pure subroutine db1sqad(tx, bcoef, nx, kx, x1, x2, f, iflag, w0)

Arguments

public subroutine db1fqad(fun, tx, bcoef, nx, kx, idx, x1, x2, tol, f, iflag, w0)

Arguments

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

Arguments

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

Arguments

public pure subroutine db3ink(x, nx, y, ny, z, nz, fcn, kx, ky, kz, iknot, tx, ty, tz, bcoef, iflag)

Arguments

public pure subroutine db3val(xval, yval, zval, idx, idy, idz, tx, ty, tz, nx, ny, nz, kx, ky, kz, bcoef, f, iflag, inbvx, inbvy, inbvz, iloy, iloz, w2, w1, w0, extrap)

Arguments

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

Arguments

public pure subroutine db4val(xval, yval, zval, qval, idx, idy, idz, idq, tx, ty, tz, tq, nx, ny, nz, nq, kx, ky, kz, kq, bcoef, f, iflag, inbvx, inbvy, inbvz, inbvq, iloy, iloz, iloq, w3, w2, w1, w0, extrap)

Arguments

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

Arguments

public pure subroutine db5val(xval, yval, zval, qval, rval, idx, idy, idz, idq, idr, tx, ty, tz, tq, tr, nx, ny, nz, nq, nr, kx, ky, kz, kq, kr, bcoef, f, iflag, inbvx, inbvy, inbvz, inbvq, inbvr, iloy, iloz, iloq, ilor, w4, w3, w2, w1, w0, extrap)

Arguments

public pure subroutine db6ink(x, nx, y, ny, z, nz, q, nq, r, nr, s, ns, fcn, kx, ky, kz, kq, kr, ks, iknot, tx, ty, tz, tq, tr, ts, bcoef, iflag)