defc and dfc procedures and support routines from [SLATEC](https: For fitting Bsplines polynomials to discrete 1D data.
For a description of the Bsplines and usage instructions to evaluate them, see:
For further discussion of (constrained) curve fitting using Bsplines, see reference 2.
Note
This module does not support the userdefined ip
integer kind.
It only uses the default integer kind.
Todo
add iflag
outputs to be consistent with the rest of the library.
Type  Visibility  Attributes  Name  Initial  

real(kind=wp),  private,  parameter  ::  drelpr  =  epsilon(1.0_wp) 
machine precision ( 
To test independence of incoming column.
Type  Intent  Optional  Attributes  Name  

integer  ::  me  
integer  ::  mend  
integer  ::  ir  
real(kind=wp)  ::  factor  
real(kind=wp)  ::  tau  
real(kind=wp)  ::  scale(*)  
real(kind=wp)  ::  wic(*) 
dcv is a companion function subprogram for dfc. The documentation for dfc has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(in)  ::  xval 
The point where the variance is desired 

integer,  intent(in)  ::  ndata 
The number of discrete (X,Y) pairs for which dfc calculated a piecewise polynomial curve. 

integer,  intent(in)  ::  nconst 
The number of conditions that constrained the Bspline in dfc. 

integer,  intent(in)  ::  nord 
The order of the Bspline used in dfc. The value of NORD must satisfy 1 < NORD < 20 . 

integer,  intent(in)  ::  nbkpt 
The number of knots in the array BKPT(). The value of NBKPT must satisfy NBKPT .GE. 2NORD. 

real(kind=wp),  intent(in)  ::  bkpt(*) 
The array of knots. Normally the problem data interval will be included between the limits BKPT(NORD) and BKPT(NBKPTNORD+1). The additional end knots BKPT(I),I=1,...,NORD1 and I=NBKPTNORD+2,...,NBKPT, are required by dfc to compute the functions used to fit the data. 

real(kind=wp)  ::  w(*) 
Real work array as used in dfc. See dfc for the required length of W(). The contents of W() must not be modified by the user if the variance function is desired. 
This subprogram fits a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as Bsplines. The fitting is done in a weighted least squares sense.
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 
Bspline 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
Bspline 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 

integer,  intent(in)  ::  Lw 
The amount of working storage actually
allocated for the working array 

real(kind=wp)  ::  w(*) 
Working Array.
Its length is specified as an input parameter
in 
This is a companion subprogram to DEFC. This subprogram does weighted least squares fitting of data by Bspline curves. The documentation for DEFC has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

integer  ::  Ndata  
real(kind=wp)  ::  Xdata(*)  
real(kind=wp)  ::  Ydata(*)  
real(kind=wp)  ::  Sddata(*)  
integer  ::  Nord  
integer  ::  Nbkpt  
real(kind=wp)  ::  Bkptin(*)  
integer  ::  Mdein  
integer  ::  Mdeout  
real(kind=wp)  ::  Coeff(*)  
real(kind=wp)  ::  Bf(Nord,*)  
real(kind=wp)  ::  Xtemp(*)  
real(kind=wp)  ::  Ptemp(*)  
real(kind=wp)  ::  Bkpt(*)  
real(kind=wp)  ::  g(Mdg,*)  
integer  ::  Mdg  
real(kind=wp)  ::  w(Mdw,*)  
integer  ::  Mdw  
integer  ::  Lw 
These subroutines solve the least squares problem Ax = b
for
banded matrices A using sequential accumulation of rows of the
data matrix. Exactly one righthand side vector is permitted.
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(inout)  ::  g(Mdg,*) 


integer,  intent(in)  ::  Mdg 
The number of rows in the working array


integer,  intent(in)  ::  Nb 
The bandwidth of the data matrix 

integer,  intent(inout)  ::  Ip 
Input Set by the user to the value 1 before the first call to DBNDAC. Its subsequent value is controlled by DBNDAC to set up for the next call to DBNDAC. 

integer,  intent(inout)  ::  Ir 
Input
Index of the row of 

integer,  intent(in)  ::  Mt 
Set by the user to indicate the number of new rows of data in the block 

integer,  intent(in)  ::  Jt 
Set by the user to indicate
the index of the first nonzero column in that
set of rows 
These subroutines solve the least squares problem Ax = b
for
banded matrices A using sequential accumulation of rows of the
data matrix. Exactly one righthand side vector is permitted.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  Mode 
Set by the user to one of the values 1, 2, or
3. These values respectively indicate that
the solution of 

real(kind=wp),  intent(in)  ::  g(Mdg,*) 


integer,  intent(in)  ::  Mdg 
The number of rows in the working array


integer,  intent(in)  ::  Nb 
This argument has the same meaning and contents as following the last call to DBNDAC. 

integer,  intent(in)  ::  Ip 
This argument has the same meaning and contents as following the last call to DBNDAC. 

integer,  intent(in)  ::  Ir 
This argument has the same meaning and contents as following the last call to DBNDAC. 

real(kind=wp),  intent(inout)  ::  x(*) 


integer,  intent(in)  ::  n 
The number of variables in the solution
vector. If any of the 

real(kind=wp),  intent(out)  ::  Rnorm 
If 
Calculates the value of all possibly nonzero Bsplines at X
of
order MAX(JHIGH,(J+1)(INDEX1))
on T
.
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(in)  ::  t(*)  
integer,  intent(in)  ::  Jhigh  
integer,  intent(in)  ::  Index  
real(kind=wp),  intent(in)  ::  x  
integer,  intent(in)  ::  Ileft  
real(kind=wp)  ::  Vnikx(*)  
integer,  intent(inout)  ::  j 
JW : added 

real(kind=wp),  intent(inout),  dimension(20)  ::  deltam 
JW : added 

real(kind=wp),  intent(inout),  dimension(20)  ::  deltap 
JW : added 
Construction and/or application of a single
Householder transformation. Q = I + U*(U**T)/B
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  Mode 
1 or 2 to select algorithm H1 or H2 . 

integer,  intent(in)  ::  Lpivot 
the index of the pivot element. 

integer,  intent(in)  ::  l1 
If 

integer,  intent(in)  ::  m 
see 

real(kind=wp),  intent(inout)  ::  u(Iue,*) 
On entry to H1 

integer,  intent(in)  ::  Iue 
the storage increment between elements of 

real(kind=wp),  intent(inout)  ::  Up 
see 

real(kind=wp),  intent(inout)  ::  c(*) 
On entry to H1 or H2 

integer,  intent(in)  ::  Ice 
Storage increment between elements of vectors in 

integer,  intent(in)  ::  Icv 
Storage increment between vectors in 

integer,  intent(in)  ::  Ncv 
Number of vectors in 
Sort an array and optionally make the same interchanges in an auxiliary array. The array may be sorted in increasing or decreasing order.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  n 
number of values in array DX to be sorted 

integer,  intent(in)  ::  Kflag 
control parameter: * Kflag < 0 : sort DX in decreasing order and optionally carry DY along. * Kflag > 0 : sort DX in increasing order and optionally carry DY along. 

real(kind=wp),  intent(inout),  dimension(*)  ::  Dx 
array of values to be sorted (usually abscissas) 

real(kind=wp),  intent(inout),  optional,  dimension(*)  ::  Dy 
array to be (optionally) carried along 
Recursive quicksoft. Modified to also carry along a second array.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  n  
real(kind=wp),  intent(inout),  dimension(*)  ::  dx 
array of values to be sorted 

real(kind=wp),  intent(inout),  optional,  dimension(*)  ::  dy 
array to be (optionally) carried along 
This subprogram fits a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as Bsplines. The fitting is done in a weighted least squares sense. Equality and inequality constraints can be imposed on the fitted curve.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  ndata 
number of points (size of 

real(kind=wp),  intent(in)  ::  xdata(*) 
X data array. No sorting of 

real(kind=wp),  intent(in)  ::  ydata(*) 
Y data array. 

real(kind=wp),  intent(in)  ::  sddata(*) 
Y value standard deviation or uncertainty.
A zero value for any entry of


integer,  intent(in)  ::  nord 
Bspline 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
Bspline package convention. For example,


integer,  intent(in)  ::  Nbkpt 
The value of 

real(kind=wp),  intent(in),  dimension(*)  ::  Bkpt 


integer,  intent(in)  ::  nconst 
The number of conditions that constrain the Bspline is NCONST. A constraint is specified by an (X,Y) pair in the arrays XCONST() and YCONST(), and by the type of constraint and derivative value encoded in the array NDERIV(*). 

real(kind=wp),  intent(in)  ::  xconst(*) 
X value of constraint. No sorting of XCONST(*) is required. 

real(kind=wp),  intent(in)  ::  yconst(*) 
Y value of constraint 

integer,  intent(in)  ::  nderiv(*) 
The value of NDERIV(*) is determined as follows. Suppose the Ith constraint applies to the Jth derivative of the Bspline. (Any nonnegative value of J < NORD is permitted. In particular the value J=0 refers to the Bspline itself.) For this Ith constraint, set 

integer,  intent(inout)  ::  mode 
Input 

real(kind=wp),  intent(out)  ::  coeff(*) 
If the output value of MODE=0 or 1, this array contains the unknowns obtained from the least squares fitting process. These N=NBKPTNORD parameters are the Bspline coefficients. For MODE=1, the equality constraints are contradictory. To make the fitting process more robust, the equality constraints are satisfied in a least squares sense. In this case the array COEFF() contains Bspline coefficients for this extended concept of a solution. If MODE=1,2 or 3 on output, the array COEFF() is undefined. 

real(kind=wp)  ::  w(*) 
real work array of length 

integer  ::  iw(*) 
integer work array of length 
This is a companion subprogram to DFC. The documentation for DFC has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

integer  ::  ndata  
real(kind=wp)  ::  xdata(*)  
real(kind=wp)  ::  ydata(*)  
real(kind=wp)  ::  sddata(*)  
integer  ::  nord  
integer  ::  nbkpt  
real(kind=wp)  ::  bkptin(*)  
integer  ::  nconst  
real(kind=wp)  ::  xconst(*)  
real(kind=wp)  ::  yconst(*)  
integer  ::  nderiv(*)  
integer  ::  mode  
real(kind=wp)  ::  coeff(*)  
real(kind=wp)  ::  bf(nord,*)  
real(kind=wp)  ::  xtemp(*)  
real(kind=wp)  ::  ptemp(*)  
real(kind=wp)  ::  bkpt(*)  
real(kind=wp)  ::  g(mdg,*)  
integer  ::  mdg  
real(kind=wp)  ::  w(mdw,*)  
integer  ::  mdw  
real(kind=wp)  ::  work(*)  
integer  ::  iwork(*) 
Calculates value and derivs of all Bsplines which do not vanish at X
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  t(*)  
integer  ::  k  
real(kind=wp)  ::  x  
integer  ::  ileft  
real(kind=wp)  ::  vnikx(k,*)  
integer  ::  nderiv 
Solve a least squares problem for banded matrices using sequential accumulation of rows of the data matrix. Exactly one righthand side vector is permitted.
Type  Intent  Optional  Attributes  Name  

real(kind=wp),  intent(inout)  ::  a(mda,*) 


integer,  intent(in)  ::  mda 
actual leading dimension of 

integer,  intent(in)  ::  m  
integer,  intent(in)  ::  n  
real(kind=wp),  intent(inout)  ::  b(mdb,*) 


integer,  intent(in)  ::  mdb 
actual leading dimension of 

integer,  intent(in)  ::  nb  
real(kind=wp),  intent(in)  ::  tau 
Absolute tolerance parameter provided by user for pseudorank determination. 

integer,  intent(out)  ::  krank 
Set by the subroutine to indicate the pseudorank of A. 

real(kind=wp),  intent(out)  ::  rnorm(*) 


real(kind=wp)  ::  h(*) 


real(kind=wp)  ::  g(*) 


integer  ::  ip(*) 

Determine an N1vector W, and an N2vector Z which minimizes the Euclidean length of W subject to GW+HZ >= Y. This is the least projected distance problem, LPDP. The matrices G and H are of respective dimensions M by N1 and M by N2.
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  a(mda,*) 


integer,  intent(in)  ::  mda  
integer  ::  m  
integer,  intent(in)  ::  n1  
integer,  intent(in)  ::  n2  
real(kind=wp)  ::  prgopt(*)  
real(kind=wp)  ::  x(*) 


real(kind=wp)  ::  wnorm  
integer,  intent(out)  ::  mode 
The value of MODE indicates the status of the computation after returning to the user. 

real(kind=wp)  ::  ws(*) 


integer  ::  is(*) 

This subprogram solves a linearly constrained least squares problem with both equality and inequality constraints, and, if the user requests, obtains a covariance matrix of the solution parameters.
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  w(mdw,*)  
integer,  intent(in)  ::  mdw  
integer  ::  me  
integer  ::  ma  
integer  ::  mg  
integer  ::  n  
real(kind=wp)  ::  prgopt(*)  
real(kind=wp)  ::  x(*)  
real(kind=wp)  ::  rnorme  
real(kind=wp)  ::  rnorml  
integer  ::  mode  
real(kind=wp)  ::  ws(*)  
integer  ::  ip(3) 
This is a companion subprogram to DLSEI. The documentation for DLSEI has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  w(mdw,*) 


integer,  intent(in)  ::  mdw 
contain (resp) var. dimension of 

integer,  intent(in)  ::  ma 
contain (resp) var. dimension of 

integer,  intent(in)  ::  mg 
contain (resp) var. dimension of 

integer,  intent(in)  ::  n 
contain (resp) var. dimension of 

real(kind=wp),  intent(in)  ::  prgopt(*) 
Program option vector. 

real(kind=wp),  intent(out)  ::  x(*) 
Solution vector(unless MODE=2) 

real(kind=wp),  intent(out)  ::  rnorm 
length of AXB. 

integer,  intent(out)  ::  mode 


real(kind=wp)  ::  ws(*) 
Working storage of dimension 

integer  ::  ip(*) 

This is a companion subprogram to DWNNLS. The documentation for DWNNLS has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  w(mdw,*)  
integer  ::  mdw  
integer  ::  m  
integer  ::  n  
integer  ::  l  
integer  ::  ipivot(*)  
integer  ::  itype(*)  
real(kind=wp)  ::  h(*)  
real(kind=wp)  ::  scale(*)  
real(kind=wp)  ::  rnorm  
integer  ::  idope(*)  
real(kind=wp)  ::  dope(*)  
logical  ::  done 
This is a companion subprogram to DWNNLS. The documentation for DWNNLS has complete usage instructions.
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  w(mdw,*)  
integer  ::  mdw  
integer  ::  mme  
integer  ::  ma  
integer  ::  n  
integer  ::  l  
real(kind=wp)  ::  prgopt(*)  
real(kind=wp)  ::  x(*)  
real(kind=wp)  ::  rnorm  
integer  ::  mode  
integer  ::  ipivot(*)  
integer  ::  itype(*)  
real(kind=wp)  ::  wd(*)  
real(kind=wp)  ::  h(*)  
real(kind=wp)  ::  scale(*)  
real(kind=wp)  ::  z(*)  
real(kind=wp)  ::  temp(*)  
real(kind=wp)  ::  d(*) 
To update the column Sum Of Squares and find the pivot column. The column Sum of Squares Vector will be updated at each step. When numerically necessary, these values will be recomputed.
Type  Intent  Optional  Attributes  Name  

integer  ::  i  
integer  ::  lend  
integer  ::  mend  
integer  ::  ir  
integer  ::  mdw  
logical  ::  recalc  
integer  ::  imax  
real(kind=wp)  ::  hbar  
real(kind=wp)  ::  h(*)  
real(kind=wp)  ::  scale(*)  
real(kind=wp)  ::  w(mdw,*) 
Perform column interchange. Exchange elements of permuted index vector and perform column interchanges.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  i  
integer,  intent(in)  ::  imax  
integer,  intent(in)  ::  m  
integer,  intent(in)  ::  mdw  
integer,  intent(inout)  ::  ipivot(*)  
real(kind=wp),  intent(inout)  ::  h(*)  
real(kind=wp),  intent(inout)  ::  w(mdw,*) 
This subprogram solves a linearly constrained least squares
problem. Suppose there are given matrices E
and A
of
respective dimensions ME
by N
and MA
by N
, and vectors F
and B
of respective lengths ME
and MA
. This subroutine
solves the problem
Type  Intent  Optional  Attributes  Name  

real(kind=wp)  ::  w(mdw,*)  
integer  ::  mdw  
integer  ::  me  
integer  ::  ma  
integer  ::  n  
integer  ::  l  
real(kind=wp)  ::  prgopt(*)  
real(kind=wp)  ::  x(*)  
real(kind=wp)  ::  rnorm  
integer  ::  mode  
integer  ::  iwork(*)  
real(kind=wp)  ::  work(*) 