This routine does basis function computations for natural
splines. This routine is called by routines SPLCW and SPLDE
to compute ICOL and BASM, which are defined as follows:
MDIM indices in IB (defined in splpak_type) determine
a specific node in the node grid (see routine SPLCC for a
description of the node grid). Every node is associated
with an MDIM-dimensional basis function and a corresponding
column in the least squares matrix (or element of the
coefficient vector). The column index (which may be thought
of as a linear address for the MDIM-dimensional node grid)
corresponding to the specified node is computed as ICOL. The
associated basis function evaluated at X (an MDIM-vector) is
computed as BASM (a scalar).In case NDERIV is not all zero, BASM will not be the value of
the basis function but rather a partial derivative of that
function as follows:
IDIM coordinate is NDERIV(IDIM) (for IDIM <= MDIM). This
routine will compute incorrect values if NDERIV(IDIM) is not
in the range 0 to 2.The technique of this routine is to transform the independent
variable in each dimension such that the nodes fall on
suitably chosen integers. On this transformed space, the
1-dimensional basis functions and their derivatives have a
particularly simple form. The desired MDIM-dimensional basis
function (or any of its partial derivatives) is computed as
a product of such 1-dimensional functions (tensor product
method of defining multi-dimensional splines). The values
which determine the location of the nodes, and hence the
above transform, are passed through common and the argument
list.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(splpak_type), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in) | :: | x(:) | |||
| integer, | intent(in) | :: | nderiv(:) | |||
| real(kind=wp), | intent(in) | :: | xmin(:) | |||
| integer, | intent(in) | :: | nodes(:) | |||
| integer, | intent(out) | :: | icol | |||
| real(kind=wp), | intent(out) | :: | basm |
subroutine bascmp(me,x,nderiv,xmin,nodes,icol,basm) class(splpak_type),intent(inout) :: me real(wp),intent(in) :: x(:) integer,intent(in) :: nderiv(:) real(wp),intent(in) :: xmin(:) integer,intent(in) :: nodes(:) integer,intent(out) :: icol real(wp),intent(out) :: basm real(wp) :: xb,bas1,z,fact,z1 integer :: idim,mdmid,ntyp,ngo ! ICOL will be a linear address corresponding to the indices in IB. icol = 0 ! BASM will be M-dimensional basis function evaluated at X. basm = 1.0_wp do idim = 1,me%mdim ! Compute ICOL by Horner's method. mdmid = me%mdim + 1 - idim icol = nodes(mdmid)*icol + me%ib(mdmid) ! NGO depends upon function type and NDERIV. ntyp = 1 ! Function type 1 (left linear) for IB = 0 or 1. if (me%ib(idim)>1) then ntyp = 2 ! Function type 2 (chapeau function) for 2 LT IB LT NODES-2. if (me%ib(idim)>=nodes(idim)-2) then ntyp = 3 end if end if ! Function type 3 (right linear) for IB = NODES-2 or NODES-1. ngo = 3*ntyp + nderiv(idim) - 2 ! XB is X value of node IB (center of basis function). xb = xmin(idim) + real(me%ib(idim),wp)*me%dx(idim) ! BAS1 will be the 1-dimensional basis function evaluated at X. bas1 = 0.0_wp select case (ngo) case(4) ! Function type 2 (chapeau function). ! ! Transform so that XB is at the origin and the other nodes are at ! the integers. z = abs(me%dxin(idim)* (x(idim)-xb)) - 2.0_wp ! This chapeau function is then that unique cubic spline which is ! identically zero for ABS(Z) GE 2 and is 1 at the origin. This ! function is the general interior node basis function. if (z<0.0_wp) then bas1 = -0.25_wp*z**3 z = z + 1.0_wp if (z<0.0_wp) then bas1 = bas1 + z**3 end if end if case(5) ! 1st derivative. z = x(idim) - xb fact = me%dxin(idim) if (z<0.0_wp) fact = -fact z = fact*z - 2.0_wp if (z<0.0_wp) then bas1 = -0.75_wp*z**2 z = z + 1.0_wp if (z<0.0_wp) then bas1 = bas1 + 3.0_wp*z**2 end if bas1 = fact*bas1 end if case(6) ! 108 ! 2nd derivative. fact = me%dxin(idim) z = fact*abs(x(idim)-xb) - 2.0_wp if (z<0.0_wp) then bas1 = -1.5_wp*z z = z + 1.0_wp if (z<0.0_wp) then bas1 = bas1 + 6.0_wp*z end if bas1 = (fact**2)*bas1 end if case(2,8) ! 1st derivative. if (ngo==2) then fact = -me%dxin(idim) else if (ngo==8) then fact = me%dxin(idim) end if z = fact* (x(idim)-xb) + 2.0_wp if (z<0.0_wp) then if (z<2.0_wp) then bas1 = 1.5_wp*z**2 z = z - 1.0_wp if (z>0.0_wp) then bas1 = bas1 - 3.0_wp*z**2 end if bas1 = fact*bas1 else bas1 = 3.0_wp*fact end if end if case(3,9) ! 2nd derivative. if (ngo==3) then fact = -me%dxin(idim) else if (ngo==9) then fact = me%dxin(idim) end if z = fact* (x(idim)-xb) + 2.0_wp z1 = z - 1.0_wp if (abs(z1)<1.0_wp) then bas1 = 3.0_wp*z if (z1>0.0_wp) then bas1 = bas1 - 6.0_wp*z1 end if bas1 = (fact**2)*bas1 end if case default ! case(1,7) ! or ngo some other value (does that ever happen?) ! (due to the computed goto in the original code) if (ngo/=7) then ! Function type 1 (left linear) is mirror image of function type 3. ! ! Transform so that XB is at 2 and the other nodes are at the integers ! (with ordering reversed to form a mirror image). z = me%dxin(idim)* (xb-x(idim)) + 2.0_wp else ! Function type 3 (right linear). ! ! Transform so that XB is at 2 and the other nodes are at the integers. z = me%dxin(idim)* (x(idim)-xb) + 2.0_wp end if ! This right linear function is defined to be that unique cubic spline ! which is identically zero for Z <= 0 and is 3*Z-3 for Z GE 2. ! This function (obviously having zero 2nd derivative for ! Z GE 2) is used for the two nodes nearest an edge in order ! to generate natural splines, which must by definition have ! zero 2nd derivative at the boundary. ! ! Note that this method of generating natural splines also provides ! a linear extrapolation which has 2nd order continuity with ! the interior splines at the boundary. if (z>0.0_wp) then if (z<2.0_wp) then bas1 = 0.5_wp*z**3 z = z - 1.0_wp if (z>0.0_wp) then bas1 = bas1 - z**3 end if else bas1 = 3.0_wp*z - 3.0_wp end if end if end select basm = basm*bas1 end do icol = icol + 1 end subroutine bascmp