C from: http://www.netlib.org/gcv/index.html
C "for     B-spline data smoothing using generalized cross-validation
C          and mean squared prediction or explicit user smoothing"
C "by      H.J. Woltring,  University of Nijmegen,
C          Philips Medical Systems, Eindhoven (The Netherlands)"
C
C***********************************************************************
C
C GCVSPL.FOR, 1986-05-12
C
C***********************************************************************
C
C SUBROUTINE GCVSPL (REAL*8)
C
C Purpose:
C *******
C
C       Natural B-spline data smoothing subroutine, using the Generali-
C       zed Cross-Validation and Mean-Squared Prediction Error Criteria
C       of Craven & Wahba (1979). Alternatively, the amount of smoothing
C       can be given explicitly, or it can be based on the effective
C       number of degrees of freedom in the smoothing process as defined
C       by Wahba (1980). The model assumes uncorrelated, additive noise
C       and essentially smooth, underlying functions. The noise may be
C       non-stationary, and the independent co-ordinates may be spaced
C       non-equidistantly. Multiple datasets, with common independent
C       variables and weight factors are accomodated.
C
C
C Calling convention:
C ******************
C
C       CALL GCVSPL ( X, Y, NY, WX, WY, M, N, K, MD, VAL, C, NC, WK, IER )
C
C Meaning of parameters:
C *********************
C
C       X(N)    ( I )   Independent variables: strictly increasing knot
C                       sequence, with X(I-1).lt.X(I), I=2,...,N.
C       Y(NY,K) ( I )   Input data to be smoothed (or interpolated).
C       NY      ( I )   First dimension of array Y(NY,K), with NY.ge.N.
C       WX(N)   ( I )   Weight factor array; WX(I) corresponds with
C                       the relative inverse variance of point Y(I,*).
C                       If no relative weighting information is
C                       available, the WX(I) should be set to ONE.
C                       All WX(I).gt.ZERO, I=1,...,N.
C       WY(K)   ( I )   Weight factor array; WY(J) corresponds with
C                       the relative inverse variance of point Y(*,J).
C                       If no relative weighting information is
C                       available, the WY(J) should be set to ONE.
C                       All WY(J).gt.ZERO, J=1,...,K.
C                       NB: The effective weight for point Y(I,J) is
C                       equal to WX(I)*WY(J).
C       M       ( I )   Half order of the required B-splines (spline
C                       degree 2*M-1), with M.gt.0. The values M =
C                       1,2,3,4 correspond to linear, cubic, quintic,
C                       and heptic splines, respectively.
C       N       ( I )   Number of observations per dataset, with N.ge.2*M.
C       K       ( I )   Number of datasets, with K.ge.1.
C       MD      ( I )   Optimization mode switch:
C                       |MD| = 1: Prior given value for p in VAL
C                                 (VAL.ge.ZERO). This is the fastest
C                                 use of GCVSPL, since no iteration
C                                 is performed in p.
C                       |MD| = 2: Generalized cross validation.
C                       |MD| = 3: True predicted mean-squared error,
C                                 with prior given variance in VAL.
C                       |MD| = 4: Prior given number of degrees of
C                                 freedom in VAL (ZERO.le.VAL.le.N-M).
C                        MD  < 0: It is assumed that the contents of
C                                 X, W, M, N, and WK have not been
C                                 modified since the previous invoca-
C                                 tion of GCVSPL. If MD < -1, WK(4)
C                                 is used as an initial estimate for
C                                 the smoothing parameter p.
C                       Other values for |MD|, and inappropriate values
C                       for VAL will result in an error condition, or
C                       cause a default value for VAL to be selected.
C                       After return from MD.ne.1, the same number of
C                       degrees of freedom can be obtained, for identical
C                       weight factors and knot positions, by selecting
C                       |MD|=1, and by copying the value of p from WK(4)
C                       into VAL. In this way, no iterative optimization
C                       is required when processing other data in Y.
C       VAL     ( I )   Mode value, as described above under MD.
C       C(NC,K) ( O )   Spline coefficients, to be used in conjunction
C                       with function SPLDER. NB: the dimensions of C
C                       in GCVSPL and in SPLDER are different
C                       only a single column of C(N,K) is needed, and the
C                       proper column C(1,J), with J=1...K should be used
C                       when calling SPLDER.
C       NC       ( I )  First dimension of array C(NC,K), NC.ge.N.
C       WK(IWK) (I/W/O) Work vector, with length IWK.ge.6*(N*M+1)+N.
C                       On normal exit, the first 6 values of WK are
C                       assigned as follows:
C
C                       WK(1) = Generalized Cross Validation value
C                       WK(2) = Mean Squared Residual.
C                       WK(3) = Estimate of the number of degrees of
C                               freedom of the residual sum of squares
C                               per dataset, with 0.lt.WK(3).lt.N-M.
C                       WK(4) = Smoothing parameter p, multiplicative
C                               with the splines' derivative constraint.
C                       WK(5) = Estimate of the true mean squared error
C                               (different formula for |MD| = 3).
C                       WK(6) = Gauss-Markov error variance.
C
C                       If WK(4) -->  0 , WK(3) -->  0 , and an inter-
C                       polating spline is fitted to the data (p --> 0).
C                       A very small value > 0 is used for p, in order
C                       to avoid division by zero in the GCV function.
C
C                       If WK(4) --> inf, WK(3) --> N-M, and a least-
C                       squares polynomial of order M (degree M-1) is
C                       fitted to the data (p --> inf). For numerical
C                       reasons, a very high value is used for p.
C
C                       Upon return, the contents of WK can be used for
C                       covariance propagation in terms of the matrices
C                       B and WE: see the source listings. The variance
C                       estimate for dataset J follows as WK(6)/WY(J).
C
C       IER     ( O )   Error parameter:
C
C                       IER = 0:        Normal exit
C                       IER = 1:        M.le.0 .or. N.lt.2*M
C                       IER = 2:        Knot sequence is not strictly
C                                       increasing, or some weight
C                                       factor is not positive.
C                       IER = 3:        Wrong mode  parameter or value.
C
C Remarks:
C *******
C
C       (1) GCVSPL calculates a natural spline of order 2*M (degree
C       2*M-1) which smoothes or interpolates a given set of data
C       points, using statistical considerations to determine the
C       amount of smoothing required (Craven & Wahba, 1979). If the
C       error variance is a priori known, it should be supplied to
C       the routine in VAL, for |MD|=3. The degree of smoothing is
C       then determined to minimize an unbiased estimate of the true
C       mean squared error. On the other hand, if the error variance
C       is not known, one may select |MD|=2. The routine then deter-
C       mines the degree of smoothing to minimize the generalized
C       cross validation function. This is asymptotically the same
C       as minimizing the true predicted mean squared error (Craven &
C       Wahba, 1979). If the estimates from |MD|=2 or 3 do not appear
C       suitable to the user (as apparent from the smoothness of the
C       M-th derivative or from the effective number of degrees of
C       freedom returned in WK(3) ), the user may select an other
C       value for the noise variance if |MD|=3, or a reasonably large
C       number of degrees of freedom if |MD|=4. If |MD|=1, the proce-
C       dure is non-iterative, and returns a spline for the given
C       value of the smoothing parameter p as entered in VAL.
C
C       (2) The number of arithmetic operations and the amount of
C       storage required are both proportional to N, so very large
C       datasets may be accomodated. The data points do not have
C       to be equidistant in the independant variable X or uniformly
C       weighted in the dependant variable Y. However, the data
C       points in X must be strictly increasing. Multiple dataset
C       processing (K.gt.1) is numerically more efficient dan
C       separate processing of the individual datasets (K.eq.1).
C
C       (3) If |MD|=3 (a priori known noise variance), any value of
C       N.ge.2*M is acceptable. However, it is advisable for N-2*M
C       be rather large (at least 20) if |MD|=2 (GCV).
C
C       (4) For |MD| > 1, GCVSPL tries to iteratively minimize the
C       selected criterion function. This minimum is unique for |MD|
C       = 4, but not necessarily for |MD| = 2 or 3. Consequently, 
C       local optima rather that the global optimum might be found,
C       and some actual findings suggest that local optima might
C       yield more meaningful results than the global optimum if N
C       is small. Therefore, the user has some control over the
C       search procedure. If MD > 1, the iterative search starts
C       from a value which yields a number of degrees of freedom
C       which is approximately equal to N/2, until the first (local)
C       minimum is found via a golden section search procedure
C       (Utreras, 1980). If MD < -1, the value for p contained in
C       WK(4) is used instead. Thus, if MD = 2 or 3 yield too noisy
C       an estimate, the user might try |MD| = 1 or 4, for suitably
C       selected values for p or for the number of degrees of
C       freedom, and then run GCVSPL with MD = -2 or -3. The con-
C       tents of N, M, K, X, WX, WY, and WK are assumed unchanged
C       if MD < 0.
C
C       (5) GCVSPL calculates the spline coefficient array C(N,K);
C       this array can be used to calculate the spline function
C       value and any of its derivatives up to the degree 2*M-1
C       at any argument T within the knot range, using subrou-
C       tines SPLDER and SEARCH, and the knot array X(N). Since
C       the splines are constrained at their Mth derivative, only
C       the lower spline derivatives will tend to be reliable
C       estimates of the underlying, true signal derivatives.
C
C       (6) GCVSPL combines elements of subroutine CRVO5 by Utre-
C       ras (1980), subroutine SMOOTH by Lyche et al. (1983), and
C       subroutine CUBGCV by Hutchinson (1985). The trace of the
C       influence matrix is assessed in a similar way as described
C       by Hutchinson & de Hoog (1985). The major difference is
C       that the present approach utilizes non-symmetrical B-spline
C       design matrices as described by Lyche et al. (1983); there-
C       fore, the original algorithm by Erisman & Tinney (1975) has
C       been used, rather than the symmetrical version adopted by
C       Hutchinson & de Hoog.
C
C References:
C **********
C
C       P. Craven & G. Wahba (1979), Smoothing noisy data with
C       spline functions. Numerische Mathematik 31, 377-403.
C
C       A.M. Erisman & W.F. Tinney (1975), On computing certain
C       elements of the inverse of a sparse matrix. Communications
C       of the ACM 18(3), 177-179.
C
C       M.F. Hutchinson & F.R. de Hoog (1985), Smoothing noisy data
C       with spline functions. Numerische Mathematik 47(1), 99-106.
C
C       M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO Division of
C       Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601,
C       Australia.
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori (1983), Fortran
C       subroutines for computing smoothing and interpolating natural
C       splines. Advances in Engineering Software 5(1), 2-5.
C
C       F. Utreras (1980), Un paquete de programas para ajustar curvas
C       mediante funciones spline. Informe Tecnico MA-80-B-209, Depar-
C       tamento de Matematicas, Faculdad de Ciencias Fisicas y Matema-
C       ticas, Universidad de Chile, Santiago.
C
C       Wahba, G. (1980). Numerical and statistical methods for mildly,
C       moderately and severely ill-posed problems with noisy data.
C       Technical report nr. 595 (February 1980). Department of Statis-
C       tics, University of Madison (WI), U.S.A.
C
C Subprograms required:
C ********************
C
C       BASIS, PREP, SPLC, BANDET, BANSOL, TRINV
C
C***********************************************************************
C
      SUBROUTINE GCVSPL ( X, Y, NY, WX, WY, M, N, K, MD, VAL, C, NC,
     1                   WK, IER )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( RATIO=2D0, TAU=1.618033983D0, IBWE=7,
     1           ZERO=0D0, HALF=5D-1 , ONE=1D0, TOL=1D-6,
     2           EPS=1D-15, EPSINV=ONE/EPS )
      DIMENSION X(N), Y(NY,K), WX(N), WY(K), C(NC,K), WK(N+6*(N*M+1))
      SAVE M2, NM1, EL
      DATA M2, NM1, EL / 2*0, 0D0 /
