V3FIT
spleen.f
1 C-------------------------------------------
2 C derived from Forsyth/Malcolm/Moler SPLINE.FOR -- different BC:
3 C y'=0 at RHS of data interval (high index side)
4 C
5 CCCCCCCCCCCCCCC
6  SUBROUTINE spleen (N, X, Y, B, C, D)
7  INTEGER N
8  REAL X(N), Y(N), B(N), C(N), D(N)
9 C
10 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
11 C FOR A CUBIC INTERPOLATING SPLINE FOR WHICH S'(XN)=0.
12 C
13 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
14 C
15 C FOR X(I) .LE. X .LE. X(I+1)
16 C
17 C INPUT..
18 C
19 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
20 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
21 C Y = THE ORDINATES OF THE KNOTS
22 C
23 C OUTPUT..
24 C
25 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
26 C
27 C USING P TO DENOTE DIFFERENTIATION,
28 C
29 C Y(I) = S(X(I))
30 C B(I) = SP(X(I))
31 C C(I) = SPP(X(I))/2
32 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
33 C
34 CCCCCCCCCCCCCCC
35 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
36 C TO EVALUATE THE SPLINE.
37 C
38 C
39  INTEGER NM1, IB, I
40  REAL T
41 C
42  nm1 = n-1
43  IF ( n .LT. 2 ) RETURN
44  IF ( n .LT. 3 ) GO TO 50
45 C
46 C SET UP TRIDIAGONAL SYSTEM
47 C
48 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
49 C
50  d(1) = x(2) - x(1)
51  c(2) = (y(2) - y(1))/d(1)
52  DO 10 i = 2, nm1
53  d(i) = x(i+1) - x(i)
54  b(i) = 2.*(d(i-1) + d(i))
55  c(i+1) = (y(i+1) - y(i))/d(i)
56  c(i) = c(i+1) - c(i)
57  10 CONTINUE
58 C
59 C END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N)
60 C OBTAINED FROM DIVIDED DIFFERENCES
61 C
62  b(1) = -d(1)
63  b(n) = +2.*d(n-1)
64  c(1) = 0.
65  c(n) = 0.
66  IF ( n .EQ. 3 ) GO TO 15
67  c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
68  c(1) = c(1)*d(1)**2/(x(4)-x(1))
69  c(n)= -(y(n)-y(n-1))/d(n-1)
70 C
71 C FORWARD ELIMINATION
72 C
73  15 DO 20 i = 2, n
74  t = d(i-1)/b(i-1)
75  b(i) = b(i) - t*d(i-1)
76  c(i) = c(i) - t*c(i-1)
77  20 CONTINUE
78 C
79 C BACK SUBSTITUTION
80 C
81  c(n) = c(n)/b(n)
82  DO 30 ib = 1, nm1
83  i = n-ib
84  c(i) = (c(i) - d(i)*c(i+1))/b(i)
85  30 CONTINUE
86 C
87 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
88 C
89 C COMPUTE POLYNOMIAL COEFFICIENTS
90 C
91  b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n))
92  DO 40 i = 1, nm1
93  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i))
94  d(i) = (c(i+1) - c(i))/d(i)
95  c(i) = 3.*c(i)
96  40 CONTINUE
97  c(n) = 3.*c(n)
98  d(n) = d(n-1)
99  RETURN
100 C
101  50 b(1) = (y(2)-y(1))/(x(2)-x(1))
102  c(1) = 0.
103  d(1) = 0.
104  b(2) = b(1)
105  c(2) = 0.
106  d(2) = 0.
107  RETURN
108  END