1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
|
subroutine setupq ( x, dx, y, npoint, v, qty )
c from * a practical guide to splines * by c. de boor
c to be called in s m o o t h
c put delx = x(.+1) - x(.) into v(.,4),
c put the three bands of q-transp*d into v(.,1-3), and
c put the three bands of (d*q)-transp*(d*q) at and above the diagonal
c into v(.,5-7) .
c here, q is the tridiagonal matrix of order (npoint-2,npoint)
c with general row 1/delx(i) , -1/delx(i) - 1/delx(i+1) , 1/delx(i+1)
c and d is the diagonal matrix with general row dx(i) .
integer npoint, i,npm1
real dx(npoint),qty(npoint),v(npoint,7),x(npoint),y(npoint),
* diff, prev
npm1 = npoint - 1
v(1,4) = x(2) - x(1)
do 11 i=2,npm1
v(i,4) = x(i+1) - x(i)
v(i,1) = dx(i-1)/v(i-1,4)
v(i,2) = - dx(i)/v(i,4) - dx(i)/v(i-1,4)
11 v(i,3) = dx(i+1)/v(i,4)
v(npoint,1) = 0.
do 12 i=2,npm1
12 v(i,5) = v(i,1)**2 + v(i,2)**2 + v(i,3)**2
if (npm1 .lt. 3) go to 14
do 13 i=3,npm1
13 v(i-1,6) = v(i-1,2)*v(i,1) + v(i-1,3)*v(i,2)
14 v(npm1,6) = 0.
if (npm1 .lt. 4) go to 16
do 15 i=4,npm1
15 v(i-2,7) = v(i-2,3)*v(i,1)
16 v(npm1-1,7) = 0.
v(npm1,7) = 0.
construct q-transp. * y in qty.
prev = (y(2) - y(1))/v(1,4)
do 21 i=2,npm1
diff = (y(i+1)-y(i))/v(i,4)
qty(i) = diff - prev
21 prev = diff
return
end
|
