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