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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
subroutine fpintb(t,n,bint,nk1,x,y)
c subroutine fpintb calculates integrals of the normalized b-splines
c nj,k+1(x) of degree k, defined on the set of knots t(j),j=1,2,...n.
c it makes use of the formulae of gaffney for the calculation of
c indefinite integrals of b-splines.
c
c calling sequence:
c call fpintb(t,n,bint,nk1,x,y)
c
c input parameters:
c t : real array,length n, containing the position of the knots.
c n : integer value, giving the number of knots.
c nk1 : integer value, giving the number of b-splines of degree k,
c defined on the set of knots ,i.e. nk1 = n-k-1.
c x,y : real values, containing the end points of the integration
c interval.
c output parameter:
c bint : array,length nk1, containing the integrals of the b-splines.
c ..
c ..scalars arguments..
integer n,nk1
real*8 x,y
c ..array arguments..
real*8 t(n),bint(nk1)
c ..local scalars..
integer i,ia,ib,it,j,j1,k,k1,l,li,lj,lk,l0,min
real*8 a,ak,arg,b,f,one
c ..local arrays..
real*8 aint(6),h(6),h1(6)
c initialization.
one = 0.1d+01
k1 = n-nk1
ak = k1
k = k1-1
do 10 i=1,nk1
bint(i) = 0.0d0
10 continue
c the integration limits are arranged in increasing order.
a = x
b = y
min = 0
if(a-b) 30,160,20
20 a = y
b = x
min = 1
30 if(a.lt.t(k1)) a = t(k1)
if(b.gt.t(nk1+1)) b = t(nk1+1)
c using the expression of gaffney for the indefinite integral of a
c b-spline we find that
c bint(j) = (t(j+k+1)-t(j))*(res(j,b)-res(j,a))/(k+1)
c where for t(l) <= x < t(l+1)
c res(j,x) = 0, j=1,2,...,l-k-1
c = 1, j=l+1,l+2,...,nk1
c = aint(j+k-l+1), j=l-k,l-k+1,...,l
c = sumi((x-t(j+i))*nj+i,k+1-i(x)/(t(j+k+1)-t(j+i)))
c i=0,1,...,k
l = k1
l0 = l+1
c set arg = a.
arg = a
do 90 it=1,2
c search for the knot interval t(l) <= arg < t(l+1).
40 if(arg.lt.t(l0) .or. l.eq.nk1) go to 50
l = l0
l0 = l+1
go to 40
c calculation of aint(j), j=1,2,...,k+1.
c initialization.
50 do 55 j=1,k1
aint(j) = 0.0d0
55 continue
aint(1) = (arg-t(l))/(t(l+1)-t(l))
h1(1) = one
do 70 j=1,k
c evaluation of the non-zero b-splines of degree j at arg,i.e.
c h(i+1) = nl-j+i,j(arg), i=0,1,...,j.
h(1) = 0.0d0
do 60 i=1,j
li = l+i
lj = li-j
f = h1(i)/(t(li)-t(lj))
h(i) = h(i)+f*(t(li)-arg)
h(i+1) = f*(arg-t(lj))
60 continue
c updating of the integrals aint.
j1 = j+1
do 70 i=1,j1
li = l+i
lj = li-j1
aint(i) = aint(i)+h(i)*(arg-t(lj))/(t(li)-t(lj))
h1(i) = h(i)
70 continue
if(it.eq.2) go to 100
c updating of the integrals bint
lk = l-k
ia = lk
do 80 i=1,k1
bint(lk) = -aint(i)
lk = lk+1
80 continue
c set arg = b.
arg = b
90 continue
c updating of the integrals bint.
100 lk = l-k
ib = lk-1
do 110 i=1,k1
bint(lk) = bint(lk)+aint(i)
lk = lk+1
110 continue
if(ib.lt.ia) go to 130
do 120 i=ia,ib
bint(i) = bint(i)+one
120 continue
c the scaling factors are taken into account.
130 f = one/ak
do 140 i=1,nk1
j = i+k1
bint(i) = bint(i)*(t(j)-t(i))*f
140 continue
c the order of the integration limits is taken into account.
if(min.eq.0) go to 160
do 150 i=1,nk1
bint(i) = -bint(i)
150 continue
160 return
end
|
