      SUBROUTINE BSTRP (TAU, GTAU, T, N, K, BCOEF, Q, IFLAG)
C
C   THIS ROUTINE PRODUCES THE B-SPLINE COEFF.S BCOEF OF THE PIECEWISE
C   POLYNOMIAL OF ORDER K WITH KNOTS T(I) (I=1,...,N+K) WHICH HAS THE
C   VALUE GTAU(I) AT TAU(I) FOR I=1,...,N.
C
C******  I N P U T  ******
C  TAU.....ARRAY OF LENGTH  N , CONTAINING DATA POINT ABSCISSAE.
C    A S S U M P T I O N . . .  TAU  IS STRICTLY INCREASING
C  GTAU.....CORRESPONDING ARRAY OF LENGTH  N , CONTAINING DATA POINT
C        ORDINATES.
C  T.....KNOT SEQUENCE, OF LENGTH  N+K
C  N.....NUMBER OF DATA POINTS AND DIMENSION OF SPLINE SPACE  S(K,T)
C  K.....ORDER OF THE PIECEWISE POLYNOMIAL
C  IFLAG.....ON AN INITIAL CALL TO THE ROUTINE, IFLAG MAY BE ASSIGNED
C        ANY VALUE EXCEPT 0. THE ROUTINE MAY BE RECALLED WHEN ONLY GTAU
C        IS MODIFIED. IFLAG=0 WHEN THIS IS DONE.
C
C******  O U T P U T  ******
C  BCOEF.....THE B-COEFFICIENTS OF THE INTERPOLANT, OF LENGTH  N
C  Q.....ARRAY OF SIZE  (2*K-1)*N , CONTAINING THE TRIANGULAR FACTORIZ-
C        ATION OF THE COEFFICIENT MATRIX OF THE LINEAR SYSTEM FOR THE B-
C        COEFFICIENTS OF THE SPLINE INTERPOLANT.
C  IFLAG.....AN INTEGER INDICATING SUCCESS (= 0)  OR FAILURE (= 1)
C        THE LINEAR SYSTEM TO BE SOLVED IS (THEORETICALLY) INVERTIBLE IF
C        AND ONLY IF
C              B(I)(TAU(I)) .NE. 0   FOR ALL I.
C        VIOLATION OF THIS CONDITION IS CERTAIN TO LEAD TO IFLAG = 1.
C
C******  M E T H O D  ******
C     THE I-TH EQUATION OF THE LINEAR SYSTEM  A*BCOEF = B  FOR THE B-CO-
C  EFFS OF THE INTERPOLANT ENFORCES INTERPOLATION AT  TAU(I), I=1,...,N.
C  HENCE,  B(I) = GTAU(I), ALL I, AND  A  IS A BAND MATRIX WITH  2K-1
C  BANDS (IF IT IS INVERTIBLE).
C     THE MATRIX  A  IS GENERATED ROW BY ROW AND STORED, DIAGONAL BY DI-
C  AGONAL, IN THE  R O W S  OF THE ARRAY  Q , WITH THE MAIN DIAGONAL GO-
C  ING INTO ROW  K .  SEE COMMENTS IN THE PROGRAM BELOW.
C     THE BANDED SYSTEM IS THEN SOLVED BY A CALL TO BANFAC (WHICH CON-
C  STRUCTS THE TRIANGULAR FACTORIZATION FOR  A  AND STORES IT IN  Q),
C  FOLLOWED BY A CALL TO BANSLV (WHICH THEN OBTAINS THE SOLUTION BCOEF
C  BY SUBSTITUTION).
C     BANFAC PERFORMS NO PIVOTING SINCE THE TOTAL POSITIVITY OF THE
C  MATRIX A MAKES THIS UNNECESSARY.
C
      INTEGER IFLAG,K,N,   I,ILP1MX,J,JJ,KM1,KPKM2,LEFT,LENQ,NP1,NPK
      REAL BCOEF(N),GTAU(N),Q(1),T(1),TAU(N),TAUI
      KM1 = K - 1
      IF (IFLAG .EQ. 0) GO TO 40
      NP1 = N + 1
      NPK = N + K
      KPKM2 = 2*KM1
      LEFT = K
