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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
|
subroutine chol1d ( p, v, qty, npoint, ncol, u, qu )
c from * a practical guide to splines * by c. de boor
c from * a practical guide to splines * by c. de boor
c to be called in s m o o t h
constructs the upper three diags. in v(i,j), i=2,,npoint-1, j=1,3, of
c the matrix 6*(1-p)*q-transp.*(d**2)*q + p*r, then computes its
c l*l-transp. decomposition and stores it also in v, then applies
c forward and backsubstitution to the right side q-transp.*y in qty
c to obtain the solution in u .
integer ncol,npoint, i,npm1,npm2
real p,qty(npoint),qu(npoint),u(npoint),v(npoint,7), prev,ratio
* ,six1mp,twop
npm1 = npoint - 1
c construct 6*(1-p)*q-transp.*(d**2)*q + p*r
six1mp = 6.*(1.-p)
twop = 2.*p
do 2 i=2,npm1
v(i,1) = six1mp*v(i,5) + twop*(v(i-1,4)+v(i,4))
v(i,2) = six1mp*v(i,6) + p*v(i,4)
2 v(i,3) = six1mp*v(i,7)
npm2 = npoint - 2
if (npm2 .ge. 2) go to 10
u(1) = 0.
u(2) = qty(2)/v(2,1)
u(3) = 0.
go to 41
c factorization
10 do 20 i=2,npm2
ratio = v(i,2)/v(i,1)
v(i+1,1) = v(i+1,1) - ratio*v(i,2)
v(i+1,2) = v(i+1,2) - ratio*v(i,3)
v(i,2) = ratio
ratio = v(i,3)/v(i,1)
v(i+2,1) = v(i+2,1) - ratio*v(i,3)
20 v(i,3) = ratio
c
c forward substitution
u(1) = 0.
v(1,3) = 0.
u(2) = qty(2)
do 30 i=2,npm2
30 u(i+1) = qty(i+1) - v(i,2)*u(i) - v(i-1,3)*u(i-1)
c back substitution
u(npoint) = 0.
u(npm1) = u(npm1)/v(npm1,1)
i = npm2
40 u(i) = u(i)/v(i,1)-u(i+1)*v(i,2)-u(i+2)*v(i,3)
i = i - 1
if (i .gt. 1) go to 40
c construct q*u
41 prev = 0.
do 50 i=2,npoint
qu(i) = (u(i) - u(i-1))/v(i-1,4)
qu(i-1) = qu(i) - prev
50 prev = qu(i)
qu(npoint) = -qu(npoint)
return
end
|
