* spline.f - interpolation routines by various authors
*        all routines in this file are for unlimited distribution

*-----------------------------------------------------------------------
*                      CUBIC SPLINE INTERPOLATION
* adapted by Alfred Morris (NSWC), from a program by Carl de Boor (U Wisconsin)
*-----------------------------------------------------------------------
      SUBROUTINE CBSPL (X, Y, A, B, C, N, IBEG, IEND, ALPHA, BETA, IERR)
      REAL X(N), Y(N), A(N), B(N), C(N)

      IF (N < 3) GO TO 200
 
*-----------------------------------------------------------------------
*     A tridiagonal linear system for the unknown slopes S(I) of
*     F at X(I), I=1,...,N, is generated and then solved by Gauss
*     elimination, with S(I) ending up in A(I) for all I.  A, B, C
*     are used initially for work spaces.
*-----------------------------------------------------------------------
      DO 10 M = 2,N
         B(M) = X(M) - X(M-1)
         IF (B(M) .LE. 0.0) GO TO 210
         C(M) = (Y(M) - Y(M-1))/B(M)
  10  CONTINUE
      IERR = 0

*   Construct the first equation from the boundary condition, of the form
*             C(1)*S(1) + B(1)*S(2) = A(1)
      IF (IBEG - 1) 20,30,40

*     No condition at left end.
  20  C(1) = B(3)
      B(1) = X(3) - X(1)
      A(1) = ((B(2) + 2.0*B(1))*B(3)*C(2) + B(2)*B(2)*C(3))/B(1)
      GO TO 50

*     Slope prescribed at left end.
  30  C(1) = 1.0
      B(1) = 0.0
      A(1) = ALPHA
      GO TO 50

*     Second derivative prescribed at left end.
  40  C(1) = 2.0
      B(1) = 1.0
      A(1) = 3.0*C(2) - 0.5*ALPHA*B(2)

*-----------------------------------------------------------------------
*     For the interior knots, generate the corresponding equations and
*     carry out the forward pass of Gauss elimination, after which the
*     M-th equation reads  C(M)*S(M) + B(M)*S(M+1) = A(M).
*-----------------------------------------------------------------------
  50  NM1 = N - 1
      DO 51 M = 2,NM1
         T = -B(M+1)/C(M-1)
         A(M) = T*A(M-1) + 3.0*(B(M)*C(M+1) + B(M+1)*C(M))
         C(M) = T*B(M-1) + 2.0*(B(M) + B(M+1))
  51  CONTINUE

*-----------------------------------------------------------------------
*     If the slope at the right end is given, then set A(N) to the
*     slope and go to back substitution.  Otherwise, construct the
*     last equation from the second boundary condition, of the form
*               R*S(N-1) + C(N)*S(N) = A(N)
*-----------------------------------------------------------------------
      IF (IEND - 1) 60, 80, 90
  60  IF (N .EQ. 3 .AND. IBEG .EQ. 0) GO TO 70

*-----------------------------------------------------------------------
*  No condition at the right end.  Either N .GE. 4 or
*  there is a condition at the left end.
*-----------------------------------------------------------------------
      R = X(N) - X(N-2)
      DEL  = (Y(NM1) - Y(N-2))/B(NM1)
      A(N) = ((B(N) + 2.0*R)*B(NM1)*C(N) + B(N)*B(N)*DEL)/R
      C(N) = B(NM1)
      GO TO 100

*-----------------------------------------------------------------------
*     No conditions at the end points and N = 3.  In this case,
*     the second boundary condition does not provide us with a 
*     new equation.  For convenience, we use the following...
*-----------------------------------------------------------------------
  70  A(N) = 2.0*C(N)
      C(N) = 1.0
      R = 1.0
      GO TO 100

*     Slope prescribed at right end.
  80  A(N) = BETA
      GO TO 110

*     Second derivative prescribed at right end.
  90  A(N) = 3.0*C(N) + 0.5*BETA*B(N)
      C(N) = 2.0
      R = 1.0

*     Complete forward pass of Gauss elimination.
 100  T = -R/C(NM1)
      A(N) = (T*A(NM1) + A(N))/(T*B(NM1) + C(N))

