/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:43 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "dpfit.h"
#include <float.h>
#include <stdlib.h>
/* PARAMETER translations */
#define FAC 1.01e0
#define HALF .5e0
#define ONE 1.e0
#define TEN 10.e0
#define TWO 2.e0
#define ZERO 0.e0
/* end of PARAMETER translations */
void /*FUNCTION*/ dpfit(
long m,
double x[],
double y[],
double sig[],
long nmax,
LOGICAL32 seekn,
LOGICAL32 comtrn,
LOGICAL32 chbbas,
double p[],
long *nfit,
double *sigfac,
double *w)
{
#define W(I_,J_) (*(w+(I_)*(nmax + 3)+(J_)))
long int _l0, i, idata, idim, ii, irank, irow, j, limit, lrow,
n, np1, np2, np3, nsolve;
double cmin, denom, param[5], s, s2, sigma, size, t,
temp, teneps, xmax, xmin;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const P = &p[0] - 1;
double *const Param = ¶m[0] - 1;
double *const Sig = &sig[0] - 1;
double *const X = &x[0] - 1;
double *const Y = &y[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 1994-10-20 DPFIT Krogh Changes to use M77CON
*>> 1990-08-13 DPFIT CLL Fixed to assure NFIT .ge. 0 in OK cases.
*>> 1987-12-09 DPFIT Lawson Initial code.
* Least squares polynomial fit to discrete data.
* Uses either the Monomial or the Chebyshev basis.
* ------------------------------------------------------------------
* M [integer, in] No. of data points.
*
* (X(I),I=1,M) [float, in] Abcissas of data.
*
* (Y(I),I=1,M) [float, in] Ordinates of data.
*
* (SIG(I),I=1,M) [float, in] Standard deviations of data Y().
* If SIG(1) .lt. 0., the subr will funcrtion as though all
* SIG(I) ate equal to abs(SIG(1)).
* In this latter case SIG() can be dimensioned 1 rather than M.
*
* NMAX [integer, in] Specifies the highest degree polynomial to be
* considered.
*
* SEEKN [logical, in] If .true. this subr will determine the
* optimal degree NFIT in the range [0, NMAX] for fitting the
* given data and compute that fit.
* If .false. this subr will do the fit with NFIT = NMAX
* unless this produces a near-singular problem, in which case
* NFIT will be reduced.
*
* COMTRN [logical, in] If .true. this subr will compute
* transformation parameters, P(1) and P(2), so that the
* transformed abcissa variable ranges from -1.0 to +1.0.
* If .false., will use P(1) and P(2) as given on entry.
*
* CHBBAS [logical, in] If .true. this subr will use the Chebyshev
* basis. If .false., will use the Monomial basis.
*
* (P(J),J=1,NMAX+3) [float, inout] P(1) and P(2) define a
* transformation of the abcissa variable as follows:
*
* S = ( X - P(1) ) / P(2)
*
* (P(I+3),I=0,...,NMAX) are polynomial coefficients computed
* by this subr. P(I+3) is the coeff of S**I, if the Monomial
* bases is used and of the I-th degree Chebyshev polynomial if
* the Chebyshev basis is used.
* If this subr sets NFIT < NMAX then it will also set the
* coefficients P(I+3) for I = NFIT+1, ..., NMAX to zero.
*
* NFIT [integer, out] On a successful return NFIT will be the
* degree of the fitted polynomial. It will be set as
* described above in the specification of SEEKN.
* If input values are inconsistent, NFIT will be set to -1
* to indicate an error condition, and no fit will be done.
*
* SIGFAC [float,out] Factor by which the given SIG() values should
* be multiplied to improve consistency with the fit, i.e.,
* SIGFAC*SIG(i) is the a posteriori estimate of the standard
* deviation of the error in Y(i). In particular,
* if all SIG(i) = 1.0, then SIGFAC is an estimate
* of the standard deviation of the data error.
* Let SUMSQ = the sum from 1 to M of
* ((Residual at ith point)/SIG(i))**2
* Then SIGFAC = sqrt(SUMSQ / max(M-NFIT-1, 1)).
*
* W() [float, work] Working space. Must be dimensioned at least
* (NMAX+3)*(NMAX+3). W() will be used in this subr as a
* 2-dim array: W(NDIM+3,NDIM+3).
* ------------------------------------------------------------------
* C.L.LAWSON, JPL, 1969 DEC 10
* C.L.L., JPL, 1970 JAN 12 Calling sequence changed.
* C.L.L., JPL, 1982 Aug 24:
* To improve portability we are replacing use of
* subrs BHSLR1 and BHSLR2 by SROTG and SROT.
* This uses Givens rotations instead of Householder
* transformations. Accuracy should be the same.
* Execution time will be greater. Storage required
* for the work array W() will be less.
* Name changed from PFIT to LSPOL2.
* C.L.Lawson, JPL, 1984 March 6. Adapted to Fortran 77.
* See type declaration for W(,).
* Counting on the Fortran 77 rule that a DO-loop will
* be skipped if the values of the control parameters
* imply no iterations.
* Name changed from LSPOL2 to [S/D]PFIT,
* 1984 APR 18 Using 'Modified Givens' rather than standard
* Givens to reduce execution time. Requires one more
* column in W().
* 1990-08-10 CLL. In cases of the Y() values lying randomly around
* zero with no polynomial trend, and SEEKN = .true., the preferred
* fit according to our degree determination test may be with no
* coefficients at all.
* In this case the subr formerly set NSOLVE = 0 and NFIT = -1.
* Setting NFIT = -1 is a bad choice on two counts. It indicates an
* error, whereas this is not an error condition. Also it may be
* incompatible with subsequent programs that expect a valid
* polynomial degree for use in polynomial evaluation.
* Changed to set NFIT = max(NSOLVE-1, 0) so we can still set
* NSOLVE = 0 in this case but NFIT will not be set less than 0. */
/* Having NSOLVE = 0 causes all polynomial coefficients to be set to
* zero on return.
* ------------------------------------------------------------------
*--D replaces "?": ?PFIT, ?ERV1, ?ERM1, ?ROTMG, ?ROTM
* Both versions use IERM1, IERV1
* Generic intrinsic functions referenced: SQRT, MAX, MIN, ABS
* ------------------------------------------------------------------
* The dimensions of X(), Y(), and SIG() must be at least M.
* */
/* ------------------------------------------------------------------
*
* N = NMAX = MAX DEGREE TO BE CONSIDERED.
* NP1 = N+1 = NO. OF COEFFS IN A POLY OF DEGREE N
* NP2 = N+2 = COL INDEX FOR y DATA AND ROW INDEX FOR RESIDUAL NORM.
* NP3 = N+3. Col NP3 holds the scale factors for the modified
* Givens method. Row NP3 is used for the entry of additional
* data after the first NP2 data points have been entered.
* Total array space used used in W() is NP3 rows by NP3 cols.
* */
n = nmax;
np1 = n + 1;
np2 = np1 + 1;
np3 = np2 + 1;
if ((n < 0 || m <= 0) || Sig[1] == ZERO)
{
ierm1( "DPFIT", 1, 0, "No fit done. Require NMAX .ge. 0, M > 0, and SIG(1) .ne. 0."
, "NMAX", n, ',' );
ierv1( "M", m, ',' );
derv1( "SIG(1)", Sig[1], '.' );
*nfit = -1;
return;
}
/* D1MACH(4) is the smallest no. that can be added to 1.0
* and will give a no. larger than 1.0 in storage.
* */
teneps = TEN*DBL_EPSILON;
idim = np3;
sigma = fabs( Sig[1] );
/* ------------------------------------------------------------------
* COMPUTE P(1),P(2) IF REQUESTED
*
* CHANGE OF INDEPENDENT VARIABLE IS GIVEN BY S=(X-P(1))/P(2)
* OR X=P(1) + P(2)*S
* */
if (comtrn)
{
xmin = X[1];
xmax = xmin;
for (i = 2; i <= m; i++)
{
xmin = fmin( xmin, X[i] );
xmax = fmax( xmax, X[i] );
}
P[1] = (xmax + xmin)*HALF;
P[2] = (xmax - xmin)*HALF;
if (P[2] == ZERO)
P[2] = ONE;
}
else
{
if (P[2] == ZERO)
{
derm1( "DPFIT", 2, 0, "No fit done. With COMTRN = .FALSE. require P(2) .ne. 0."
, "P(2)", P[2], '.' );
*nfit = -1;
return;
}
}
/* ------------------------------------------------------------------
*
* ACCUMULATION LOOP BEGINS HERE
* IDATA COUNTS THE TOTAL NO. OF DATA POINTS ACCUMULATED.
* LROW is the index of the row of W(,) in which the
* new row of data will be placed.
* */
for (idata = 1; idata <= m; idata++)
{
lrow = min( idata, np3 );
s = (X[idata] - P[1])/P[2];
if (Sig[1] > ZERO)
{
sigma = Sig[idata];
if (sigma <= ZERO)
{
derm1( "DPFIT", 3, 0, "No fit done. With SIG(1) > 0. require all SIG(I) > 0."
, "SIG(1)", Sig[1], ',' );
ierv1( "I", idata, ',' );
derv1( "SIG(I)", sigma, '.' );
*nfit = -1;
return;
}
}
W(0,lrow - 1) = ONE/sigma;
W(np2 - 1,lrow - 1) = Y[idata]/sigma;
W(np3 - 1,lrow - 1) = ONE;
if (n > 0)
{
if (chbbas)
{
/* Chebyshev basis */
W(1,lrow - 1) = s/sigma;
s2 = TWO*s;
for (j = 3; j <= np1; j++)
{
W(j - 1,lrow - 1) = s2*W(j - 2,lrow - 1) - W(j - 3,lrow - 1);
}
}
else
{
/* Monomial basis */
for (j = 2; j <= np1; j++)
{
W(j - 1,lrow - 1) = s*W(j - 2,lrow - 1);
}
}
}
/* Accumulate new data row into triangular array.
* */
for (irow = 1; irow <= (lrow - 1); irow++)
{
drotmg( &W(np3 - 1,irow - 1), &W(np3 - 1,lrow - 1), &W(irow - 1,irow - 1),
W(irow - 1,lrow - 1), param );
if (irow < np2)
drotm( np2 - irow, &W(irow,irow - 1), idim, &W(irow,lrow - 1),
idim, param );
}
}
/* END OF ACCUMULATION LOOP
* ------------------------------------------------------------------
*
* Replace Modified Givens weights by their square roots.
* */
for (i = 1; i <= min( np2, m ); i++)
{
W(np3 - 1,i - 1) = sqrt( W(np3 - 1,i - 1) );
}
/* ------------------------------------------------------------------
*
* Set IRANK = no. of leading diagonal elements of
* the triangular matrix that are not extremely small
* relative to the other elements in the same column.
* */
irank = min( np1, m );
/* Fortran 77 skips this loop if IRANK .le. 1. */
for (j = 2; j <= irank; j++)
{
size = ZERO;
for (ii = 1; ii <= (j - 1); ii++)
{
size = fmax( size, fabs( W(j - 1,ii - 1)*W(np3 - 1,ii - 1) ) );
}
if (fabs( W(j - 1,j - 1)*W(np3 - 1,j - 1) ) < teneps*size)
{
irank = j - 1;
goto L_260;
}
}
L_260:
;
/* ------------------------------------------------------------------
*
* TEMPORARILY COPY RT-SIDE VECTOR INTO P()
* */
for (i = 1; i <= irank; i++)
{
P[i + 2] = W(np2 - 1,i - 1);
}
/* ------------------------------------------------------------------
*
* Now we deal with 3 possible cases.
* We must determine NSOLVE and SIGFAC for each of these cases.
* NSOLVE is the no. of coefficients that will be computed.
* We initially set NSOLVE = IRANK, but NSOLVE may be
* reset to a smaller value in Case 3 below.
* 1. SEEKN = .false. and IRANK = NMAX+1
* This is the simple case.
* 2. SEEKN = .false. and IRANK .lt. NMAX+1
* Requires extra work to compute SIGFAC.
* 3. SEEKN = .true.
* This requires more work to determine NSOLVE and SIGFAC.
* */
nsolve = irank;
if (!seekn && irank == np1)
{
/* Here for Case 1.
* */
temp = m - np1;
if (temp == ZERO)
{
*sigfac = ZERO;
}
else
{
/* Here M .gt. NP1. */
*sigfac = fabs( W(np2 - 1,np2 - 1)*W(np3 - 1,np2 - 1) )/
sqrt( temp );
}
goto L_235;
}
/* ------------------------------------------------------------------
*
* Here for Cases 2 and 3.
*
* Set SIZE = max abs value of elts in col NP2.
* SIZE is used to scale quantities to avoid possible
* trouble due to overflow or underflow.
* */
size = ZERO;
for (i = 1; i <= min( np2, m ); i++)
{
W(np2 - 1,i - 1) = fabs( W(np2 - 1,i - 1)*W(np3 - 1,i - 1) );
size = fmax( size, W(np2 - 1,i - 1) );
}
/* SIZE will be zero if and only if all of
* the given Y() data is zero.
* */
if (size == ZERO)
{
*sigfac = ZERO;
nsolve = 1;
goto L_235;
}
/* ------------------------------------------------------------------
*
* Col NP2 now contains data from which sums of squares of
* residuals can be computed for various possible settings
* of NSOLVE. We will set LIMIT = the largest row index to
* be used in Col NP2 in analyzing residual norms.
* In the usual case of M .ge. NP2, we set LIMIT = NP2.
* Otherwise, when M .le. NP1, we set LIMIT = M+1 and set
* W(LIMIT,NP2) = 0. This reflects the fact that there is the
* possibility of reducing the residual norm to zero by
* exact interpolation when M .le. NP1. The subr will do this
* unless it must set NSOLVE smaller than M due to IRANK being
* less than M.
* Transform col NP2 so the (i+1)-st elt is
* Sum(i) divided by SIZE**2, where Sum(i)
* is the sum of squares of weighted residuals that
* would be obtained if only the first i coefficients
* were computed.
* */
if (m >= np2)
{
limit = np2;
}
else
{
limit = m + 1;
W(np2 - 1,limit - 1) = ZERO;
}
W(np2 - 1,limit - 1) = powi(W(np2 - 1,limit - 1)/size,2);
for (i = limit - 1; i >= 1; i--)
{
W(np2 - 1,i - 1) = powi(W(np2 - 1,i - 1)/size,2) + W(np2 - 1,i);
}
/* ------------------------------------------------------------------
*
* > Do Case 3 if SEEKN is true, and Case 2 if SEEKN is false.
* */
if (seekn)
{
/* > Divide each W(i,NP2) by the no. of degrees of
* > freedom which is M - (i-1).
* > Then set CMIN = smallest of these quotients.
* > Then set NSOLVE.
* */
denom = m;
W(np2 - 1,0) /= denom;
cmin = W(np2 - 1,0);
for (i = 2; i <= (irank + 1); i++)
{
denom = fmax( denom - 1, ONE );
W(np2 - 1,i - 1) /= denom;
cmin = fmin( cmin, W(np2 - 1,i - 1) );
}
temp = FAC*cmin;
for (i = 1; i <= (irank + 1); i++)
{
if (W(np2 - 1,i - 1) <= temp)
{
nsolve = i - 1;
goto L_232;
}
}
L_232:
;
}
else
{
denom = max( m - nsolve, 1 );
W(np2 - 1,nsolve) /= denom;
}
*sigfac = size*sqrt( W(np2 - 1,nsolve) );
/* ------------------------------------------------------------------
*
* Solve for NSOLVE coefficients.
* */
L_235:
;
for (i = nsolve; i >= 1; i--)
{
t = P[i + 2];
/* Fortran 77 will skip this loop when I .EQ. NSOLVE. */
for (j = i + 1; j <= nsolve; j++)
{
t += -W(j - 1,i - 1)*P[j + 2];
}
P[i + 2] = t/W(i - 1,i - 1);
}
/* ------------------------------------------------------------------
*
* Set missing high order coeffs to zero.
*
* Counting on Fortran 77 skipping following loop if
* NSOLVE .GE. NP1.
* */
for (i = nsolve + 1; i <= np1; i++)
{
P[i + 2] = ZERO;
}
*nfit = max( nsolve - 1, 0 );
return;
#undef W
} /* end of function */