C
C***  Parameter check and work array initialization
C
      IER = 0
C***  Check on mode parameter
      IF ((IABS(MD).GT.4) .OR.(  MD.EQ. 0  ) .OR.
     1  ((IABS(MD).EQ.1).AND.( VAL.LT.ZERO)).OR.
     2  ((IABS(MD).EQ.3).AND.( VAL.LT.ZERO)).OR.
     3  ((IABS(MD).EQ.4).AND.((VAL.LT.ZERO) .OR.(VAL.GT.N-M)))) THEN
         IER = 3      
         RETURN
      ENDIF
C***  Check on M and N
      IF (MD.GT.0) THEN
         M2  = 2 * M
         NM1 = N - 1
      ELSE
         IF ((M2.NE.2*M).OR.(NM1.NE.N-1)) THEN
            IER = 3      
            RETURN
         ENDIF
      ENDIF
      IF ((M.LE.0).OR.(N.LT.M2)) THEN
         IER = 1      
         RETURN
      ENDIF
C***  Check on knot sequence and weights
      IF (WX(1).LE.ZERO) IER = 2
      DO 10 I=2,N
         IF ((WX(I).LE.ZERO).OR.(X(I-1).GE.X(I))) IER = 2
         IF (IER.NE.0) RETURN
   10 CONTINUE
      DO 15 J=1,K
         IF (WY(J).LE.ZERO) IER = 2
         IF (IER.NE.0) RETURN
   15 CONTINUE