*     Carry out back substitution.
 110  DO 120 I = 1,NM1
         J = N - I
         A(J) = (A(J) - B(J)*A(J+1))/C(J)
 120  CONTINUE

*     Generate the cubic coefficients B(I) and C(I).
      DO 130 I = 1,NM1
         H = B(I+1)
         DEL = (Y(I+1) - Y(I))/H
         T = A(I) + A(I+1) - 2.0*DEL
         B(I) = (DEL - A(I) - T)/H
         C(I) = (T/H)/H
 130  CONTINUE
      RETURN

*     error return
 200  IERR = 1
      RETURN
  210 IERR = 2
      RETURN
      END ! of CBSPL


*-----------------------------------------------------------------------
*                       CUBIC SPLINE EVALUATION
*  adapted by Alfred Morris (NSWC), from a program by Rondall Jones (Sandia NL)
*-----------------------------------------------------------------------
*         SCOMP evaluates a cubic spline at the abscissas in XI.
*         It is assumed that the coefficients of the polynomials
*         which form the spline are provided.
*-----------------------------------------------------------------------
*
*     DESCRIPTION OF ARGUMENTS
*
*       --INPUT--
*
*         X   - array of the first N abscissas (in increasing order)
*               that define the spline.
*         Y   - array of the first N ordinates that define the spline.
*         A,B,C arrays that contain the coefficients of the polynomials
*               which form the spline.  If I = 1,...,N then the spline
*               has the value
*                    Y(I) + A(I)*DX + B(I)*DX**2 + C(I)*DX**3
*               for X(I) .LE. XX .LE. X(I+1).  Here DX = XX-X(I).
*         N   - the number of polynomials that define the spline.
*               the arrays X, Y, A, B, C must be dimensioned at 
*               least N.  N must be greater than or equal to 1.
*         XI  - the abscissa or array of abscissas (in arbitrary order)
*               at which the spline is to be evaluated.
*         NI  - the number of abscissas at which the spline is to be
*               evaluated.  If NI is greater than 1 then XI and YI
*               must be arrays of dimension NI or larger.
*               It is assumed that NI is greater than or equal to 1.
*
*       --OUTPUT--
*
*         YI  - array of values of the spline (ordinates) at XI.
*         IERR- status code
*               0  the spline was evaluated at each abscissa in XI.
*               1  input error - NI is not positive.
*
*-----------------------------------------------------------------------
      SUBROUTINE SCOMP (X,Y,A,B,C,N,XI,YI,NI,IERR)
      REAL X(N), Y(N), A(N), B(N), C(N), XI(NI), YI(NI)

*     check input
      IF (NI > 0) GO TO 1
         IERR = 1
         RETURN
   1  IERR = 0

*     K is index on value of XI being worked on.  XX is that value.
*     I is current index into X array.
      K  = 1
      XX = XI(1)
      IF (XX < X(1)) GO TO 90
      IF (XX .GE. X(N)) GO TO 80
      IL = 1
      IR = N

*     bisection search
  10  I = (IL + IR)/2
      IF (I .EQ. IL) GO TO 100
      IF (XX - X(I)) 20,100,30
  20  IR = I
      GO TO 10
  30  IL = I
      GO TO 10

*     linear forward search
  50  IF (XX < X(I+1)) GO TO 100
      I = I + 1
      GO TO 50

*     XX is greater than X(N) or less than X(1)
  80  I = N
      GO TO 100
  90  I = 1

*     evaluation
 100  DX = XX - X(I)
      YI(K) = Y(I) + DX*(A(I) + DX*(B(I) + DX*C(I)))

*     next point
      IF (K .GE. NI) RETURN
      K = K + 1
      XX = XI(K)
      IF (XX < X(1)) GO TO 90
      IF (XX .GE. X(N)) GO TO 80
      IF (XX - XI(K-1)) 110, 100, 50
 110  IL = 1
      IR = MIN0(I+1,N)
      GO TO 10
      END ! of SCOMP


