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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include "bspln.h"
define ILLEGAL_ORDER 2 #error returns
define NO_DEG_FREEDOM 3
.help splsqv 2 "IRAF Math Library"
Adapted from L2APPR, * A Practical Guide To Splines * by C. DeBoor.
Eliminated data entry via the common /data/ (D.Tody, 4-july-82). Calls
subprograms BSPLVB, BCHFAC, BCHSLV.
SPLSQV constructs the (weighted discrete) l2-approximation by splines of
order K with knot sequence T[1], ..., t[nt+k] to given data points
(TAU[i], GTAU[i]), i=1,...,ntau. The b-spline coefficients BCOEF of the
approximating spline are determined from the normal equations using
Cholesky'smethod.
SPLSQV (tau, gtau, weight, ntau, t, nt, k, work, diag, bcoef, ier)
INPUT -----
ntau Number of data points
tau[ntau] x-value of the data points.
gtau[ntau] y-value of the data points.
weight[ntau] The corresponding weights.
t[nt] The knot sequence
nt The dimension of the space of splines of order k with knots t.
k The order of the b-spline to be fitted.
WORK ARRAYS -----
work[k,nt] A work array of size (at least) K*NT. its first K rows are used
for the K lower diagonals of the gramian matrix C.
diag[nt] A work array of length NT used in BCHFAC.
OUTPUT -----
bcoef[nt] The b-spline coefficients of the l2-approximation.
ier Error code: zero if ok, else small integer error code.
METHOD -----
The b-spline coefficients of the l2-appr. are determined as the solution
of the normal equations
sum ((B[i],B[j]) * BCOEF[j] : j=1:nt) = (B[i],G), i = 1 to nt.
where, B[i] denotes the i-th b-spline, G denotes the function to be
approximated, and the INNER PRODUCT of two functions F and G is given by
(F,G) := sum (F(TAU[i]) * G(TAU[i]) * WEIGHT[i] : i=1:ntau).
The relevant function values of the b-splines B[i], i=1:nt, are supplied
by the subprogram BSPLVB. The coeff. matrix C, with
C(i,j) := (B[i], B[j]), i,j=1:n,
of the normal equations is symmetric and (2*k-1)-banded, therefore can be
specified by giving its k bands at or below the diagonal. for i=1,...,n,
we store
(B[i],B[j]) = B[i,j] in WORK[i-j+1,j], j=i,...,min0(i+k-1,nt)
and the right side
(B[i], G) in BCOEF[i].
since b-spline values are most efficiently generated by finding simultaneously
the value of EVERY nonzero b-spline at one point, the entries of C (i.e., of
WORK), are generated by computing, for each ll, all the terms involving
TAU(ll) simultaneously and adding them to all relevant entries.
.endhelp _______________________________________________________________
procedure splsqv (tau, gtau, weight, ntau, t, nt, k, work, diag, bcoef, ier)
real tau[ntau], gtau[ntau], weight[ntau]
real t[nt], work[k,nt], diag[nt], bcoef[nt]
int ntau, nt, k, ier
int i, j, jj, left, leftmk, ll, mm
real biatx[KMAX], dw
begin
ier = 0
if (k < 1 || k > KMAX) {
ier = ILLEGAL_ORDER
return
}
if (nt <= k || nt >= ntau+k) {
ier = NO_DEG_FREEDOM
return
}
do j = 1, nt {
bcoef[j] = 0.
do i = 1, k
work[i,j] = 0.
}
left = k
leftmk = 0
do ll = 1, ntau {
while (left != nt) {
if (tau[ll] < t[left+1])
break
left = left + 1
leftmk = leftmk + 1
}
call bsplvb (t, k, 1, tau[ll], left, biatx)
# BIATX[mm] contains the value of B[left-k+mm] at TAU[ll].
# hence, with DW := BIATX[mm] * WEIGHT[ll], the number DW * GTAU[ll]
# is a summand in the inner product
# (B[left-k+mm],G) which goes into BCOEF[left-k+mm]
# and the number BIATX[jj] * dw is a summand in the inner product
# (B[left-k+jj], B[left-k+mm]), into WORK[jj-mm+1,left-k+mm]
# since (left-k+jj) - (left-k+mm) + 1 = jj - mm + 1.
do mm = 1, k {
dw = biatx[mm] * weight[ll]
j = leftmk + mm
bcoef[j] = dw * gtau[ll] + bcoef[j]
i = 1
do jj = mm, k {
work[i,j] = biatx[jj] * dw + work[i,j]
i = i + 1
}
}
}
# construct cholesky factorization for C in WORK, then use it to solve
# the normal equations
# c * x = bcoef
# for X, and store X in BCOEF.
call bchfac (work, k, nt, diag)
call bchslv (work, k, nt, bcoef)
return
end
|
