/*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 #include "fcrt.h" #include "dpfit.h" #include #include /* 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 */