C
C***  Work array parameters (address information for covariance 
C***  propagation by means of the matrices STAT, B, and WE). NB:
C***  BWE cannot be used since it is modified by function TRINV.
C
      NM2P1 = N*(M2+1)
      NM2M1 = N*(M2-1)
C     ISTAT = 1            
C     IBWE  = ISTAT + 6      
      IB    = IBWE  + NM2P1      
      IWE   = IB    + NM2M1      
C     IWK   = IWE   + NM2P1      
C
C***  Compute the design matrices B and WE, the ratio
C***  of their L1-norms, and check for iterative mode.
C
      IF (MD.GT.0) THEN
         CALL BASIS ( M, N, X, WK(IB), R1, WK(IBWE) )
         CALL PREP  ( M, N, X, WX, WK(IWE), EL )
         EL = EL / R1      
      ENDIF
      IF (IABS(MD).NE.1) GO TO 20
C***     Prior given value for p
         R1 = VAL
         GO TO 100
C
C***  Iterate to minimize the GCV function (|MD|=2),
C***  the MSE function (|MD|=3), or to obtain the prior
C***  given number of degrees of freedom (|MD|=4).
C
   20 IF (MD.LT.-1) THEN
         R1 = WK(4)      
      ELSE
         R1 = ONE / EL      
      ENDIF      
      R2 = R1 * RATIO
      GF2 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R2,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
   40 GF1 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R1,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
      IF (GF1.GT.GF2) GO TO 50
         IF (WK(4).LE.ZERO) GO TO 100            
         R2  = R1
         GF2 = GF1
         R1  = R1 / RATIO
         GO TO 40
   50 R3 = R2 * RATIO
   60 GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
      IF (GF3.GT.GF2) GO TO 70
         IF (WK(4).GE.EPSINV) GO TO 100      
         R2  = R3      
         GF2 = GF3
         R3  = R3 * RATIO
         GO TO 60
   70 R2  = R3
      GF2 = GF3
      ALPHA = (R2-R1) / TAU
      R4 = R1 + ALPHA
      R3 = R2 - ALPHA
      GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
      GF4 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R4,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
   80 IF (GF3.LE.GF4) THEN
         R2  = R4
         GF2 = GF4
         ERR = (R2-R1) / (R1+R2)
         IF ((ERR*ERR+ONE.EQ.ONE).OR.(ERR.LE.TOL)) GO TO 90
         R4  = R3
         GF4 = GF3
         ALPHA = ALPHA / TAU
         R3  = R2 - ALPHA
         GF3 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R3,EPS,C,NC,
     1             WK,WK(IB),WK(IWE),EL,WK(IBWE))
      ELSE
         R1  = R3
         GF1 = GF3
         ERR = (R2-R1) / (R1+R2)
         IF ((ERR*ERR+ONE.EQ.ONE).OR.(ERR.LE.TOL)) GO TO 90
         R3  = R4
         GF3 = GF4
         ALPHA = ALPHA / TAU
         R4 = R1 + ALPHA
         GF4 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R4,EPS,C,NC,
     1             WK,WK(IB),WK(IWE),EL,WK(IBWE))
      ENDIF
      GO TO 80
   90 R1 = HALF * (R1+R2)
C
C***  Calculate final spline coefficients
C
  100 GF1 = SPLC(M,N,K,Y,NY,WX,WY,MD,VAL,R1,EPS,C,NC,
     1          WK,WK(IB),WK(IWE),EL,WK(IBWE))
C
C***  Ready
C
      RETURN
      END
C BASIS.FOR, 1985-06-03
C
C***********************************************************************
C
C SUBROUTINE BASIS (REAL*8)
C
C Purpose:
C *******
C
C       Subroutine to assess a B-spline tableau, stored in vectorized
C       form.
C
C Calling convention:
C ******************
C
C       CALL BASIS ( M, N, X, B, BL, Q )
C
C Meaning of parameters:
C *********************
C
C       M               ( I )   Half order of the spline (degree 2*M-1),
C                               M > 0.
C       N               ( I )   Number of knots, N >= 2*M.
C       X(N)            ( I )   Knot sequence, X(I-1) < X(I), I=2,N.
C       B(1-M:M-1,N)    ( O )   Output tableau. Element B(J,I) of array
C                               B corresponds with element b(i,i+j) of
C                               the tableau matrix B.
C       BL              ( O )   L1-norm of B.
C       Q(1-M:M)        ( W )   Internal work array.
C
C Remark:
C ******
C
C       This subroutine is an adaptation of subroutine BASIS from the
C       paper by Lyche et al. (1983). No checking is performed on the
C       validity of M and N. If the knot sequence is not strictly in-
C       creasing, division by zero may occur.
C
C Reference:
C *********
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      SUBROUTINE BASIS ( M, N, X, B, BL, Q )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( ZERO=0D0, ONE=1D0 )
      DIMENSION X(N), B(1-M:M-1,N), Q(1-M:M)
C
      IF (M.EQ.1) THEN
C***         Linear spline
         DO 3 I=1,N
            B(0,I) = ONE
    3    CONTINUE
         BL = ONE
         RETURN
      ENDIF
C
C***  General splines
C
      MM1 = M - 1
      MP1 = M + 1
      M2  = 2 * M
      DO 15 L=1,N
