dpchbs – PCHIP

public subroutine dpchbs(n, x, f, d, incfd, knotyp, nknots, t, bcoef, ndim, kord, ierr)

Piecewise Cubic Hermite to B-Spline converter.

DPCHBS computes the B-spline representation of the PCH function determined by N,X,F,D. To be compatible with the rest of PCHIP, DPCHBS includes INCFD, the increment between successive values of the F- and D-arrays.

The output is the B-representation for the function: NKNOTS, T, BCOEF, NDIM, KORD.

Caution:

Since it is assumed that the input PCH function has been computed by one of the other routines in the package PCHIP, input arguments N, X, INCFD are not checked for validity.

Restrictions/assumptions:

N>=2 . (not checked)
X(i)<X(i+1), i=1,...,N . (not checked)
INCFD>0 . (not checked)
KNOTYP<=2 . (error return if not)
NKNOTS = NDIM+4 = 2*N+4 . (error return if not)
T(2k+1) = T(2k) = X(k), k=1,...,N . (not checked)

5 an 6 apply only if KNOTYP<0 .

Portability:

Argument INCFD is used only to cause the compiler to generate efficient code for the subscript expressions (1+(I-1)*INCFD) . The normal usage, in which DPCHBS is called with one-dimensional arrays F and D, is probably non-Fortran 77, in the strict sense, but it works on all systems on which DPCHBS has been tested.

References

F. N. Fritsch, "Representations for parametric cubic splines," Computer Aided Geometric Design 6 (1989), pp.79-82.

Arguments

Type	Intent		Name
integer,	intent(in)	::	n	the number of data points, N>=2 . (not checked)
real(kind=wp),	intent(in)	::	x(*)	array of independent variable values. The elements of X must be strictly increasing: X(I-1) < X(I), I = 2(1)N. (not checked) nmax, the dimension of X, must be >=N.
real(kind=wp),	intent(in)	::	f(incfd,*)	array of dependent variable values. F(1+(I-1)*INCFD) is the value corresponding to X(I). nmax, the second dimension of F, must be >=N.
real(kind=wp),	intent(in)	::	d(incfd,*)	array of derivative values at the data points. D(1+(I-1)*INCFD) is the value corresponding to X(I). nmax, the second dimension of D, must be >=N.
integer,	intent(in)	::	incfd	the increment between successive values in F and D. This argument is provided primarily for 2-D applications. It may have the value 1 for one-dimensional applications, in which case F and D may be singly-subscripted arrays.
integer,	intent(in)	::	knotyp	a flag to control the knot sequence. The knot sequence T is normally computed from X by putting a double knot at each X and setting the end knot pairs according to the value of KNOTYP: KNOTYP = 0: Quadruple knots at X(1) and X(N). (default) KNOTYP = 1: Replicate lengths of extreme subintervals: T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1)) ; T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)). KNOTYP = 2: Periodic placement of boundary knots: T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1)); T(M+4) = T(M+3) = X(N) + (X(2)-X(1)) . Here M=NDIM=2*N. If the input value of KNOTYP is negative, however, it is assumed that NKNOTS and T were set in a previous call. This option is provided for improved efficiency when used in a parametric setting.
integer,	intent(inout)	::	nknots	the number of knots. If KNOTYP>=0, then NKNOTS will be set to NDIM+4. If KNOTYP<0, then NKNOTS is an input variable, and an error return will be taken if it is not equal to NDIM+4.
real(kind=wp),	intent(inout)	::	t(*)	array of 2N+4 knots for the B-representation. If KNOTYP>=0, T will be returned by DPCHBS with the interior double knots equal to the X-values and the boundary knots set as indicated above. If KNOTYP<0, it is assumed that T was set by a previous call to DPCHBS. (This routine does not* verify that T forms a legitimate knot sequence.)
real(kind=wp),	intent(out)	::	bcoef(*)	array of 2*N B-spline coefficients.
integer,	intent(out)	::	ndim	the dimension of the B-spline space. (Set to 2*N.)
integer,	intent(out)	::	kord	the order of the B-spline. (Set to 4.)
integer,	intent(out)	::	ierr	an error flag: Normal return: IERR = 0 (no errors). "Recoverable" errors: IERR = -4 if KNOTYP>2 . IERR = -5 if KNOTYP<0 and NKNOTS/=(2*N+4).

Calls

Help

