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