C***     1st row
         DO 5 J=-MM1,M
            Q(J) = ZERO
    5    CONTINUE
         Q(MM1) = ONE
         IF ((L.NE.1).AND.(L.NE.N))
     1      Q(MM1) = ONE / ( X(L+1) - X(L-1) )
C***     Successive rows
         ARG = X(L)
         DO 13 I=3,M2
            IR = MP1 - I
            V  = Q(IR)
            IF (L.LT.I) THEN
C***               Left-hand B-splines
               DO 6 J=L+1,I
                  U     = V
                  V     = Q(IR+1)
                  Q(IR) = U + (X(J)-ARG)*V
                  IR    = IR + 1
    6          CONTINUE
            ENDIF
            J1 = MAX0(L-I+1,1)
            J2 = MIN0(L-1,N-I)
            IF (J1.LE.J2) THEN
C***               Ordinary B-splines
               IF (I.LT.M2) THEN
                  DO 8 J=J1,J2
                     Y     = X(I+J)
                     U     = V
                     V     = Q(IR+1)
                     Q(IR) = U + (V-U)*(Y-ARG)/(Y-X(J))
                     IR = IR + 1
    8             CONTINUE
               ELSE
                  DO 10 J=J1,J2
                     U     = V
                     V     = Q(IR+1)
                     Q(IR) = (ARG-X(J))*U + (X(I+J)-ARG)*V
                     IR    = IR + 1
   10             CONTINUE
               ENDIF
            ENDIF
            NMIP1 = N - I + 1
            IF (NMIP1.LT.L) THEN
C***           Right-hand B-splines
               DO 12 J=NMIP1,L-1
                  U     = V
                  V     = Q(IR+1)
                  Q(IR) = (ARG-X(J))*U + V
                  IR    = IR + 1
   12          CONTINUE
            ENDIF
   13    CONTINUE
         DO 14 J=-MM1,MM1
            B(J,L) = Q(J)
   14    CONTINUE
   15 CONTINUE
C
C***  Zero unused parts of B
C
      DO 17 I=1,MM1
         DO 16 K=I,MM1
            B(-K,    I) = ZERO
            B( K,N+1-I) = ZERO
   16    CONTINUE
   17 CONTINUE
C
C***  Assess L1-norm of B
C
      BL = 0D0
      DO 19 I=1,N
         DO 18 K=-MM1,MM1
            BL = BL + ABS(B(K,I))
   18    CONTINUE
   19 CONTINUE
      BL = BL / N
C
C***  Ready
C
      RETURN
      END
C PREP.FOR, 1985-07-04
C
C***********************************************************************
C
C SUBROUTINE PREP (REAL*8)
C
C Purpose:
C *******
C
C       To compute the matrix WE of weighted divided difference coeffi-
C       cients needed to set up a linear system of equations for sol-
C       ving B-spline smoothing problems, and its L1-norm EL. The matrix
C       WE is stored in vectorized form.
C
C Calling convention:
C ******************
C
C       CALL PREP ( M, N, X, W, WE, EL )
C
C Meaning of parameters:
C *********************
C
C       M               ( I )   Half order of the B-spline (degree
C                               2*M-1), with M > 0.
C       N               ( I )   Number of knots, with N >= 2*M.
C       X(N)            ( I )   Strictly increasing knot array, with
C                               X(I-1) < X(I), I=2,N.
C       W(N)            ( I )   Weight matrix (diagonal), with
C                               W(I).gt.0.0, I=1,N.
C       WE(-M:M,N)      ( O )   Array containing the weighted divided
C                               difference terms in vectorized format.
C                               Element WE(J,I) of array E corresponds
C                               with element e(i,i+j) of the matrix
C                               W**-1 * E.
C       EL              ( O )   L1-norm of WE.
C
C Remark:
C ******
C
C       This subroutine is an adaptation of subroutine PREP from the paper
C       by Lyche et al. (1983). No checking is performed on the validity
C       of M and N. Division by zero may occur if the knot sequence is
C       not strictly increasing.
C
C Reference:
C *********
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      SUBROUTINE PREP ( M, N, X, W, WE, EL )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( ZERO=0D0, ONE=1D0 )
      DIMENSION X(N), W(N), WE((2*M+1)*N)      
C
C***  Calculate the factor F1
C
      M2   = 2 * M
      MP1  = M + 1
      M2M1 = M2 - 1
      M2P1 = M2 + 1
      NM   = N - M
      F1   = -ONE
      IF (M.NE.1) THEN
         DO 5 I=2,M
            F1 = -F1 * I
    5    CONTINUE
         DO 6 I=MP1,M2M1
            F1 = F1 * I
    6    CONTINUE
      END IF
C
C***  Columnwise evaluation of the unweighted design matrix E
C
      I1 = 1
      I2 = M
      JM = MP1
      DO 17 J=1,N
         INC = M2P1
         IF (J.GT.NM) THEN
            F1 = -F1
            F  =  F1
         ELSE
            IF (J.LT.MP1) THEN
                INC = 1
                F   = F1
            ELSE
                F   = F1 * (X(J+M)-X(J-M))
            END IF
         END IF
         IF ( J.GT.MP1) I1 = I1 + 1
         IF (I2.LT.  N) I2 = I2 + 1
         JJ = JM
C***     Loop for divided difference coefficients
         FF = F
         Y = X(I1)
         I1P1 = I1 + 1
         DO 11 I=I1P1,I2
            FF = FF / (Y-X(I))
   11    CONTINUE
         WE(JJ) = FF
         JJ = JJ + M2
         I2M1 = I2 - 1
         IF (I1P1.LE.I2M1) THEN
            DO 14 L=I1P1,I2M1
               FF = F
               Y  = X(L)
               DO 12 I=I1,L-1
                  FF = FF / (Y-X(I))
   12          CONTINUE
               DO 13 I=L+1,I2
                  FF = FF / (Y-X(I))
   13          CONTINUE
               WE(JJ) = FF
               JJ = JJ + M2
   14       CONTINUE
         END IF
         FF = F
         Y = X(I2)
         DO 16 I=I1,I2M1
            FF = FF / (Y-X(I))
   16    CONTINUE
         WE(JJ) = FF
         JJ = JJ + M2
         JM = JM + INC
   17 CONTINUE