*-------------------------------------------------------------------------
*                       CUBIC  SPLINE  EVALUTATION
* adapted by Alfred Morris (NSWC), from a program by Rondall Jones (Sandia NL)
*-------------------------------------------------------------------------
*
*     ABSTRACT
*
*         SEVAL evaluates a cubic spline and its first
*         derivative at the abscissas in XI.  It is assumed that
*         the coefficients of the polynomials which form the spline
*         are provided.
*
*     DESCRIPTION OF ARGUMENTS
*
*       --INPUT--
*
*         X   - array of the first N abscissas (in increasing order)
*               that define the spline.
*         Y   - array of the first N ordinates that define the spline.
*         A,B,C arrays that contain the coefficients of the polynomials
*               which form the spline.  If I = 1,...,N  then the spline
*               has the value
*                    Y(I) + A(I)*DX + B(I)*DX**2 + C(I)*DX**3
*               for X(I) .LE. XX .LE. X(I+1).  Here DX = XX-X(I).
*         N   - the number of polynomials that define the spline.
*               the arrays X, Y, A, C, C  must be dimensioned at least N.
*               N must be greater than or equal to 1.
*         XI  - the abscissa or array of abscissas (in arbitrary order)
*               at which the spline is to be evaluated.
*         NI  - the number of abscissas at which the spline is to be
*               evaluated.  If NI is greater than 1, then XI, YI,
*               and YPI must be arrays dimensioned at least NI.
*               NI must be greater than or equal to 1.
*
*       --OUTPUT--
*
*         YI  - array of values of the spline (ordinates) at XI.
*         YPI - array of values of the first derivative of spline at XI.
*         IERR- status code
*               0  the spline was evaluated at each abscissa in XI.
*               1  input error - NI is not positive.
*
*-------------------------------------------------------------------------
      SUBROUTINE SEVAL (X,Y,A,B,C,N,XI,YI,YPI,NI,IERR)

      REAL X(N),Y(N),A(N),B(N),C(N),XI(NI),YI(NI),YPI(NI)

*     check input
      IF (NI .GT. 0) GO TO 1
         IERR = 1
         RETURN
  1   IERR = 0

*     K is index on value of XI being worked on.  XX is that value.
*     I is current index into X array.
      K  = 1
      XX = XI(1)
      IF (XX .LT. X(1)) GO TO 90
      IF (XX .GE. X(N)) GO TO 80
      IL = 1
      IR = N

*     bisection search
  10  I = (IL + IR)/2
      IF (I .EQ. IL) GO TO 100
      IF (XX - X(I)) 20,100,30
  20  IR = I
      GO TO 10
  30  IL = I
      GO TO 10

*     linear forward search
  50  IF (XX .LT. X(I+1)) GO TO 100
      I = I + 1
      GO TO 50

*     XX is greater than X(N) or less than X(1)
  80  I = N
      GO TO 100
  90  I = 1

*     evaluation
 100  DX = XX - X(I)
      YI(K) = Y(I) + DX*(A(I) + DX*(B(I) + DX*C(I)))
      BI = B(I) + B(I)
      CI = 3.0*C(I)
      YPI(K) = A(I) + DX*(BI + DX*CI)

*     next point
      IF (K .GE. NI) RETURN
      K = K + 1
      XX = XI(K)
      IF (XX .LT. X(1)) GO TO 90
      IF (XX .GE. X(N)) GO TO 80
      IF (XX - XI(K-1)) 110,100,50
 110  IL = 1
      IR = MIN0(I+1,N)
      GO TO 10
      END


*    The following routines are from the Naval Surface Warfare Center library
*      by Alfred Morris unless otherwise noted.  For unlimited distribution.

*-----------------------------------------------------------------------
*            WEIGHTED LEAST SQUARES CUBIC SPLINE FITTING
*-----------------------------------------------------------------------
      SUBROUTINE SPFIT (X, Y, WGT, M, BREAK, L, Z, A, B, C, WK, IERR)

      REAL X(M), Y(M), WGT(M), BREAK(L)
      REAL Z(*), A(*), B(*), C(*), WK(*)
      REAL TEMP(20)
*---------------------
*     REAL Z(L-1), A(L-1), B(L-1), C(L-1), WK(7*L + 18)
*---------------------
      IF (L .LT. 2) GO TO 100
      N = L + 2

