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
|
subroutine splev(t,n,c,k,x,y,m,ier)
c subroutine splev evaluates in a number of points x(i),i=1,2,...,m
c a spline s(x) of degree k, given in its b-spline representation.
c
c calling sequence:
c call splev(t,n,c,k,x,y,m,ier)
c
c input parameters:
c t : array,length n, which contains the position of the knots.
c n : integer, giving the total number of knots of s(x).
c c : array,length n, which contains the b-spline coefficients.
c k : integer, giving the degree of s(x).
c x : array,length m, which contains the points where s(x) must
c be evaluated.
c m : integer, giving the number of points where s(x) must be
c evaluated.
c
c output parameter:
c y : array,length m, giving the value of s(x) at the different
c points.
c ier : error flag
c ier = 0 : normal return
c ier =10 : invalid input data (see restrictions)
c
c restrictions:
c m >= 1
c-- t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
c
c other subroutines required: fpbspl.
c
c references :
c de boor c : on calculating with b-splines, j. approximation theory
c 6 (1972) 50-62.
c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
c applics 10 (1972) 134-149.
c dierckx p. : curve and surface fitting with splines, monographs on
c numerical analysis, oxford university press, 1993.
c
c author :
c p.dierckx
c dept. computer science, k.u.leuven
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c latest update : march 1987
c
c++ pearu: 11 aug 2003
c++ - disabled cliping x values to interval [min(t),max(t)]
c++ - removed the restriction of the orderness of x values
c++ - fixed initialization of sp to double precision value
c
c ..scalar arguments..
integer n,k,m,ier
c ..array arguments..
real*8 t(n),c(n),x(m),y(m)
c ..local scalars..
integer i,j,k1,l,ll,l1,nk1
c++..
integer k2
c..++
real*8 arg,sp,tb,te
c ..local array..
real*8 h(6)
c ..
c before starting computations a data check is made. if the input data
c are invalid control is immediately repassed to the calling program.
ier = 10
c-- if(m-1) 100,30,10
c++..
if(m.lt.1) go to 100
c..++
c-- 10 do 20 i=2,m
c-- if(x(i).lt.x(i-1)) go to 100
c-- 20 continue
30 ier = 0
c fetch tb and te, the boundaries of the approximation interval.
k1 = k+1
c++..
k2 = k1+1
c..++
nk1 = n-k1
tb = t(k1)
te = t(nk1+1)
l = k1
l1 = l+1
c main loop for the different points.
do 80 i=1,m
c fetch a new x-value arg.
arg = x(i)
c-- if(arg.lt.tb) arg = tb
c-- if(arg.gt.te) arg = te
c search for knot interval t(l) <= arg < t(l+1)
c++..
35 if(arg.ge.t(l) .or. l1.eq.k2) go to 40
l1 = l
l = l-1
go to 35
c..++
40 if(arg.lt.t(l1) .or. l.eq.nk1) go to 50
l = l1
l1 = l+1
go to 40
c evaluate the non-zero b-splines at arg.
50 call fpbspl(t,n,k,arg,l,h)
c find the value of s(x) at x=arg.
sp = 0.0d0
ll = l-k1
do 60 j=1,k1
ll = ll+1
sp = sp+c(ll)*h(j)
60 continue
y(i) = sp
80 continue
100 return
end
|