C
C***  Zero the upper left and lower right corners of E
C
      KL = 1
      N2M = M2P1*N + 1
      DO 19 I=1,M
         KU = KL + M - I
         DO 18 K=KL,KU
            WE(    K) = ZERO
            WE(N2M-K) = ZERO
   18    CONTINUE
         KL = KL + M2P1
   19 CONTINUE
C
C***  Weighted matrix WE = W**-1 * E and its L1-norm
C
   20 JJ = 0
      EL = 0D0
      DO 22 I=1,N
         WI = W(I)
         DO 21 J=1,M2P1
            JJ     = JJ + 1
            WE(JJ) = WE(JJ) / WI
            EL     = EL + ABS(WE(JJ))
   21    CONTINUE
   22 CONTINUE
      EL = EL / N
C
C***  Ready
C
      RETURN
      END
C SPLC.FOR, 1985-12-12
C
C Author: H.J. Woltring
C
C Organizations: University of Nijmegen, and
C                Philips Medical Systems, Eindhoven
C                (The Netherlands)
C
C***********************************************************************
C
C FUNCTION SPLC (REAL*8)
C
C Purpose:
C *******
C
C       To assess the coefficients of a B-spline and various statistical
C       parameters, for a given value of the regularization parameter p.
C
C Calling convention:
C ******************
C
C       FV = SPLC ( M, N, K, Y, NY, WX, WY, MODE, VAL, P, EPS, C, NC,
C       1           STAT, B, WE, EL, BWE)
C
C Meaning of parameters:
C *********************
C
C       SPLC            ( O )   GCV function value if |MODE|.eq.2,
C                               MSE value if |MODE|.eq.3, and absolute
C                               difference with the prior given number of
C                               degrees of freedom if |MODE|.eq.4.
C       M               ( I )   Half order of the B-spline (degree 2*M-1),
C                               with M > 0.
C       N               ( I )   Number of observations, with N >= 2*M.
C       K               ( I )   Number of datasets, with K >= 1.
C       Y(NY,K)         ( I )   Observed measurements.
C       NY              ( I )   First dimension of Y(NY,K), with NY.ge.N.
C       WX(N)           ( I )   Weight factors, corresponding to the
C                               relative inverse variance of each measure-
C                               ment, with WX(I) > 0.0.
C       WY(K)           ( I )   Weight factors, corresponding to the
C                               relative inverse variance of each dataset,
C                               with WY(J) > 0.0.
C       MODE            ( I )   Mode switch, as described in GCVSPL.
C       VAL             ( I )   Prior variance if |MODE|.eq.3, and
C                               prior number of degrees of freedom if
C                               |MODE|.eq.4. For other values of MODE,
C                               VAL is not used.
C       P               ( I )   Smoothing parameter, with P >= 0.0. If
C                               P.eq.0.0, an interpolating spline is
C                               calculated.
C       EPS             ( I )   Relative rounding tolerance*10.0. EPS is
C                               the smallest positive number such that
C                               EPS/10.0 + 1.0 .ne. 1.0.
C       C(NC,K)         ( O )   Calculated spline coefficient arrays. NB:
C                               the dimensions of in GCVSPL and in SPLDER
C                               are different
C                               column of C(N,K) is needed, and the proper
C                               column C(1,J), with J=1...K, should be used
C                               when calling SPLDER.
C       NC              ( I )   First dimension of C(NC,K), with NC.ge.N.
C       STAT(6)         ( O )   Statistics array. See the description in
C                               subroutine GCVSPL.
C       B (1-M:M-1,N)   ( I )   B-spline tableau as evaluated by subroutine
C                               BASIS.
C       WE( -M:M  ,N)   ( I )   Weighted B-spline tableau (W**-1 * E) as
C                               evaluated by subroutine PREP.
C       EL              ( I )   L1-norm of the matrix WE as evaluated by
C                               subroutine PREP.
C       BWE(-M:M,N)     ( O )   Central 2*M+1 bands of the inverted
C                               matrix ( B  +  p * W**-1 * E )**-1
C
C Remarks:
C *******
C
C       This subroutine combines elements of subroutine SPLC0 from the
C       paper by Lyche et al. (1983), and of subroutine SPFIT1 by
C       Hutchinson (1985).
C
C References:
C **********
C
C       M.F. Hutchinson (1985), Subroutine CUBGCV. CSIRO division of
C       Mathematics and Statistics, P.O. Box 1965, Canberra, ACT 2601,
C       Australia.
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      FUNCTION SPLC( M, N, K, Y, NY, WX, WY, MODE, VAL, P, EPS, C, NC,
     1              STAT, B, WE, EL, BWE)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( ZERO=0D0, ONE=1D0, TWO=2D0 )
      DIMENSION Y(NY,K), WX(N), WY(K), C(NC,K), STAT(6),
     1         B(1-M:M-1,N), WE(-M:M,N), BWE(-M:M,N)
C
C***  Check on p-value
C
      DP = P
      STAT(4) = P
      PEL = P * EL
C***  Pseudo-interpolation if p is too small
      IF (PEL.LT.EPS) THEN
         DP = EPS / EL
         STAT(4) = ZERO
      ENDIF
C***  Pseudo least-squares polynomial if p is too large
      IF (PEL*EPS.GT.ONE) THEN
         DP = ONE / (EL*EPS)
         STAT(4) = DP
      ENDIF
C
C***  Calculate  BWE  =  B  +  p * W**-1 * E
C
      DO 40 I=1,N
         KM = -MIN0(M,I-1)
         KP =  MIN0(M,N-I)
         DO 30 L=KM,KP
            IF (IABS(L).EQ.M) THEN
               BWE(L,I) =          DP * WE(L,I)
            ELSE
               BWE(L,I) = B(L,I) + DP * WE(L,I)
            ENDIF
   30    CONTINUE
   40 CONTINUE
C
C***  Solve BWE * C = Y, and assess TRACE [ B * BWE**-1 ]
C
      CALL BANDET ( BWE, M, N )
      CALL BANSOL ( BWE, Y, NY, C, NC, M, N, K )
      STAT(3) = TRINV ( WE, BWE, M, N ) * DP      
      TRN = STAT(3) / N