*                define the knots for the B-splines
      WK(1) = BREAK(1)
      WK(2) = BREAK(1)
      WK(3) = BREAK(1)
      WK(4) = BREAK(1)
      DO 10 J = 2,L
         IF (BREAK(J - 1) .GE. BREAK(J)) GO TO 110
         WK(J + 3) = BREAK(J)
  10  CONTINUE
      WK(L + 4) = BREAK(L)
      WK(L + 5) = BREAK(L)
      WK(L + 6) = BREAK(L)

*     obtain the B-spline coefficients of the least squares fit
      LA = N + 5
      LW = LA + N
      LQ = LW + N
      CALL BSLSQ (X, Y, WGT, M, WK(1), N, 4, WK(LA),
     &            WK(LW), WK(LQ), IERR)
      IF (IERR .LT. 0) GO TO 120
      IERR = 0

*     obtain the coefficients of the fit in Taylor series form
      CALL BSPP (WK(1), WK(LA), N, 4, BREAK,
     &             WK(LQ), LM1, TEMP)
      K = LQ
      DO 20 J = 1,LM1
         Z(J) = WK(K)
         A(J) = WK(K + 1)
         B(J) = WK(K + 2)
         C(J) = WK(K + 3)
         K = K + 4
  20  CONTINUE
      RETURN

*                       error return
 100  IERR = 1
      RETURN
 110  IERR = 2
      RETURN
 120  IERR = 3
      RETURN
      END


*-----------------------------------------------------------------------
*              CONVERSION FROM B-SPLINE REPRESENTATION
*              TO PIECEWISE POLYNOMIAL REPRESENTATION
*
*
*     INPUT ...
*
*       T     knot sequence of length N+K
*       A     B-spline coefficient sequence of length n
*       N     length of A
*       K     order of the B-splines
*
*     OUTPUT ...
*
*       BREAK breakpoint sequence, of length L+1, containing
*             (in increasing order) the distinct points of the
*             sequence T(K),...,T(N+1).
*       C     KXL matrix where C(I,J) = (I-1)st right derivative
*             of the PP at BREAK(J) divided by factorial(I-1).
*       L     number of polynomials which form the PP
*
*     WORK AREA ...
*
*       WK    2-dimensional array of dimension (K,K+1)
*
*-----------------------------------------------------------------------
      SUBROUTINE BSPP (T, A, N, K, BREAK, C, L, WK)

      REAL T(*), A(N), BREAK(*), C(K,*), WK(K,*)

      L = 0
      BREAK(1) = T(K)
      IF (K .EQ. 1) GO TO 100
      KM1 = K - 1
      KP1 = K + 1

*          general K-th order case
      DO 60 LEFT = K,N
         IF (T(LEFT) .EQ. T(LEFT + 1)) GO TO 60
         L = L + 1
         BREAK(L + 1) = T(LEFT + 1)
         DO 10 J = 1,K
            JJ = LEFT - K + J
            WK(J,1) = A(JJ)
  10     CONTINUE

         DO 21 J = 1,KM1
            JP1 = J + 1
            KMJ = K - J
            DO 20 I = 1,KMJ
               IL = I + LEFT
               ILKJ = IL - KMJ
               DIFF = T(IL) - T(ILKJ)
               WK(I,JP1) = (WK(I+1,J) - WK(I,J))/DIFF
  20        CONTINUE
  21     CONTINUE

         WK(1,KP1) = 1.0
         X = T(LEFT)
         C(K,L) = WK(1,K)
         R = 1.0
         DO 50 J = 1,KM1
            JP1 = J + 1
            S = 0.0
            DO 30 I = 1,J
               IL = I + LEFT
               ILJ = IL - J
               TERM = WK(I,KP1)/(T(IL) - T(ILJ))
               WK(I,KP1) = S + (T(IL) - X)*TERM
               S = (X - T(ILJ))*TERM
  30        CONTINUE
            WK(JP1,KP1) = S

            S = 0.0
            KMJ = K - J
            DO 40 I = 1,JP1
               S = S + WK(I,KMJ)*WK(I,KP1)
  40        CONTINUE
            R = (R*FLOAT(KMJ))/FLOAT(J)
            C(KMJ,L) = R*S
  50     CONTINUE
  60  CONTINUE
      RETURN

