bspline Module

  1. mod_bspline.f90
  2. bspline

Determines the parameters of a function that interpolates History

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

Evaluates the tensor product piecewise polynomial History

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

Determines the parameters of a function that interpolates 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.

Evaluates the tensor product piecewise polynomial 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.

Check the validity of the inputs to the "ink" routines. Notes

The code is new, but the logic is based on the original logic in the CMLIB routines db2ink and db3ink.

History

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

dbknot chooses a knot sequence for interpolation of order k at the History

  • Jacob Williams, 2/24/2015 : Refactored this routine.

dbtpcf computes b-spline interpolation coefficients for nf sets History

  • Jacob Williams, 2/24/2015 : Refactored this routine.

dbintk produces the b-spline coefficients, bcoef, of the Error conditions

  • improper input
  • singular system of equations

History

  • splint written by carl de boor [5]
  • dbintk author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • 000330 modified array declarations. (jec)
  • Jacob Williams, 5/10/2015 : converted to free-form Fortran.

Returns in w the LU-factorization (without pivoting) of the banded Work array

Input

    w array of size (nroww,nrow) contains the interesting
    part of a banded matrix  a , with the diagonals or bands of  a
    stored in the rows of  w , while columns of  a  correspond to
    columns of  w . this is the storage mode used in  linpack  and
    results in efficient innermost loops.
       explicitly,  a  has  nbandl  bands below the diagonal
                        +     1     (main) diagonal
                        +   nbandu  bands above the diagonal
    and thus, with    middle = nbandu + 1,
      a(i+j,j)  is in  w(i+middle,j)  for i=-nbandu,...,nbandl
                                          j=1,...,nrow .
    for example, the interesting entries of a (1,2)-banded matrix
    of order  9  would appear in the first  1+1+2 = 4  rows of  w
    as follows.
                      13 24 35 46 57 68 79
                   12 23 34 45 56 67 78 89
                11 22 33 44 55 66 77 88 99
                21 32 43 54 65 76 87 98

    all other entries of  w  not identified in this way with an en-
    try of  a  are never referenced .

Output

  • if iflag = 1, then w contains the lu-factorization of a into a unit lower triangu- lar matrix l and an upper triangular matrix u (both banded) and stored in customary fashion over the corresponding entries of a . this makes it possible to solve any particular linear system a*x = b for x by a call dbnslv ( w, nroww, nrow, nbandl, nbandu, b ) with the solution x contained in b on return .
  • if iflag = 2, then one of nrow-1, nbandl,nbandu failed to be nonnegative, or else one of the potential pivots was found to be zero indicating that a does not have an lu-factorization. this implies that a is singular in case it is totally positive .

History

  • banfac written by carl de boor [5]
  • dbnfac from CMLIB [1]
  • Jacob Williams, 5/10/2015 : converted to free-form Fortran.

Companion routine to dbnfac. it returns the solution x of the History

  • banslv written by carl de boor [5]
  • dbnslv from SLATEC library [1]
  • Jacob Williams, 5/10/2015 : converted to free-form Fortran.

Calculates the value of all (possibly) nonzero basis Error Conditions

  • improper input

History

  • bsplvn written by carl de boor [5]
  • dbspvn author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • 000330 modified array declarations. (jec)
  • Jacob Williams, 2/24/2015 : extensive refactoring of CMLIB routine.

Evaluates the b-representation (t,a,n,k) of a b-spline Error Conditions

  • improper input

History

  • bvalue written by carl de boor [5]
  • dbvalu author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • 000330 modified array declarations. (jec)
  • Jacob Williams, 2/24/2015 : extensive refactoring of CMLIB routine.

Computes the largest integer ileft in 1 <= ileft <= lxt History

  • interv written by carl de boor [5]
  • dintrv author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • Jacob Williams, 2/24/2015 : updated to free-form Fortran.

Uses

  • module~~bspline~~UsesGraph module~bspline bspline iso_fortran_env iso_fortran_env module~bspline->iso_fortran_env

Used by

  • module~~bspline~~UsedByGraph module~bspline bspline module~ms_film ms_film module~ms_film->module~bspline module~test_musst test_musst module~test_musst->module~ms_film module~inout_files inout_files module~test_musst->module~inout_files module~inout_files->module~ms_film program~main main program~main->module~test_musst

Contents

Variables

wp

Functions

Subroutines


Variables

TypeVisibility AttributesNameInitial
integer, private, parameter:: wp =real64

Real precision

Functions

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

Arguments

Type IntentOptional AttributesName
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

Return Value real(kind=wp)


Subroutines

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

Arguments

Type IntentOptional AttributesName
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)

Arguments

Type IntentOptional AttributesName
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)

Arguments

Type IntentOptional AttributesName
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)

Arguments

Type IntentOptional AttributesName
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 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

Type IntentOptional AttributesName
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

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

Arguments

Type IntentOptional AttributesName
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 dbtpcf(x, n, fcn, ldf, nf, t, k, bcoef, work, iflag)

Arguments

Type IntentOptional AttributesName
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

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

Arguments

Type IntentOptional AttributesName
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

corresponding basis function and the system is singular. 104: the system of solver detects a singular system. although the theoretical conditions for a solution were satisfied.

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

Arguments

Type IntentOptional AttributesName
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)

Arguments

Type IntentOptional AttributesName
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)

Arguments

Type IntentOptional AttributesName
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

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

Arguments

Type IntentOptional AttributesName
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