C
C***  Compute mean-squared weighted residual
C
      ESN = ZERO
      DO 70 J=1,K
         DO 60 I=1,N
            DT = -Y(I,J)
            KM = -MIN0(M-1,I-1)
            KP =  MIN0(M-1,N-I)
            DO 50 L=KM,KP
               DT = DT + B(L,I)*C(I+L,J)
   50       CONTINUE
            ESN = ESN + DT*DT*WX(I)*WY(J)
   60    CONTINUE
   70 CONTINUE
      ESN = ESN / (N*K)
C
C***  Calculate statistics and function value
C
      STAT(6) = ESN / TRN             
      STAT(1) = STAT(6) / TRN         
      STAT(2) = ESN                   
C     STAT(3) = trace [p*B * BWE**-1] 
C     STAT(4) = P                     
      IF (IABS(MODE).NE.3) THEN
C***     Unknown variance: GCV
         STAT(5) = STAT(6) - ESN
         IF (IABS(MODE).EQ.1) SPLC = ZERO
         IF (IABS(MODE).EQ.2) SPLC = STAT(1)
         IF (IABS(MODE).EQ.4) SPLC = DABS( STAT(3) - VAL )
      ELSE
C***     Known variance: estimated mean squared error
         STAT(5) = ESN - VAL*(TWO*TRN - ONE)
         SPLC = STAT(5)
      ENDIF
C
      RETURN
      END
C BANDET.FOR, 1985-06-03
C
C***********************************************************************
C
C SUBROUTINE BANDET (REAL*8)
C
C Purpose:
C *******
C
C       This subroutine computes the LU decomposition of an N*N matrix
C       E. It is assumed that E has M bands above and M bands below the
C       diagonal. The decomposition is returned in E. It is assumed that
C       E can be decomposed without pivoting. The matrix E is stored in
C       vectorized form in the array E(-M:M,N), where element E(J,I) of
C       the array E corresponds with element e(i,i+j) of the matrix E.
C
C Calling convention:
C ******************
C
C       CALL BANDET ( E, M, N )
C
C Meaning of parameters:
C *********************
C
C       E(-M:M,N)       (I/O)   Matrix to be decomposed.
C       M, N            ( I )   Matrix dimensioning parameters,
C                               M >= 0, N >= 2*M.
C
C Remark:
C ******
C
C       No checking on the validity of the input data is performed.
C       If (M.le.0), no action is taken.
C
C***********************************************************************
C
      SUBROUTINE BANDET ( E, M, N )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION E(-M:M,N)
C
      IF (M.LE.0) RETURN
      DO 40 I=1,N
         DI = E(0,I)
         MI = MIN0(M,I-1)
         IF (MI.GE.1) THEN
            DO 10 K=1,MI
               DI = DI - E(-K,I)*E(K,I-K)
   10       CONTINUE
            E(0,I) = DI
         ENDIF
         LM = MIN0(M,N-I)
         IF (LM.GE.1) THEN
            DO 30 L=1,LM
               DL = E(-L,I+L)
               KM = MIN0(M-L,I-1)
               IF (KM.GE.1) THEN
                  DU = E(L,I)
                  DO 20 K=1,KM
                     DU = DU - E(  -K,  I)*E(L+K,I-K)
                     DL = DL - E(-L-K,L+I)*E(  K,I-K)
   20             CONTINUE
                  E(L,I) = DU
               ENDIF
               E(-L,I+L) = DL / DI
   30       CONTINUE
         ENDIF
   40 CONTINUE
C
C***  Ready
C
      RETURN
      END
C BANSOL.FOR, 1985-12-12
C
C***********************************************************************
C
C SUBROUTINE BANSOL (REAL*8)
C
C Purpose:
C *******
C
C       This subroutine solves systems of linear equations given an LU
C       decomposition of the design matrix. Such a decomposition is pro-
C       vided by subroutine BANDET, in vectorized form. It is assumed
C       that the design matrix is not singular. 
C
C Calling convention:
C ******************
C
C       CALL BANSOL ( E, Y, NY, C, NC, M, N, K )
C
C Meaning of parameters:
C *********************
C
C       E(-M:M,N)       ( I )   Input design matrix, in LU-decomposed,
C                               vectorized form. Element E(J,I) of the
C                               array E corresponds with element
C                               e(i,i+j) of the N*N design matrix E.
C       Y(NY,K)         ( I )   Right hand side vectors.
C       C(NC,K)         ( O )   Solution vectors.
C       NY, NC, M, N, K ( I )   Dimensioning parameters, with M >= 0,
C                               N > 2*M, and K >= 1.
C
C Remark:
C ******
C
C       This subroutine is an adaptation of subroutine BANSOL from the
C       paper by Lyche et al. (1983). No checking is performed on the
C       validity of the input parameters and data. Division by zero may
C       occur if the system is singular.
C
C Reference:
C *********
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      SUBROUTINE BANSOL ( E, Y, NY, C, NC, M, N, K )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION E(-M:M,N), Y(NY,K), C(NC,K)
C
C***  Check on special cases: M=0, M=1, M>1
C
      NM1 = N - 1
      IF (M-1) 10,40,80
C
C***  M = 0: Diagonal system
C
   10 DO 30 I=1,N
         DO 20 J=1,K
            C(I,J) = Y(I,J) / E(0,I)
   20    CONTINUE
   30 CONTINUE
      RETURN
C
C***  M = 1: Tridiagonal system
C
   40 DO 70 J=1,K
         C(1,J) = Y(1,J)
         DO 50 I=2,N            
            C(I,J) =  Y(I,J) - E(-1,I)*C(I-1,J)
   50      CONTINUE
         C(N,J) = C(N,J) / E(0,N)
         DO 60 I=NM1,1,-1      
            C(I,J) = (C(I,J) - E( 1,I)*C(I+1,J)) / E(0,I)
   60    CONTINUE
   70 CONTINUE
      RETURN
C
C***  M > 1: General system
C
   80 DO 130 J=1,K
         C(1,J) = Y(1,J)
         DO 100 I=2,N            
            MI = MIN0(M,I-1)
            D  = Y(I,J)
            DO 90 L=1,MI
               D = D - E(-L,I)*C(I-L,J)
   90       CONTINUE
            C(I,J) = D
  100    CONTINUE
         C(N,J) = C(N,J) / E(0,N)
         DO 120 I=NM1,1,-1      
            MI = MIN0(M,N-I)
            D  = C(I,J)
            DO 110 L=1,MI
               D = D - E( L,I)*C(I+L,J)
  110       CONTINUE
            C(I,J) = D / E(0,I)
  120    CONTINUE
  130 CONTINUE
      RETURN
