* 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 ++++++++++++++++++++++++