*          piecewise constant case
 100  DO 110 LEFT = K,N
         IF (T(LEFT) .EQ. T(LEFT + 1)) GO TO 110
         L = L + 1
         BREAK(L + 1) = T(LEFT + 1)
         C(1,L) = A(LEFT)
 110  CONTINUE
      RETURN
      END


*-----------------------------------------------------------------------
*
*        BSLSQ PRODUCES THE B-SPLINE COEFFICIENTS OF A PIECEWISE
*              POLYNOMIAL P(X) OF ORDER K WHICH MINIMIZES
*
*                SUM (WGT(J)*(P(TAU(J)) - GTAU(J))**2).
*
*
*     INPUT ...
*
*       TAU   array of length ntau containing data point abscissae.
*       GTAU  array of length ntau containing data point ordinates.
*       WGT   array of length ntau containing the weights.
*       NTAU  number of data points to be fitted.
*       T     knot sequence of length N + K.
*       N     dimension of the piecewise polynomial space.
*       K     order of the B-splines.
*
*     OUTPUT ...
*
*       A     array of length N containing the B-spline coefficients
*             of the L2 approximation.
*
*       IERR  integer reporting the status of the results ...
*
*             0  the coefficient matrix is nonsingular.  the
*                unique least squares solution was obtained.
*             1  the coefficient matrix is singular.  a
*                least squares solution was obtained.
*            -1  input errors were detected.
*
*-----------------------------------------------------------------------
      SUBROUTINE BSLSQ (TAU, GTAU, WGT, NTAU, T, N, K, A, WK, Q, IERR)
      REAL TAU(NTAU), GTAU(NTAU), WGT(NTAU)
      REAL T(*), A(N), WK(N), Q(K,N)

      IF (NTAU .LT. MAX0(2,K)) GO TO 100
      IF (TAU(1) .LT. T(K) .OR. TAU(NTAU) .GT. T(N + 1)) GO TO 100

      DO 10 I = 2,NTAU
         IF (TAU(I - 1) .GT. TAU(I)) GO TO 100
  10  CONTINUE

      DO 21 J = 1,N
         A(J) = 0.0
         DO 20 I = 1,K
            Q(I,J) = 0.0
  20     CONTINUE
  21  CONTINUE

      LEFT = K
      DO 70 L = 1,NTAU

*        *** find the index left such that
*            T(LEFT) .LE. TAU(L) .LT. T(LEFT+1)
  30     IF (LEFT .EQ. N) GO TO 40
            IF (TAU(L) .LT. T(LEFT+1)) GO TO 40
            LEFT = LEFT + 1
            GO TO 30

  40     JJ = 0
         CALL BSPVB (T, K, K, JJ, TAU(L), LEFT, WK)

         LEFTMK = LEFT - K
         DO 61 MM = 1,K
            DW = WK(MM)*WGT(L)
            J = LEFTMK + MM
            A(J) = DW*GTAU(L) + A(J)
            I = 1
            DO 60 JJ = MM,K
               Q(I,J) = WK(JJ)*DW + Q(I,J)
               I = I + 1
  60        CONTINUE
  61     CONTINUE
  70  CONTINUE

*        solve the normal equations
      CALL BCHFAC (Q, K, N, WK, IERR)
      CALL BCHSLV (Q, K, N, A)
      RETURN

*             error return
  100 IERR = -1
      RETURN
      END