Source Code

    subroutine dpchbs (n, x, f, d, incfd, knotyp, nknots, t, bcoef, ndim, kord, ierr)

    integer,intent(in)     :: n               !! the number of data points, N>=2 .  (not checked)
    integer,intent(in)     :: incfd           !! the increment between successive values in F and D.
                                              !! This argument is provided primarily for 2-D applications.
                                              !! It may have the value 1 for one-dimensional applications,
                                              !! in which case F and D may be singly-subscripted arrays.
    integer,intent(in)     :: knotyp          !! a flag to control the knot sequence.
                                              !! The knot sequence T is normally computed from X by putting
                                              !! a double knot at each X and setting the end knot pairs
                                              !! according to the value of KNOTYP:
                                              !!
                                              !!  * KNOTYP = 0:  Quadruple knots at X(1) and X(N).  (default)
                                              !!  * KNOTYP = 1:  Replicate lengths of extreme subintervals:
                                              !!                 T( 1 ) = T( 2 ) = X(1) - (X(2)-X(1))  ;
                                              !!                 T(M+4) = T(M+3) = X(N) + (X(N)-X(N-1)).
                                              !!  * KNOTYP = 2:  Periodic placement of boundary knots:
                                              !!                 T( 1 ) = T( 2 ) = X(1) - (X(N)-X(N-1));
                                              !!                 T(M+4) = T(M+3) = X(N) + (X(2)-X(1))  .
                                              !! Here M=NDIM=2*N.
                                              !!
                                              !! If the input value of KNOTYP is negative, however, it is
                                              !! assumed that NKNOTS and T were set in a previous call.
                                              !! This option is provided for improved efficiency when used
                                              !! in a parametric setting.
    integer,intent(inout)  :: nknots          !! the number of knots.
                                              !!
                                              !!  * If KNOTYP>=0, then NKNOTS will be set to NDIM+4.
                                              !!  * If KNOTYP<0, then NKNOTS is an input variable, and an
                                              !!    error return will be taken if it is not equal to NDIM+4.
    integer,intent(out)    :: ndim            !! the dimension of the B-spline space.  (Set to 2*N.)
    integer,intent(out)    :: kord            !! the order of the B-spline.  (Set to 4.)
    integer,intent(out)    :: ierr            !! an error flag:
                                              !!
                                              !! * Normal return:
                                              !!    * IERR = 0  (no errors).
                                              !! * "Recoverable" errors:
                                              !!    * IERR = -4  if KNOTYP>2 .
                                              !!    * IERR = -5  if KNOTYP<0 and NKNOTS/=(2*N+4).
    real(wp),intent(in)    :: x(*)            !! array of independent variable values.  The
                                              !! elements of X must be strictly increasing:
                                              !!      X(I-1) < X(I),  I = 2(1)N.   (not checked)
                                              !! nmax, the dimension of X, must be >=N.
    real(wp),intent(in)    :: f(incfd,*)      !! array of dependent variable values.
                                              !! F(1+(I-1)*INCFD) is the value corresponding to X(I).
                                              !! nmax, the second dimension of F, must be >=N.
    real(wp),intent(in)    :: d(incfd,*)      !! array of derivative values at the data points.
                                              !! D(1+(I-1)*INCFD) is the value corresponding to X(I).
                                              !! nmax, the second dimension of D, must be >=N.
    real(wp),intent(inout) :: t(*)            !! array of 2*N+4 knots for the B-representation.
                                              !!
                                              !! * If KNOTYP>=0, T will be returned by DPCHBS with the
                                              !!   interior double knots equal to the X-values and the
                                              !!   boundary knots set as indicated above.
                                              !! * If KNOTYP<0, it is assumed that T was set by a
                                              !!   previous call to DPCHBS.  (This routine does **not**
                                              !!   verify that T forms a legitimate knot sequence.)
    real(wp),intent(out)   :: bcoef(*)        !! array of 2*N B-spline coefficients.

    integer :: k, kk
    real(wp) :: dov3, hnew, hold

    ! Initialize.
    ndim = 2*n
    kord = 4
    ierr = 0

    ! check argument validity.  set up knot sequence if ok.
    if ( knotyp>2 ) then
        ierr = -1
        call xermsg ('PCHIP', 'dpchbs', 'knotyp greater than 2', ierr, 1)
        return
    endif
    if ( knotyp<0 ) then
        if ( nknots/=ndim+4 ) then
            ierr = -2
            call xermsg ('PCHIP', 'dpchbs', 'knotyp<0 and nknots/=(2*n+4)', ierr, 1)
            return
        endif
    else
        ! set up knot sequence.
        nknots = ndim + 4
        call dpchkt (n, x, knotyp, t)
    endif

    ! Compute B-spline coefficients.
    hnew = t(3) - t(1)
    do k = 1, n
        kk = 2*k
        hold = hnew
        ! The following requires mixed mode arithmetic.
        dov3 = d(1,k)/3
        bcoef(kk-1) = f(1,k) - hold*dov3
        ! The following assumes T(2*K+1) = X(K).
        hnew = t(kk+3) - t(kk+1)
        bcoef(kk) = f(1,k) + hnew*dov3
    end do

    end subroutine dpchbs

dpchbs Subroutine

Contents

Source Code

public subroutine dpchbs(n, x, f, d, incfd, knotyp, nknots, t, bcoef, ndim, kord, ierr)

Caution:

Restrictions/assumptions:

Portability:

See Also

References

Arguments

Calls

Source Code