C                ZERO OUT ALL ENTRIES OF Q
      LENQ = N*(K+KM1)
      DO 5 I=1,LENQ
    5    Q(I) = 0.
C
C  ***   LOOP OVER I TO CONSTRUCT THE  N  INTERPOLATION EQUATIONS
      DO 30 I=1,N
         TAUI = TAU(I)
         ILP1MX = MIN0(I+K,NP1)
C        *** FIND  LEFT  IN THE CLOSED INTERVAL (I,I+K-1) SUCH THAT
C                T(LEFT) .LE. TAU(I) .LT. T(LEFT+1)
C        MATRIX IS SINGULAR IF THIS IS NOT POSSIBLE
         LEFT = MAX0(LEFT,I)
         IF (TAUI .LT. T(LEFT))         GO TO 998
   15       IF (TAUI .LT. T(LEFT+1))    GO TO 16
            LEFT = LEFT + 1
            IF (LEFT .LT. ILP1MX)       GO TO 15
         IF (LEFT .EQ. I+K)             GO TO 998
         IF (I .LT. N)                  GO TO 998
         IF (TAUI .GT. T(LEFT))         GO TO 998
         LEFT = N
         IF (T(LEFT) .GE. T(LEFT+1))    GO TO 998
C        *** THE I-TH EQUATION ENFORCES INTERPOLATION AT TAUI, HENCE
C        A(I,J) = B(J,K,T)(TAUI), ALL J. ONLY THE  K  ENTRIES WITH  J =
C        LEFT-K+1,...,LEFT ACTUALLY MIGHT BE NONZERO. THESE  K  NUMBERS
C        ARE RETURNED, IN  BCOEF (USED FOR TEMP.STORAGE HERE), BY THE
C        FOLLOWING
   16    CALL BSPEV(T,NPK,TAUI,LEFT,0,K,BCOEF,IFLAG)
C        LET Q DENOTE A TWO-DIMENSIONAL ARRAY OF DIMENSION (2*K-1,N).
C        WE THEREFORE WANT  BCOEF(J) = B(LEFT-K+J)(TAUI) TO GO INTO
C        A(I,LEFT-K+J), I.E., INTO  Q(I-(LEFT+J)+2*K,(LEFT+J)-K) SINCE
C        A(I+J,J)  IS TO GO INTO  Q(I+K,J), ALL I,J.  IN THE CURRENT
C        ROUTINE WE TREAT Q AS AN EQUIVALENT ONE-DIMENSIONAL ARRAY.
C        THUS WE WANT BCOEF(J) TO BE INSERTED INTO ENTRY
C            I -(LEFT+J) + 2*K + ((LEFT+J) - K-1)*(2*K-1)
C                   =  I-LEFT+1 + (LEFT -K)*(2*K-1) + (2*K-2)*J
C        OF  Q .
         JJ = I-LEFT+1 + (LEFT-K)*(K+KM1)
         DO 30 J=1,K
            JJ = JJ+KPKM2
   30       Q(JJ) = BCOEF(J)
C
C     ***OBTAIN FACTORIZATION OF  A  , STORED AGAIN IN  Q.
      CALL BANFAC ( Q, K+KM1, N, KM1, KM1, IFLAG )
      IFLAG = IFLAG - 1
      IF (IFLAG .NE. 0) RETURN
C     *** SOLVE  A*BCOEF = GTAU  BY BACKSUBSTITUTION
   40 DO 41 I=1,N
   41    BCOEF(I) = GTAU(I)
      CALL BANSLV ( Q, K+KM1, N, KM1, KM1, BCOEF )
      RETURN
C
C     *** ERROR RETURN
  998 IFLAG = 1
      RETURN
      END
