V3FIT
r8spline.f
1 C******************** START FILE SPLINE.FOR ; GROUP TRKRLIB ******************
2 CC
3 CCCCCCCCCCCCCCCC
4 C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:
5 C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR
6 C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.
7 C
8 C THE CODES (SPLEEN,SPLAAN & SPEVAL) ARE ADAPTATIONS
9 C BY R.M. WIELAND FOR SPECIAL CASES ... SEE COMMENTS
10 
11 
12  SUBROUTINE r8spline (N, X, Y, B, C, D)
13 !============
14 ! idecl: explicitize implicit INTEGER declarations:
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  INTEGER ip1,ii
18 !============
19 ! idecl: explicitize implicit REAL declarations:
20  real*8 zbomb
21 !============
22  INTEGER N
23  real*8 x(n), y(n), b(n), c(n), d(n)
24 C
25 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
26 C FOR A CUBIC INTERPOLATING SPLINE
27 C
28 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
29 C
30 C FOR X(I) .LE. X .LE. X(I+1)
31 C
32 C INPUT..
33 C
34 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
35 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
36 C Y = THE ORDINATES OF THE KNOTS
37 C
38 C OUTPUT..
39 C
40 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
41 C
42 C USING P TO DENOTE DIFFERENTIATION,
43 C
44 C Y(I) = S(X(I))
45 C B(I) = SP(X(I))
46 C C(I) = SPP(X(I))/2
47 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
48 C
49 CCCCCCCCCCCCCCC
50 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
51 C TO EVALUATE THE SPLINE.
52 C
53 C
54  INTEGER NM1, IB, I
55  real*8 t
56 C
57  nm1 = n-1
58  IF ( n .LT. 2 ) RETURN
59 C
60 C
61 C DMC - CHECK ORDINATES
62  DO 5 i=1,nm1
63  ip1=i+1
64  IF(x(ip1).LE.x(i)) then
65  write(6,9000)
66  write(6,9001) (x(ii),ii=1,n)
67  zbomb=sqrt(-1.0_r8/(x(1)*x(1)))
68  write(6,*) zbomb
69  write(6,*) ' spline.for: code will stop.'
70  call exit(1)
71  9000 FORMAT(/
72  >' ? UNORDERED ORDINATE ARRAY PASSED TO SPLINE SUBROUTINE:')
73  9001 FORMAT(1x,5(1x,1pe12.5))
74  endif
75  5 CONTINUE
76 C
77  IF ( n .LT. 3 ) GO TO 50
78 C
79 C SET UP TRIDIAGONAL SYSTEM
80 C
81 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
82 C
83  d(1) = x(2) - x(1)
84  c(2) = (y(2) - y(1))/d(1)
85  DO 10 i = 2, nm1
86  d(i) = x(i+1) - x(i)
87  b(i) = 2._r8*(d(i-1) + d(i))
88  c(i+1) = (y(i+1) - y(i))/d(i)
89  c(i) = c(i+1) - c(i)
90  10 CONTINUE
91 C
92 C END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N)
93 C OBTAINED FROM DIVIDED DIFFERENCES
94 C
95  b(1) = -d(1)
96  b(n) = -d(n-1)
97  c(1) = 0._r8
98  c(n) = 0._r8
99  IF ( n .EQ. 3 ) GO TO 15
100  c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1))
101  c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
102  c(1) = c(1)*d(1)**2/(x(4)-x(1))
103  c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
104 C
105 C FORWARD ELIMINATION
106 C
107  15 DO 20 i = 2, n
108  t = d(i-1)/b(i-1)
109  b(i) = b(i) - t*d(i-1)
110  c(i) = c(i) - t*c(i-1)
111  20 CONTINUE
112 C
113 C BACK SUBSTITUTION
114 C
115  c(n) = c(n)/b(n)
116  DO 30 ib = 1, nm1
117  i = n-ib
118  c(i) = (c(i) - d(i)*c(i+1))/b(i)
119  30 CONTINUE
120 C
121 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
122 C
123 C COMPUTE POLYNOMIAL COEFFICIENTS
124 C
125  b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2._r8*c(n))
126  DO 40 i = 1, nm1
127  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2._r8*c(i))
128  d(i) = (c(i+1) - c(i))/d(i)
129  c(i) = 3._r8*c(i)
130  40 CONTINUE
131  c(n) = 3._r8*c(n)
132  d(n) = d(n-1)
133  RETURN
134 C
135  50 b(1) = (y(2)-y(1))/(x(2)-x(1))
136  c(1) = 0._r8
137  d(1) = 0._r8
138  b(2) = b(1)
139  c(2) = 0._r8
140  d(2) = 0._r8
141  RETURN
142  END
143 C******************** END FILE SPLINE.FOR ; GROUP TRKRLIB ******************