C
      END
C TRINV.FOR, 1985-06-03
C
C***********************************************************************
C
C FUNCTION TRINV (REAL*8)
C
C Purpose:
C *******
C
C       To calculate TRACE [ B * E**-1 ], where B and E are N * N
C       matrices with bandwidth 2*M+1, and where E is a regular matrix
C       in LU-decomposed form. B and E are stored in vectorized form,
C       compatible with subroutines BANDET and BANSOL.
C
C Calling convention:
C ******************
C
C       TRACE = TRINV ( B, E, M, N )
C
C Meaning of parameters:
C *********************
C
C       B(-M:M,N)       ( I ) Input array for matrix B. Element B(J,I)
C                             corresponds with element b(i,i+j) of the
C                             matrix B.
C       E(-M:M,N)       (I/O) Input array for matrix E. Element E(J,I)
C                             corresponds with element e(i,i+j) of the
C                             matrix E. This matrix is stored in LU-
C                             decomposed form, with L unit lower tri-
C                             angular, and U upper triangular. The unit
C                             diagonal of L is not stored. Upon return,
C                             the array E holds the central 2*M+1 bands
C                             of the inverse E**-1, in similar ordering.
C       M, N            ( I ) Array and matrix dimensioning parameters
C                             (M.gt.0, N.ge.2*M+1).
C       TRINV           ( O ) Output function value TRACE [ B * E**-1 ]
C
C Reference:
C *********
C
C       A.M. Erisman & W.F. Tinney, On computing certain elements of the
C       inverse of a sparse matrix. Communications of the ACM 18(1975),
C       nr. 3, pp. 177-179.
C
C***********************************************************************
C
      REAL*8 FUNCTION TRINV ( B, E, M, N )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( ZERO=0D0, ONE=1D0 )
      DIMENSION B(-M:M,N), E(-M:M,N)
C
C***  Assess central 2*M+1 bands of E**-1 and store in array E
C
      E(0,N) = ONE / E(0,N)      
      DO 40 I=N-1,1,-1
         MI = MIN0(M,N-I)
         DD  = ONE / E(0,I)      
C***     Save Ith column of L and Ith row of U, and normalize U row
         DO 10 K=1,MI
            E( K,N) = E( K,  I) * DD      
            E(-K,1) = E(-K,K+I)      
   10    CONTINUE
         DD = DD + DD
C***     Invert around Ith pivot
         DO 30 J=MI,1,-1
            DU = ZERO
            DL = ZERO
            DO 20 K=1,MI
               DU = DU - E( K,N)*E(J-K,I+K)
               DL = DL - E(-K,1)*E(K-J,I+J)
   20       CONTINUE
            E( J,  I) = DU
            E(-J,J+I) = DL
            DD = DD - (E(J,N)*DL + E(-J,1)*DU)
   30    CONTINUE
         E(0,I) = 5D-1 * DD
   40 CONTINUE
C
C***  Assess TRACE [ B * E**-1 ] and clear working storage
C
      DD = ZERO
      DO 60 I=1,N
         MN = -MIN0(M,I-1)
         MP =  MIN0(M,N-I)
         DO 50 K=MN,MP
            DD = DD + B(K,I)*E(-K,K+I)
   50    CONTINUE
   60 CONTINUE
      TRINV = DD
      DO 70 K=1,M
         E( K,N) = ZERO
         E(-K,1) = ZERO
   70 CONTINUE
C
C***  Ready
C
      RETURN
      END
C SPLDER.FOR, 1985-06-11
C
C***********************************************************************
C
C FUNCTION SPLDER (REAL*8)
C
C Purpose:
C *******
C
C       To produce the value of the function (IDER.eq.0) or of the
C       IDERth derivative (IDER.gt.0) of a 2M-th order B-spline at
C       the point T. The spline is described in terms of the half
C       order M, the knot sequence X(N), N.ge.2*M, and the spline
C       coefficients C(N).
C
C Calling convention:
C ******************
C
C       SVIDER = SPLDER ( IDER, M, N, T, X, C, L, Q )
C
C Meaning of parameters:
C *********************
C
C       SPLDER  ( O )   Function or derivative value.
C       IDER    ( I )   Derivative order required, with 0.le.IDER
C                       and IDER.le.2*M. If IDER.eq.0, the function
C                       value is returned; otherwise, the IDER-th
C                       derivative of the spline is returned.
C       M       ( I )   Half order of the spline, with M.gt.0.
C       N       ( I )   Number of knots and spline coefficients,
C                       with N.ge.2*M.
C       T       ( I )   Argument at which the spline or its deri-
C                       vative is to be evaluated, with X(1).le.T
C                       and T.le.X(N).
C       X(N)    ( I )   Strictly increasing knot sequence array,
C                       X(I-1).lt.X(I), I=2,...,N.
C       C(N)    ( I )   Spline coefficients, as evaluated by
C                       subroutine GVCSPL.
C       L       (I/O)   L contains an integer such that:
C                       X(L).le.T and T.lt.X(L+1) if T is within
C                       the range X(1).le.T and T.lt.X(N). If
C                       T.lt.X(1), L is set to 0, and if T.ge.X(N),
C                       L is set to N. The search for L is facili-
C                       tated if L has approximately the right
C                       value on entry.
C       Q(2*M)  ( W )   Internal work array.
C
C Remark:
C ******
C
C       This subroutine is an adaptation of subroutine SPLDER of
C       the paper by Lyche et al. (1983). No checking is performed
C       on the validity of the input parameters.
C
C Reference:
C *********
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      REAL*8 FUNCTION SPLDER ( IDER, M, N, T, X, C, L, Q )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER ( ZERO=0D0, ONE=1D0 )
      DIMENSION X(N), C(N), Q(2*M)
C
C***  Derivatives of IDER.ge.2*M are alway zero
C
      M2 =  2 * M
      K  = M2 - IDER
      IF (K.LT.1) THEN
         SPLDER = ZERO
         RETURN
      ENDIF
C
C***  Search for the interval value L
C
      CALL SEARCH ( N, X, T, L )