*-----------------------------------------------------------------------
*
*    BSPVB CALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES
*     AT X OF ORDER MAX(JHIGH,J + 1) WHERE T(K) .LE. X .LT. T(N+1).
*
*     DESCRIPTION OF ARGUMENTS
*
*         INPUT
*
*          T       - knot vector of length N + K.
*          K       - highest possible order of the B-splines.
*          JHIGH   - ORDER OF B-SPLINES (1 .LE. JHIGH .LE. K).
*          JHIGH   - order of B-splines (1 .LE. JHIGH .LE. K).
*          J       - J .LE. 0  GIVES B-SPLINES OF ORDER JHIGH.
*          J       - J .LE. 0  gives B-splines of order JHIGH.
*                    J .GE. 1  on a previous call to BSPVB the B-splines of 
*                              order J were computed and stored in BLIST.
*                              It is assumed that work has not been modified
*                              and that J .LT. K.
*          X       - argument of the B-splines.          
*          LEFT    - largest integer such that
*                    T(LEFT) .LE. X .LT. T(LEFT+1)
*
*         OUTPUT
*
*          BLIST   - vector of length K for spline values.
*          J       - B-splines of order J have been computed
*                    and stored in BLIST.
*
*-----------------------------------------------------------------------
*      written by Carl de Boor (University of Wisconsin) and modified
*         by A.H. Morris (NSWC).
*-----------------------------------------------------------------------
      SUBROUTINE BSPVB (T, K, JHIGH, J, X, LEFT, BLIST)

      REAL T(*), BLIST(K)

      IF (J .GT. 0) GO TO 10
         J = 1
         BLIST(1) = 1.0
         IF (J .GE. JHIGH) RETURN

  10  S = 0.0
      DO 20 L = 1,J
         I = LEFT + L
         IMJ = I - J
         TIMJ = T(IMJ)
         TI = T(I)
         TERM = BLIST(L)/(TI - TIMJ)
         BLIST(L) = S + (TI - X)*TERM
         S = (X - TIMJ)*TERM
  20  CONTINUE
      J = J + 1
      BLIST(J) = S
      IF (J .LT. JHIGH) GO TO 10

      RETURN
      END


*-----------------------------------------------------------------------
*  from  * A PRACTICAL GUIDE TO SPLINES *  by C. de Boor
*  constructs Cholesky factorization
*                     C  =  L * D * L-TRANSPOSE
*  with L unit lower triangular and D diagonal, for given matrix C of
*  order N, in case C is (symmetric) positive semidefinite
*  and banded, having NB diagonals at and below the main diagonal.
*
*******  INPUT  ******
*
*     N      the order of the matrix C.
*
*     NB     the bandwidth of C, i.e.,
*               C(I,J) = 0 for ABS(I-J) .GT. NB .
*
*     W      work array of size NB by N containing the NB diagonals
*            in its rows, with the main diagonal in row 1.  Precisely,
*            W(I,J) contains C(I+J-1,J), I=1,...,NB, J=1,...,N.
*            for example, the interesting entries of a seven diagonal
*            symmetric matrix C of order 9 would be stored in W as
*
*                       11 22 33 44 55 66 77 88 99
*                       21 32 43 54 65 76 87 98
*                       31 42 53 64 75 86 97
*                       41 52 63 74 85 96
*
*            All other entries of W not identified with an entry of C
*            are never referenced.          
*
*     DIAG   work array of length N.
*
*******  O U T P U T  ******
*                                                         T
*     W      contains the Cholesky factorization C = L*D*L   where
*            W(1,I) = 1/D(I,I) and W(I,J) = L(I-1+J,J) (I=2,...,NB).
*
*     IFLAG  0 if C is nonsingular and 1 if C is singular.
*
******  M E T H O D  ******
*
*   Gauss elimination, adapted to the symmetry and bandedness of C, is used.
*   Near zero pivots are handled in a special way.  The diagonal element
*     C(K,K) = W(1,K)  is saved initially in DIAG(K), all K.  At the K-th
*  elimination step, the current pivot element, viz. W(1,K), is compared
*  with its original value, DIAG(K).  If, as the result of prior
*  elimination steps, this element has been reduced by about a word
*  length, (i.e., if W(1,K)+DIAG(K) .LE. DIAG(K)), then the pivot is declared
*  to be zero, and the entire K-th row is declared to be linearly 
*  dependent on the preceding rows.  This has the effect of producing
*   X(K) = 0  when solving  C*X = B  for  X, regardless of  B.  Justification
*  for this is as follows.  In contemplated applications of this program,
*  the given equations are the normal equations for some least-squares
*  approximation problem, DIAG(K) = C(K,K) gives the norm-square 
*  of the K-th basis function, and, at this point,  W(1,K)  contains the
*  norm-square of the error in the least-squares approximation to the 
*  K-th basis function by linear combinations of the first K-1.  Having
*  W(1,K)+DIAG(K) .LE. DIAG(K) signifies that the K-th function is linearly
*  dependent to machine accuracy on the first K-1 functions, therefore
*  can safely be left out from the basis of approximating functions.
*   The solution of a linear system
*                       C*X = B
*   is effected by the succession of the following  T W O  calls ...
*     CALL BCHFAC (W, NB, N, DIAG, IFLAG)   , to get factorization
*     CALL BCHSLV (W, NB, N, B, X )            , to solve for X.
*-----------------------------------------------------------------------
      SUBROUTINE BCHFAC (W, NB, N, DIAG, IFLAG)

      REAL W(NB,N), DIAG(N)

      IF (N .GT. 1) GO TO 10
         IFLAG = 1
         IF (W(1,1) .EQ. 0.0) RETURN
         IFLAG = 0
         W(1,1) = 1.0/W(1,1)
         RETURN

