This subprogram fits a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as B-splines. The fitting is done in a weighted least squares sense.
The data can be processed in groups of modest size. The size of the group is chosen by the user. This feature may be necessary for purposes of using constrained curve fitting with subprogram DFC on a very large data set.
To evaluate derivative number IDER at XVAL,
use the function subprogram DBVALU.
f = dbvalu(bkpt,coeff,nbkpt-nord,nord,ider,xval,inbv,workb)
The output of this subprogram will not be
defined unless an output value of MDEOUT=1
was obtained from DEFC, XVAL is in the data
interval, and IDER is nonnegative and < NORD.
The first time DBVALU is called, INBV=1
must be specified. This value of INBV is the
overwritten by DBVALU. The array WORKB(*)
must be of length at least 3*NORD, and must
not be the same as the W(*) array used in the
call to DEFC.
DBVALU expects the breakpoint array BKPT(*)
to be sorted.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | Ndata |
number of points (size of |
||
| real(kind=wp), | intent(in), | dimension(ndata) | :: | Xdata |
X data array. No sorting of |
|
| real(kind=wp), | intent(in), | dimension(ndata) | :: | Ydata |
Y data array. |
|
| real(kind=wp), | intent(in), | dimension(ndata) | :: | Sddata |
Y value standard deviation or uncertainty.
A zero value for any entry of
|
|
| integer, | intent(in) | :: | Nord |
B-spline order.
(The order of the spline is one more than the
degree of the piecewise polynomial defined on
each interval. This is consistent with the
B-spline package convention. For example,
|
||
| integer, | intent(in) | :: | Nbkpt |
The value of |
||
| real(kind=wp), | intent(in), | dimension(:) | :: | Bkpt |
|
|
| integer, | intent(in) | :: | Mdein |
An integer flag, with one of two possible values (1 or 2), that directs the subprogram action with regard to new data points provided by the user:
|
||
| integer, | intent(out) | :: | Mdeout |
An output flag that indicates the status of the curve fit:
|
||
| real(kind=wp), | intent(out) | :: | Coeff(*) |
If the output value of If the user is not satisfied with the fitted
curve returned by DEFC, the constrained
least squares curve fitting subprogram DFC
may be required. The work done within DEFC
to accumulate the data can be utilized by
the user, if so desired. This involves
saving the first |
||
| integer, | intent(in) | :: | Lw |
The amount of working storage actually
allocated for the working array The length of the array |
||
| real(kind=wp) | :: | w(*) |
Working Array.
Its length is specified as an input parameter
in |
subroutine defc(Ndata, Xdata, Ydata, Sddata, Nord, Nbkpt, Bkpt, Mdein, & Mdeout, Coeff, Lw, w) integer,intent(in) :: Ndata !! number of points (size of `xdata` and `ydata`). !! Any non-negative value of `NDATA` is allowed. !! A negative value of `NDATA` is an error. real(wp),dimension(ndata), intent(in) :: Xdata !! X data array. No sorting of `XDATA(*)` is required. real(wp),dimension(ndata), intent(in) :: Ydata !! Y data array. real(wp),dimension(ndata), intent(in) :: Sddata !! Y value standard deviation or uncertainty. !! A zero value for any entry of !! `SDDATA(*)` will weight that data point as 1. !! Otherwise the weight of that data point is !! the reciprocal of this entry. integer,intent(in) :: Nord !! B-spline order. !! (The order of the spline is one more than the !! degree of the piecewise polynomial defined on !! each interval. This is consistent with the !! B-spline package convention. For example, !! `NORD=4` when we are using piecewise cubics.) !! `NORD` must be in the range `1 <= NORD <= 20`. integer,intent(in) :: Nbkpt !! The value of `NBKPT` must satisfy the condition `NBKPT >= 2*NORD`. real(wp),dimension(:),intent(in) :: Bkpt !! `NBKPT` knots of the B-spline. !! Normally the !! problem data interval will be included between !! the limits `BKPT(NORD)` and `BKPT(NBKPT-NORD+1)`. !! The additional end knots `BKPT(I),I=1,...,NORD-1` !! and `I=NBKPT-NORD+2,...,NBKPT`, are !! required to compute the functions used to fit !! the data. No sorting of `BKPT(*)` is required. !! Internal to [[DEFC]] the extreme end knots may !! be reduced and increased respectively to !! accommodate any data values that are exterior !! to the given knot values. The contents of !! `BKPT(*)` is not changed. integer,intent(in) :: Mdein !! An integer flag, with one of two possible !! values (1 or 2), that directs the subprogram !! action with regard to new data points provided !! by the user: !! !! * `= 1` The first time that [[DEFC]] has been !! entered. There are NDATA points to process. !! * `= 2` This is another entry to DEFC(). The !! subprogram [[DEFC]] has been entered with MDEIN=1 !! exactly once before for this problem. There !! are NDATA new additional points to merge and !! process with any previous points. !! (When using [[DEFC]] with MDEIN=2 it is !! important that the set of knots remain fixed at the !! same values for all entries to [[DEFC]].) integer,intent(out) :: Mdeout !! An output flag that indicates the status !! of the curve fit: !! !! * `=-1` A usage error of [[DEFC]] occurred. The !! offending condition is noted with the SLATEC !! library error processor, `XERMSG( )`. In case !! the working array `W(*)` is not long enough, the !! minimal acceptable length is printed. !! !! * `=1` The B-spline coefficients for the fitted !! curve have been returned in array `COEFF(*)`. !! !! * `=2` Not enough data has been processed to !! determine the B-spline coefficients. !! The user has one of two options. Continue !! to process more data until a unique set !! of coefficients is obtained, or use the !! subprogram [[DFC]] to obtain a specific !! set of coefficients. The user should read !! the usage instructions for [[DFC]] for further !! details if this second option is chosen. real(wp),intent(out) :: Coeff(*) !! If the output value of `MDEOUT=1`, this array !! contains the unknowns obtained from the least !! squares fitting process. These `N=NBKPT-NORD` !! parameters are the B-spline coefficients. !! For `MDEOUT=2`, not enough data was processed to !! uniquely determine the B-spline coefficients. !! In this case, and also when `MDEOUT=-1`, all !! values of `COEFF(*)` are set to zero. !! !! If the user is not satisfied with the fitted !! curve returned by [[DEFC]], the constrained !! least squares curve fitting subprogram [[DFC]] !! may be required. The work done within [[DEFC]] !! to accumulate the data can be utilized by !! the user, if so desired. This involves !! saving the first `(NBKPT-NORD+3)*(NORD+1)` !! entries of `W(*)` and providing this data !! to [[DFC]] with the "old problem" designation. !! The user should read the usage instructions !! for subprogram [[DFC]] for further details. integer,intent(in) :: Lw !! The amount of working storage actually !! allocated for the working array `W(*)`. !! This quantity is compared with the !! actual amount of storage needed in [[DEFC]]. !! Insufficient storage allocated for `W(*)` is !! an error. This feature was included in [[DEFC]] !! because misreading the storage formula !! for `W(*)` might very well lead to subtle !! and hard-to-find programming bugs. !! !! The length of the array `W(*)` must satisfy !!``` !! LW >= (NBKPT-NORD+3)*(NORD+1)+ !! (NBKPT+1)*(NORD+1)+ !! 2*MAX(NDATA,NBKPT)+NBKPT+NORD**2 !!``` real(wp) :: w(*) !! Working Array. !! Its length is specified as an input parameter !! in `LW` as noted above. The contents of `W(*)` !! must not be modified by the user between calls !! to [[DEFC]] with values of `MDEIN=1,2,2,...` . !! The first `(NBKPT-NORD+3)*(NORD+1)` entries of !! `W(*)` are acceptable as direct input to [[DFC]] !! for an "old problem" only when `MDEOUT=1` or `2`. integer :: lbf, lbkpt, lg, lptemp, lww, lxtemp, mdg, mdw ! LWW=1 USAGE IN DEFCMN( ) OF W(*).. ! LWW,...,LG-1 W(*,*) ! LG,...,LXTEMP-1 G(*,*) ! LXTEMP,...,LPTEMP-1 XTEMP(*) ! LPTEMP,...,LBKPT-1 PTEMP(*) ! LBKPT,...,LBF BKPT(*) (LOCAL TO DEFCMN( )) ! LBF,...,LBF+NORD**2 BF(*,*) mdg = Nbkpt + 1 mdw = Nbkpt - Nord + 3 lww = 1 lg = lww + mdw*(Nord + 1) lxtemp = lg + mdg*(Nord + 1) lptemp = lxtemp + max(Ndata, Nbkpt) lbkpt = lptemp + max(Ndata, Nbkpt) lbf = lbkpt + Nbkpt call defcmn(Ndata, Xdata, Ydata, Sddata, Nord, Nbkpt, Bkpt, Mdein, Mdeout, & Coeff, w(lbf), w(lxtemp), w(lptemp), w(lbkpt), w(lg), mdg, & w(lww), mdw, Lw) end subroutine defc