C
C***  Initialize parameters and the 1st row of the B-spline
C***  coefficients tableau
C
      TT   = T
      MP1  =  M + 1
      NPM  =  N + M
      M2M1 = M2 - 1
      K1   =  K - 1
      NK   =  N - K
      LK   =  L - K
      LK1  = LK + 1
      LM   =  L - M
      JL   =  L + 1
      JU   =  L + M2
      II   =  N - M2
      ML   = -L
      DO 2 J=JL,JU
         IF ((J.GE.MP1).AND.(J.LE.NPM)) THEN
            Q(J+ML) = C(J-M)
         ELSE
            Q(J+ML) = ZERO
         ENDIF
    2 CONTINUE
C
C***  The following loop computes differences of the B-spline
C***  coefficients. If the value of the spline is required,
C***  differencing is not necessary.
C
      IF (IDER.GT.0) THEN
         JL = JL - M2
         ML = ML + M2
         DO 6 I=1,IDER
            JL = JL + 1
            II = II + 1
            J1 = MAX0(1,JL)
            J2 = MIN0(L,II)
            MI = M2 - I
            J  = J2 + 1
            IF (J1.LE.J2) THEN
               DO 3 JIN=J1,J2
                  J  =  J - 1
                  JM = ML + J
                  Q(JM) = (Q(JM) - Q(JM-1)) / (X(J+MI) - X(J))
    3          CONTINUE
            ENDIF
            IF (JL.GE.1) GO TO 6
               I1 =  I + 1
               J  = ML + 1
               IF (I1.LE.ML) THEN
                  DO 5 JIN=I1,ML
                     J    =  J - 1
                     Q(J) = -Q(J-1)
    5             CONTINUE
               ENDIF
    6    CONTINUE
         DO 7 J=1,K
            Q(J) = Q(J+IDER)
    7    CONTINUE
      ENDIF
C
C***  Compute lower half of the evaluation tableau
C
      IF (K1.GE.1) THEN      
         DO 14 I=1,K1
            NKI  =  NK + I
            IR   =   K
            JJ   =   L
            KI   =   K - I
            NKI1 = NKI + 1
C***        Right-hand B-splines
            IF (L.GE.NKI1) THEN
               DO 9 J=NKI1,L
                  Q(IR) = Q(IR-1) + (TT-X(JJ))*Q(IR)
                  JJ    = JJ - 1
                  IR    = IR - 1
    9          CONTINUE
            ENDIF
C***        Middle B-splines
            LK1I = LK1 + I
            J1 = MAX0(1,LK1I)
            J2 = MIN0(L, NKI)
            IF (J1.LE.J2) THEN
               DO 11 J=J1,J2
                  XJKI  = X(JJ+KI)
                  Z     = Q(IR)
                  Q(IR) = Z + (XJKI-TT)*(Q(IR-1)-Z)/(XJKI-X(JJ))
                  IR    = IR - 1
                  JJ    = JJ - 1
   11          CONTINUE
            ENDIF
C***        Left-hand B-splines
            IF (LK1I.LE.0) THEN
               JJ    = KI
               LK1I1 =  1 - LK1I
               DO 13 J=1,LK1I1
                  Q(IR) = Q(IR) + (X(JJ)-TT)*Q(IR-1)
                  JJ    = JJ - 1
                  IR    = IR - 1
   13          CONTINUE
            ENDIF
   14    CONTINUE
      ENDIF
C
C***  Compute the return value
C
      Z = Q(K)
C***  Multiply with factorial if IDER.gt.0
      IF (IDER.GT.0) THEN
         DO 16 J=K,M2M1
            Z = Z * J
   16    CONTINUE
      ENDIF
      SPLDER = Z
C
C***  Ready
C
      RETURN
      END
C SEARCH.FOR, 1985-06-03
C
C***********************************************************************
C
C SUBROUTINE SEARCH (REAL*8)
C
C Purpose:
C *******
C
C       Given a strictly increasing knot sequence X(1) < ... < X(N),
C       where N >= 1, and a real number T, this subroutine finds the
C       value L such that X(L) <= T < X(L+1).  If T < X(1), L = 0;
C       if X(N) <= T, L = N.
C
C Calling convention:
C ******************
C
C       CALL SEARCH ( N, X, T, L )
C
C Meaning of parameters:
C *********************
C
C       N       ( I )   Knot array dimensioning parameter.
C       X(N)    ( I )   Stricly increasing knot array.
C       T       ( I )   Input argument whose knot interval is to
C                       be found.
C       L       (I/O)   Knot interval parameter. The search procedure
C                       is facilitated if L has approximately the
C                       right value on entry.
C
C Remark:
C ******
C
C       This subroutine is an adaptation of subroutine SEARCH from
C       the paper by Lyche et al. (1983). No checking is performed
C       on the input parameters and data; the algorithm may fail if
C       the input sequence is not strictly increasing.
C
C Reference:
C *********
C
C       T. Lyche, L.L. Schumaker, & K. Sepehrnoori, Fortran subroutines
C       for computing smoothing and interpolating natural splines.
C       Advances in Engineering Software 5(1983)1, pp. 2-5.
C
C***********************************************************************
C
      SUBROUTINE SEARCH ( N, X, T, L )
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION X(N)
C
      IF (T.LT.X(1)) THEN
C***     Out of range to the left
         L = 0
         RETURN
      ENDIF
      IF (T.GE.X(N)) THEN
C***     Out of range to the right
         L = N
         RETURN
      ENDIF
C***  Validate input value of L
      L = MAX0(L,1)
      IF (L.GE.N) L = N-1
C
C***  Often L will be in an interval adjoining the interval found
C***  in a previous call to search
C
      IF (T.GE.X(L)) GO TO 5
      L = L - 1
      IF (T.GE.X(L)) RETURN
C
C***  Perform bisection
C
      IL = 1
    3 IU = L
    4 L = (IL+IU) / 2
      IF (IU-IL.LE.1) RETURN
      IF (T.LT.X(L)) GO TO 3
      IL = L
      GO TO 4
    5 IF (T.LT.X(L+1)) RETURN
      L = L + 1
      IF (T.LT.X(L+1)) RETURN
      IL = L + 1
      IU = N
      GO TO 4
C
      END