N-dimensional cubic spline coefficient calculation by weighted least squares on arbitrarily located data.
The spline (or its derivatives) may then be evaluated by using function SPLFE (or SPLDE).
A grid of evenly spaced nodes in NDIM space is defined by the arguments XMIN, XMAX and NODES. A linear basis for the class of natural splines on these nodes is formed, and a set of corresponding coefficients is computed in the array COEF. These coefficients are chosen to minimize the weighted sum of squared errors between the spline and the arbitrarily located data values described by the arguments XDATA, YDATA and NDATA. The smoothness of the spline in data sparse areas is controlled by the argument XTRAP.
In order to understand the arguments of this routine, one should realize that the node grid need not bear any particular relation to the data points. In the theory of exact-fit interpolatory splines, the nodes would in fact be data locations, but in this case they serve only to define the class of splines from which the approximating function is chosen. This node grid is a rectangular arrangement of points in NDIM space, with the restriction that along any coordinate direction the nodes are equally spaced. The class of natural splines on this grid of nodes (NDIM-cubic splines whose 2nd derivatives normal to the boundaries are 0) has as many degrees of freedom as the grid has nodes. Thus the smoothness or flexibility of the splines is determined by the choice of the node grid.
An overdetermined system of linear equations is formed -- one equation for each data point plus equations for derivative constraints. This system is solved using subroutine suprls.
If there is exactly one data point in the near vicinity of each node and no extra data, the resulting spline will agree with the data values to machine accuracy. However, if the problem is overdetermined or the sparse data option is utilized, the accuracy is hard to predict. Basically, smooth functions require fewer nodes than rough ones for the same accuracy.
The execution time is roughly proportional
to NDATA*NCOF**2 where NCOF = NODES(1)*...*NODES(NDIM).
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(splpak_type), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | ndim |
The dimensionality of the problem. The
spline is a function of |
||
| real(kind=wp), | intent(in) | :: | xdata(l1xdat,ndata) |
A collection of locations for the data
values, i.e., points from the independent
variable space. This collection is a
2-dimensional array whose 1st dimension
indexes the |
||
| integer, | intent(in) | :: | l1xdat |
The length of the 1st dimension of Note:For 1-dimensional problems |
||
| real(kind=wp), | intent(in) | :: | ydata(ndata) |
A collection of data values corresponding to
the points in |
||
| real(kind=wp), | intent(in) | :: | wdata(:) |
A collection of weights. Thus, if the reliability of a data point is known to be low, the corresponding weight may be made small (relative to the other weights) so that the sum over all data points is affected less by discrepencies at the unreliable point. Data points with zero weight are completely ignored. Note:If The dimension is assumed to be |
||
| integer, | intent(in) | :: | ndata |
The number of data points mentioned in the above arguments. |
||
| real(kind=wp), | intent(in) | :: | xmin(ndim) |
A vector describing the lower extreme corner
of the node grid. A set of evenly spaced
nodes is formed along each coordinate axis
and |
||
| real(kind=wp), | intent(in) | :: | xmax(ndim) |
A vector describing the upper extreme corner
of the node grid. A set of evenly spaced
nodes is formed along each coordinate axis
and |
||
| integer, | intent(in) | :: | nodes(ndim) |
A vector of integers describing the number of
nodes along each axis. Note:The node grid is completely defined by
the arguments A node in this grid may be indexed by
an |
||
| real(kind=wp), | intent(in) | :: | xtrap |
A parameter to control extrapolation to data
sparse areas. The region described by Note:If |
||
| real(kind=wp), | intent(out) | :: | coef(ncf) |
The array of coefficients computed by this
routine. Each coefficient corresponds to a
particular basis function which in turn
corresponds to a node in the node grid. This
correspondence between the node grid and the
array |
||
| integer, | intent(in) | :: | ncf |
The length of the array |
||
| real(kind=wp) | :: | work(nwrk) |
A workspace array for solving the least
squares matrix generated by this routine.
Its required size is a function of the total
number of nodes in the node grid. This
total, |
|||
| integer, | intent(in) | :: | nwrk |
The length of the array |
||
| integer, | intent(out) | :: | ierror |
An error flag with the following meanings:
|
subroutine splcw(me,ndim,xdata,l1xdat,ydata,wdata,ndata,xmin,xmax, & nodes,xtrap,coef,ncf,work,nwrk,ierror) class(splpak_type),intent(inout) :: me integer,intent(in) :: ndim !! The dimensionality of the problem. The !! spline is a function of `NDIM` variables or !! coordinates and thus a point in the !! independent variable space is an `NDIM` vector. !! `NDIM` must be `>= 1`. integer,intent(in) :: l1xdat !! The length of the 1st dimension of `XDATA` in !! the calling program. `L1XDAT` must be `>= NDIM`. !! !!#### Note: !! For 1-dimensional problems `L1XDAT` is usually 1. integer,intent(in) :: ncf !! The length of the array `COEF` in the calling !! program. If `NCF` is `< NODES(1)*...*NODES(NDIM)`, !! a fatal error is diagnosed. integer,intent(in) :: nwrk !! The length of the array `WORK` in the calling !! program. If !! `NCOL = NODES(1)*...*NODES(NDIM)` is the total !! number of nodes, then a fatal error is !! diagnosed if `NWRK` is less than !! `NCOL*(NCOL+1)`. integer,intent(in) :: ndata !! The number of data points mentioned in the !! above arguments. real(wp),intent(in) :: xdata(l1xdat,ndata) !! A collection of locations for the data !! values, i.e., points from the independent !! variable space. This collection is a !! 2-dimensional array whose 1st dimension !! indexes the `NDIM` coordinates of a given point !! and whose 2nd dimension labels the data !! point. For example, the data point with !! label `IDATA` is located at the point !! `(XDATA(1,IDATA),...,XDATA(NDIM,IDATA))` where !! the elements of this vector are the values of !! the `NDIM` coordinates. The location, number !! and ordering of the data points is arbitrary. !! The dimension of `XDATA` is assumed to be !! `XDATA(L1XDAT,NDATA)`. real(wp),intent(in) :: ydata(ndata) !! A collection of data values corresponding to !! the points in `XDATA`. `YDATA(IDATA)` is the !! data value associated with the point !! `(XDATA(1,IDATA),...,XDATA(NDIM,IDATA))` in the !! independent variable space. The spline whose !! coefficients are computed by this routine !! approximates these data values in the least !! squares sense. The dimension is assumed to be !! `YDATA(NDATA)`. real(wp),intent(in) :: wdata(:) !! A collection of weights. `WDATA(IDATA)` is a !! weight associated with the data point !! labelled `IDATA`. It should be non-negative, !! but may be of any magnitude. The weights !! have the effect of forcing greater or lesser !! accuracy at a given point as follows: this !! routine chooses coefficients to minimize the !! sum over all data points of the quantity !!```fortran !! (WDATA(IDATA)*(YDATA(IDATA) ! spline value at XDATA(IDATA)))**2. !!``` !! Thus, if the reliability !! of a data point is known to be low, the !! corresponding weight may be made small !! (relative to the other weights) so that the !! sum over all data points is affected less by !! discrepencies at the unreliable point. Data !! points with zero weight are completely !! ignored. !! !!#### Note: !! If `WDATA(1)` is `< 0`, the other !! elements of `WDATA` are not !! referenced, and all weights are !! assumed to be unity. !! !! The dimension is assumed to be `WDATA(NDATA)` !! unless `WDATA(1) < 0.`, in which case the !! dimension is assumed to be 1. real(wp),intent(in) :: xmin(ndim) !! A vector describing the lower extreme corner !! of the node grid. A set of evenly spaced !! nodes is formed along each coordinate axis !! and `XMIN(IDIM)` is the location of the first !! node along the `IDIM` axis. The dimension is !! assumed to be `XMIN(NDIM)`. real(wp),intent(in) :: xmax(ndim) !! A vector describing the upper extreme corner !! of the node grid. A set of evenly spaced !! nodes is formed along each coordinate axis !! and `XMAX(IDIM)` is the location of the last !! node along the `IDIM` axis. The dimension is !! assumed to be `XMAX(NDIM)`. real(wp),intent(in) :: xtrap !! A parameter to control extrapolation to data !! sparse areas. The region described by `XMIN` !! and `XMAX` is divided into rectangles, the !! number of which is determined by `NODES`, and !! any rectangle containing a disproportionately !! small number of data points is considered to !! be data sparse (rectangle is used here to !! mean `NDIM`-dimensional rectangle). If `XTRAP` !! is nonzero the least squares problem is !! augmented with derivative constraints in the !! data sparse areas to prevent the matrix from !! becoming poorly conditioned. `XTRAP` serves as !! a weight for these constraints, and thus may !! be used to control smoothness in data sparse !! areas. Experience indicates that unity is a !! good first guess for this parameter. !! !!#### Note: !! If `XTRAP` is zero, substantial !! portions of the routine will be !! skipped, but a singular matrix !! can result if large portions of !! the region are without data. integer,intent(in) :: nodes(ndim) !! A vector of integers describing the number of !! nodes along each axis. `NODES(IDIM)` is the !! number of nodes (counting endpoints) along !! the `IDIM` axis and determines the flexibility !! of the spline in that coordinate direction. !! `NODES(IDIM)` must be `>= 4`, but may be as !! large as the arrays `COEF` and `WORK` allow. !! The dimension is assumed to be `NODES(NDIM)`. !! !!#### Note: !! The node grid is completely defined by !! the arguments `XMIN`, `XMAX` and `NODES`. !! The spacing of this grid in the `IDIM` !! coordinate direction is: !!```fortran !! DX(IDIM) = (XMAX(IDIM)-XMIN(IDIM)) / (NODES(IDIM)-1). !!``` !! A node in this grid may be indexed by !! an `NDIM` vector of integers !! `(IN(1),...,IN(NDIM))` where !! `1 <= IN(IDIM) <= NODES(IDIM)`. !! The location of such a node may be !! represented by an `NDIM` vector !! `(X(1),...,X(NDIM))` where !! `X(IDIM) = XMIN(IDIM) + (IN(IDIM)-1) * DX(IDIM)`. real(wp) :: work(nwrk) !! A workspace array for solving the least !! squares matrix generated by this routine. !! Its required size is a function of the total !! number of nodes in the node grid. This !! total, `NCOL = NODES(1)*...*NODES(NDIM)`, is !! also the number of columns in the least !! squares matrix. The length of the array `WORK` !! must equal or exceed `NCOL*(NCOL+1)`. real(wp),intent(out) :: coef(ncf) !! The array of coefficients computed by this !! routine. Each coefficient corresponds to a !! particular basis function which in turn !! corresponds to a node in the node grid. This !! correspondence between the node grid and the !! array `COEF` is as if `COEF` were an !! `NDIM`-dimensional Fortran array with !! dimensions `NODES(1),...,NODES(NDIM)`, i.e., to !! store the array linearly, the leftmost !! indices are incremented most frequently. !! Hence the length of the `COEF` array must equal !! or exceed the total number of nodes, which is !! `NODES(1)*...*NODES(NDIM)`. The computed array !! `COEF` may be used with function [[SPLFE]] !! (or [[SPLDE]]) to evaluate the spline (or its !! derivatives) at an arbitrary point in `NDIM` !! space. The dimension is assumed to be `COEF(NCF)`. integer,intent(out) :: ierror !! An error flag with the following meanings: !! !! * ` 0` No error. !! * `101` `NDIM` is < 1. !! * `102` `NODES(IDIM)` is < 4 for some `IDIM`. !! * `103` `XMIN(IDIM) = XMAX(IDIM)` for some `IDIM`. !! * `104` `NCF` (size of `COEF`) is `< NODES(1)*...*NODES(NDIM)`. !! * `105` `NDATA` is `< 1`. !! * `106` `NWRK` (size of `WORK`) is too small. !! * `107` [[suprls]] failure (usually insufficient !! data) -- ordinarily occurs only if !! `XTRAP` is zero or `WDATA` contains all !! zeros. real(wp),dimension(:),allocatable :: x integer,dimension(:),allocatable :: nderiv,in,inmx real(wp) :: xrng,swght,rowwt,rhs,basm,reserr,totlwt,& bump,wtprrc,expect,dcwght integer :: ncol,idim,nod,nwrk1,mdata,nwlft,irow,idata,& icol,it,lserr,iin,nrect,idimc,idm,jdm,inidim logical :: boundary real(wp),parameter :: spcrit = 0.75_wp !! SPCRIT is used to determine data sparseness as follows - !! the weights assigned to all data points are totaled into the !! variable TOTLWT. (If no weights are entered, it is set to !! NDATA.) Each node of the node network is assigned a !! rectangle (in which it is contained) and the weights of all !! data points which fall in that rectangle are totaled. If that !! total is less than SPCRIT*EXPECT (EXPECT is defined below), !! then the node is ascertained to be in a data sparse location. !! EXPECT is that fraction of TOTLWT that would be expected by !! comparing the area of the rectangle with the total area under !! consideration. ! size the arrays: call me%destroy(ndim) allocate(x(ndim)) allocate(nderiv(ndim)) allocate(in(ndim)) allocate(inmx(ndim)) ierror = 0 me%mdim = ndim if (me%mdim<1) then ierror = 101 call cfaerr(ierror, & ' splcc or splcw - NDIM is less than 1') return end if ncol = 1 do idim = 1,me%mdim nod = nodes(idim) if (nod<4) then ierror = 102 call cfaerr(ierror, & ' splcc or splcw - NODES(IDIM) is less than 4 for some IDIM') return end if ! Number of columns in least squares matrix = number of coefficients = ! product of nodes over all dimensions. ncol = ncol*nod xrng = xmax(idim) - xmin(idim) if (xrng==0.0_wp) then ierror = 103 call cfaerr(ierror, & ' splcc or splcw - XMIN(IDIM) equals XMAX(IDIM) for some IDIM') return end if ! DX(IDIM) is the node spacing along the IDIM coordinate. me%dx(idim) = xrng/real(nod-1,wp) me%dxin(idim) = 1.0_wp/me%dx(idim) nderiv(idim) = 0 end do if (ncol>ncf) then ierror = 104 call cfaerr(ierror, & ' splcc or splcw - NCF (size of COEF) is too small') return end if nwrk1 = 1 mdata = ndata if (mdata<1) then ierror = 105 call cfaerr(ierror, & ' splcc or splcw - Ndata Is less than 1') return end if ! SWGHT is a local variable = XTRAP, and can be considered a smoothing ! weight for data sparse areas. If SWGHT == 0, no smoothing ! computations are performed. swght = xtrap ! Set aside workspace for counting data points. if (swght/=0.0_wp) nwrk1 = ncol + 1 ! NWLFT is the length of the remaining workspace. nwlft = nwrk - nwrk1 + 1 if (nwlft<1) then ierror = 106 call cfaerr(ierror, & ' splcc or splcw - NWRK (size of WORK) is too small') return end if irow = 0 ! ROWWT is used to weight rows of the least squares matrix. rowwt = 1.0_wp ! Loop through all data points, computing a row for each. do idata = 1,mdata ! WDATA(1)<0 means weights have not been entered. In that case, ! ROWWT is left equal to 1. for all points. Otherwise ROWWT is ! equal to WDATA(IDATA). ! ! Every element of the row, as well as the corresponding right hand ! side, is multiplied by ROWWT. if (wdata(1)>=0.0_wp) then rowwt = wdata(idata) ! Data points with 0 weight are ignored. if (rowwt==0.0_wp) cycle end if irow = irow + 1 ! One row of the least squares matrix corresponds to each data ! point. The right hand for that row will correspond to the ! function value YDATA at that point. rhs = rowwt*ydata(idata) do idim = 1,me%mdim x(idim) = xdata(idim,idata) end do ! The COEF array serves as a row of least squares matrix. ! Its value is zero except for columns corresponding to functions ! which are nonzero at X. do icol = 1,ncol coef(icol) = 0.0_wp end do ! Compute the indices of basis functions which are nonzero at X. ! IBMN is in the range 0 to nodes-2 and IBMX is in range 1 ! to NODES-1. do idim = 1,me%mdim nod = nodes(idim) it = me%dxin(idim)* (x(idim)-xmin(idim)) me%ibmn(idim) = min(max(it-1,0),nod-2) me%ib(idim) = me%ibmn(idim) me%ibmx(idim) = max(min(it+2,nod-1),1) end do basis_index : do ! Begining of basis index loop - traverse all indices corresponding ! to basis functions which are nonzero at X. The indices are in ! IB and are passed through common to BASCMP. call me%bascmp(x,nderiv,xmin,nodes,icol,basm) ! BASCMP computes ICOL and BASM where BASM is the value at X of ! the N-dimensional basis function corresponding to column ICOL. coef(icol) = rowwt*basm ! Increment the basis indices. do idim = 1,me%mdim me%ib(idim) = me%ib(idim) + 1 if (me%ib(idim)<=me%ibmx(idim)) cycle basis_index me%ib(idim) = me%ibmn(idim) end do exit basis_index ! End of basis index loop. end do basis_index ! Send a row of the least squares matrix to the reduction routine. call me%suprls(irow,coef,ncol,rhs,work(nwrk1),nwlft,coef,reserr,lserr) if (lserr/=0) then ierror = 107 call cfaerr(ierror, ' splcc or splcw - suprls failure '//& '(this usually indicates insufficient input data)') end if end do ! Row computations for all data points are now complete. ! ! If SWGHT==0, the least squares matrix is complete and no ! smoothing rows are computed. if (swght/=0.0_wp) then ! Initialize smoothing computations for data sparse areas. ! Derivative constraints will always have zero right hand side. rhs = 0.0_wp nrect = 1 ! Initialize the node indices and compute number of rectangles ! formed by the node network. do idim = 1,me%mdim in(idim) = 0 inmx(idim) = nodes(idim) - 1 nrect = nrect*inmx(idim) end do ! Every node is assigned an element of the workspace (set aside ! previously) in which data points are counted. do iin = 1,ncol work(iin) = 0.0_wp end do ! Assign each data point to a node, total the assignments for ! each node, and save in the workspace. totlwt = 0.0_wp do idata = 1,mdata ! BUMP is the weight associated with the data point. bump = 1.0_wp if (wdata(1)>=0.0_wp) bump = wdata(idata) if (bump==0.0_wp) cycle ! Find the nearest node. iin = 0 do idimc = 1,me%mdim idim = me%mdim + 1 - idimc inidim = int(me%dxin(idim)* (xdata(idim,idata)-xmin(idim))+0.5_wp) ! Points not in range (+ or - 1/2 node spacing) are not counted. if (inidim<0 .or. inidim>inmx(idim)) cycle ! Compute linear address of node in workspace by Horner's method. iin = (inmx(idim)+1)*iin + inidim end do ! Bump counter for that node. work(iin+1) = work(iin+1) + bump totlwt = totlwt + bump end do ! Compute the expected weight per rectangle. wtprrc = totlwt/real(nrect,wp) ! IN contains indices of the node (previously initialized). ! IIN will be the linear address of the node in the workspace. iin = 0 ! Loop through all nodes, computing derivative constraint rows ! for those in data sparse locations. ! ! Begining of node index loop - traverse all node indices. ! The indices are in IN. node_index : do iin = iin + 1 expect = wtprrc ! Rectangles at edge of network are smaller and hence less weight ! should be expected. do idim = 1,me%mdim if (in(idim)==0 .or. in(idim)==inmx(idim)) expect = 0.5_wp*expect end do ! The expected weight minus the actual weight serves to define ! data sparseness and is also used to weight the derivative ! constraint rows. ! ! There is no constraint if not data sparse. if (work(iin)<spcrit*expect) then dcwght = expect - work(iin) do idim = 1,me%mdim inidim = in(idim) ! Compute the location of the node. x(idim) = xmin(idim) + real(inidim,wp)*me%dx(idim) ! Compute the indices of the basis functions which are non-zero ! at the node. me%ibmn(idim) = inidim - 1 me%ibmx(idim) = inidim + 1 ! Distinguish the boundaries. if (inidim==0) me%ibmn(idim) = 0 if (inidim==inmx(idim)) me%ibmx(idim) = inmx(idim) ! Initialize the basis indices. me%ib(idim) = me%ibmn(idim) end do ! Multiply by the extrapolation parameter (this acts as a ! smoothing weight). dcwght = swght*dcwght ! The COEF array serves as a row of the least squares matrix. ! Its value is zero except for columns corresponding to functions ! which are non-zero at the node. do icol = 1,ncol coef(icol) = 0.0_wp end do ! The 2nd derivative of a function of MDIM variables may be thought ! of as a symmetric MDIM x MDIM matrix of 2nd order partial ! derivatives. Traverse the upper triangle of this matrix and, ! for each element, compute a row of the least squares matrix. do idm = 1,me%mdim do jdm = idm,me%mdim do idim = 1,me%mdim nderiv(idim) = 0 end do boundary = .true. ! Off-diagonal elements appear twice by symmetry, so the corresponding ! row is weighted by a factor of 2. rowwt = 2.0_wp*dcwght if (jdm==idm) then ! Diagonal. rowwt = dcwght nderiv(jdm) = 2 if (in(idm)/=0 .and. in(idm)/=inmx(idm)) then boundary = .false. end if end if if (boundary) then ! Node is at boundary. ! ! Normal 2nd derivative constraint at boundary is not appropriate for ! natural splines (2nd derivative 0 by definition). Substitute ! a 1st derivative constraint. nderiv(idm) = 1 nderiv(jdm) = 1 end if irow = irow + 1 basis : do ! Begining of basis index loop - traverse all indices corresponding ! to basis functions which are non-zero at X. ! The indices are in IB and are passed through common to BASCMP. call me%bascmp(x,nderiv,xmin,nodes,icol,basm) ! BASCMP computes ICOL and BASM where BASM is the value at X of the ! N-dimensional basis function corresponding to column ICOL. coef(icol) = rowwt*basm ! Increment the basis indices. do idim = 1,me%mdim me%ib(idim) = me%ib(idim) + 1 if (me%ib(idim)<=me%ibmx(idim)) cycle basis me%ib(idim) = me%ibmn(idim) end do ! End of basis index loop. exit basis end do basis ! Send row of least squares matrix to reduction routine. call me%suprls(irow,coef,ncol,rhs,work(nwrk1),nwlft,coef,reserr,lserr) if (lserr/=0) then ierror = 107 call cfaerr(ierror, & ' splcc or splcw - suprls failure '//& '(this usually indicates insufficient input data)') end if end do end do end if ! Increment node indices. do idim = 1,me%mdim in(idim) = in(idim) + 1 if (in(idim)<=inmx(idim)) cycle node_index in(idim) = 0 end do exit node_index ! End of node index loop. end do node_index end if ! Call for least squares solution in COEF array. irow = 0 call me%suprls(irow,coef,ncol,rhs,work(nwrk1),nwlft,coef,reserr,lserr) if (lserr/=0) then ierror = 107 call cfaerr(ierror, & ' splcc or splcw - suprls failure '//& '(this usually indicates insufficient input data)') end if end subroutine splcw