*     store the diagonal of C in DIAG
  10  DO 11 K = 1,N
         DIAG(K) = W(1,K)
  11  CONTINUE

*     factorization
      IFLAG = 0
      DO 60 K = 1,N
         T = W(1,K) + DIAG(K)
         IF (T .NE. DIAG(K)) GO TO 30
            IFLAG = 1
            DO 20 J = 1,NB
               W(J,K) = 0.0
  20        CONTINUE
            GO TO 60

  30     T = 1.0/W(1,K)
         W(1,K) = T
         IMAX = MIN0(NB - 1,N - K)
         IF (IMAX .LT. 1) GO TO 60
         JMAX = IMAX
         DO 50 I = 1,IMAX
            RATIO = T*W(I+1,K)
            KPI = K + I
            DO 40 J = 1,JMAX
               IPJ = I + J
               W(J,KPI) = W(J,KPI) - W(IPJ,K)*RATIO
  40        CONTINUE
            JMAX = JMAX - 1
            W(I+1,K) = RATIO
  50     CONTINUE
  60  CONTINUE
      RETURN
      END


*-----------------------------------------------------------------------
*
*     BCHSLV solves the linear system C*X = B for X when W contains
*     the Cholesky factorization obtained by the subroutine BCHFAC
*     for the banded symmetric positive definite matrix C.
*
*     INPUT ...
*
*        N   the order of the matrix C
*        NB  the bandwidth of C
*        W   the Cholesky factorization of C
*        B   vector of length N containing the right side
*
*     OUTPUT ...
*
*        B   solution X of the linear system C*X = B
*
*                                       T
*     NOTE.  The factorization C = L*D*L  is used, where L is a
*     unit lower triangular matrix and D a diagonal matrix.
*
*-----------------------------------------------------------------------
      SUBROUTINE BCHSLV (W, NB, N, B)

      REAL W(NB,N), B(N)

      IF (N .GT. 1) GO TO 10
      B(1) = B(1)*W(1,1)
      RETURN

*     forward substitution.  Solve L*Y = B for Y and store Y in B.
  10  NBM1 = NB - 1
      DO 30 K = 1,N
         JMAX = MIN0(NBM1,N - K)
         IF (JMAX .LT. 1) GO TO 30
         DO 20 J = 1,JMAX
            JPK = J + K
            B(JPK) = B(JPK) - W(J + 1,K)*B(K)
  20     CONTINUE
  30  CONTINUE

*                              T     -1
*     backsubstitution. Solve L X = D  Y  for X and store X in B.
      K = N
  40     B(K) = B(K)*W(1,K)
         JMAX = MIN0(NBM1,N - K)
         IF (JMAX .LT. 1) GO TO 60
         DO 50 J = 1,JMAX
            JPK = J + K
            B(K) = B(K) - W(J + 1,K)*B(JPK)
  50     CONTINUE
  60     K = K - 1
         IF (K .GT. 0) GO TO 40
      RETURN
      END

*+++++++++++++++++++++++++ End of file spline.f ++++++++++++++++++++++++
