V3FIT
r8pspline.f
1 C PSPLINE -- modified from SPLINE.FOR -- dmc 3 Oct 1995
2 C
3 C THE CODES (SPLINE & SEVAL) ARE TAKEN FROM:
4 C FORSYTHE,MALCOLM AND MOLER, "COMPUTER METHODS FOR
5 C MATHEMATICAL COMPUTATIONS",PRENTICE-HALL, 1977.
6 c
7 c PSPLINE -- periodic spline -- adaptation of SPLINE.FOR by D. McCune
8 c October 1995: The function to be splined is taken as having periodic
9 C derivatives (it need not be periodic itself), so that
10 c the interpolating function satisfies
11 c y'(x(1))=y'(x(N))
12 c y''(x(1))=y''(x(N))
13 c this results in a fully specified system (no addl BC assumptions
14 c required); it is not a tridiagonal system but a modified NM1xNM1
15 c tridiagonal system with non-zero offdiagonal corner (1,NM1), (NM1,1)
16 c elements. It remains symmetric & diagonally dominant, and is
17 c solved by Gaussian elimination w/o pivoting. NM1=N-1.
18 c
19  SUBROUTINE r8pspline (N, X, Y, B, C, D, WK)
20 !============
21 ! idecl: explicitize implicit REAL declarations:
22  IMPLICIT NONE
23  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
24  real*8 q
25 !============
26  INTEGER N
27  real*8 x(n), y(n), b(n), c(n), d(n), wk(n)
28 C
29 C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
30 C FOR A periodic CUBIC INTERPOLATING SPLINE
31 C
32 C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
33 C
34 C FOR X(I) .LE. X .LE. X(I+1)
35 C
36 C WK(...) is used as a workspace during the Guassian elimination.
37 C it represents the rightmost column of the array
38 C
39 C note B(N),C(N),D(N) give coeffs for a periodic continuation into
40 C the next period, up to X(N)+(X(2)-X(1)) only.
41 C
42 C SEVAL can be used to evaluate the spline, but, it is up to the
43 C SEVAL caller to bring the argument into the normal range X(1) to X(N).
44 C
45 C INPUT..
46 C
47 C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
48 C includes a complete period of the data, the last point being
49 C a repeat of the first point.
50 C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
51 C X(N)-X(1) is the period of the periodic function represented.
52 C Y = THE ORDINATES OF THE KNOTS
53 C
54 C OUTPUT..
55 C
56 C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
57 C
58 C USING P TO DENOTE DIFFERENTIATION,
59 C
60 C Y(I) = S(X(I))
61 C B(I) = SP(X(I))
62 C C(I) = SPP(X(I))/2
63 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT)
64 C
65 CCCCCCCCCCCCCCC
66 C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED
67 C TO EVALUATE THE SPLINE.
68 C
69 C
70  INTEGER NM1, NM2, IB, I
71  real*8 t
72 C
73 C-------------------------------------
74 C
75  nm1 = n-1
76  nm2 = nm1-1
77 C
78  IF ( n .LT. 4 ) THEN
79  write(6,9901)
80  9901 format(/
81  > ' ?PSPLINE -- at least 4 pts required for periodic spline.')
82  return
83  ENDIF
84 C
85 C SET UP MODIFIED NM1 x NM1 TRIDIAGONAL SYSTEM:
86 C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
87 C WK(1:NM2) = rightmost column above diagonal
88 C WK(NM1) = lower left corner element
89 C
90  d(1) = x(2) - x(1)
91  d(nm1) = x(n) - x(nm1)
92  b(1) = 2.0_r8*(d(1)+d(nm1))
93  c(1) = (y(n) - y(nm1))/d(nm1)
94  c(2) = (y(2) - y(1))/d(1)
95  c(1) = c(2) - c(1)
96  wk(1) = d(nm1)
97  DO 10 i = 2, nm1
98  d(i) = x(i+1) - x(i)
99  b(i) = 2._r8*(d(i-1) + d(i))
100  c(i+1) = (y(i+1) - y(i))/d(i)
101  c(i) = c(i+1) - c(i)
102  wk(i) = 0.0_r8
103  10 CONTINUE
104  wk(nm2) = d(nm2)
105  wk(nm1) = d(nm1)
106 C
107 C END CONDITIONS -- implied by periodicity
108 C C(1)=C(N) B(1)=B(N) D(1)=D(N) -- no additional assumption needed.
109 C
110 C FORWARD ELIMINATION
111 C WK(1)--WK(NM2) represent the rightmost column above the
112 C diagonal; WK(NM1) represents the non-zero lower left corner element
113 C which in each step is moved one column to the right.
114 C
115  15 DO 20 i = 2, nm2
116  t = d(i-1)/b(i-1)
117  b(i) = b(i) - t*d(i-1)
118  c(i) = c(i) - t*c(i-1)
119  wk(i) = wk(i) - t*wk(i-1)
120  q = wk(nm1)/b(i-1)
121  wk(nm1) = -q*d(i-1)
122  b(nm1) = b(nm1) - q*wk(i-1)
123  c(nm1) = c(nm1) - q*c(i-1)
124  20 CONTINUE
125 C
126 C correct the (NM1,NM2) element
127 C
128  wk(nm1) = wk(nm1) + d(nm2)
129 C
130 C complete the forward elimination: now WK(NM1) and WK(NM2) are
131 C the off diagonal elements of the 2x2 at the lower right hand corner
132 C
133  t = wk(nm1)/b(nm2)
134  b(nm1) = b(nm1) - t*wk(nm2)
135  c(nm1) = c(nm1) - t*c(nm2)
136 C
137 C BACK SUBSTITUTION
138 C
139  c(nm1) = c(nm1)/b(nm1)
140  c(nm2) = (c(nm2) - wk(nm2)*c(nm1))/b(nm2)
141 C
142  DO 30 ib = 3, nm1
143  i = n-ib
144  c(i) = (c(i) - d(i)*c(i+1) - wk(i)*c(nm1))/b(i)
145  30 CONTINUE
146 C
147 C PERIODICITY:
148 C
149  c(n)=c(1)
150 C
151 C C(I) IS NOW THE SIGMA(I) OF THE TEXT
152 C
153 C COMPUTE POLYNOMIAL COEFFICIENTS
154 C
155  DO 40 i = 1, nm1
156  b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2._r8*c(i))
157  d(i) = (c(i+1) - c(i))/d(i)
158  c(i) = 3._r8*c(i)
159  40 CONTINUE
160 C periodic continuation...
161  b(n) = b(1)
162  c(n) = c(1)
163  d(n) = d(1)
164 C
165  RETURN
166  END