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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
|
subroutine splopt ( tau, n, k, scrtch, t, iflag )
c from * a practical guide to splines * by c. de boor
calls bsplvb, banfac/slv
computes the knots t for the optimal recovery scheme of order k
c for data at tau(i), i=1,...,n .
c
c****** i n p u t ******
c tau.....array of length n , containing the interpolation points.
c a s s u m e d to be nondecreasing, with tau(i).lt.tau(i+k),all i.
c n.....number of data points .
c k.....order of the optimal recovery scheme to be used .
c
c****** w o r k a r r a y *****
c scrtch.....array of length (n-k)(2k+3) + 5k + 3 . the various
c contents are specified in the text below .
c
c****** o u t p u t ******
c iflag.....integer indicating success (=1) or failure (=2) .
c if iflag = 1, then
c t.....array of length n+k containing the optimal knots ready for
c use in optimal recovery. specifically, t(1) = ... = t(k) =
c tau(1) and t(n+1) = ... = t(n+k) = tau(n) , while the n-k
c interior knots t(k+1), ..., t(n) are calculated as described
c below under *method* .
c if iflag = 2, then
c k .lt. 3, or n .lt. k, or a certain linear system was found to
c be singular.
c
c****** p r i n t e d o u t p u t ******
c a comment will be printed in case ilfag = 2 or newton iterations
c failed to converge in n e w t m x iterations .
c
c****** m e t h o d ******
c the (interior) knots t(k+1), ..., t(n) are determined by newtons
c method in such a way that the signum function which changes sign at
c t(k+1), ..., t(n) and nowhere else in (tau(1),tau(n)) is orthogon-
c al to the spline space spline( k , tau ) on that interval .
c let xi(j) be the current guess for t(k+j), j=1,...,n-k. then
c the next newton iterate is of the form
c xi(j) + (-)**(n-k-j)*x(j) , j=1,...,n-k,
c with x the solution of the linear system
c c*x = d .
c here, c(i,j) = b(i)(xi(j)), all j, with b(i) the i-th b-spline of
c order k for the knot sequence tau , all i, and d is the vector
c given by d(i) = sum( -a(j) , j=i,...,n )*(tau(i+k)-tau(i))/k, all i,
c with a(i) = sum ( (-)**(n-k-j)*b(i,k+1,tau)(xi(j)) , j=1,...,n-k )
c for i=1,...,n-1, and a(n) = -.5 .
c (see chapter xiii of text and references there for a derivation)
c the first guess for t(k+j) is (tau(j+1)+...+tau(j+k-1))/(k-1) .
c iteration terminates if max(abs(x(j))) .lt. t o l , with
c t o l = t o l r t e *(tau(n)-tau(1))/(n-k) ,
c or else after n e w t m x iterations , currently,
c newtmx, tolrte / 10, .000001
c
integer iflag,k,n, i,id,index,j,km1,kpk,kpkm1,kpn,kp1,l,left
*,leftmk,lenw,ll,llmax,llmin,na,nb,nc,nd,newtmx,newton,nmk,nmkm1,nx
real scrtch(1),t(1),tau(n), del,delmax,floatk,sign,signst,sum
* ,tol,tolrte,xij
c dimension scrtch((n-k)*(2*k+3)+5*k+3), t(n+k)
current fortran standard makes it impossible to specify the precise dim-
c ensions of scrtch and t without the introduction of otherwise
c superfluous additional arguments .
data newtmx,tolrte / 10,.000001/
nmk = n-k
if (nmk) 1,56,2
1 print 601,n,k
601 format(13h argument n =,i4,29h in splopt is less than k =,i3)
go to 999
2 if (k .gt. 2) go to 3
print 602,k
602 format(13h argument k =,i3,27h in splopt is less than 3)
go to 999
3 nmkm1 = nmk - 1
floatk = k
kpk = k+k
kp1 = k+1
km1 = k-1
kpkm1 = kpk-1
kpn = k+n
signst = -1.
if (nmk .gt. (nmk/2)*2) signst = 1.
c scrtch(i) = tau-extended(i), i=1,...,n+k+k
nx = n+kpk+1
c scrtch(i+nx) = xi(i),i=0,...,n-k+1
na = nx + nmk + 1
c scrtch(i+na) = -a(i), i=1,...,n
nd = na + n
c scrtch(i+nd) = x(i) or d(i), i=1,...,n-k
nb = nd + nmk
c scrtch(i+nb) = biatx(i),i=1,...,k+1
nc = nb + kp1
c scrtch(i+(j-1)*(2k-1) + nc) = w(i,j) = c(i-k+j,j), i=j-k,...,j+k,
c j=1,...,n-k.
lenw = kpkm1*nmk
c extend tau to a knot sequence and store in scrtch.
do 5 j=1,k
scrtch(j) = tau(1)
5 scrtch(kpn+j) = tau(n)
do 6 j=1,n
6 scrtch(k+j) = tau(j)
c first guess for scrtch (.+nx) = xi .
scrtch(nx) = tau(1)
scrtch(nmk+1+nx) = tau(n)
do 10 j=1,nmk
sum = 0.
do 9 l=1,km1
9 sum = sum + tau(j+l)
10 scrtch(j+nx) = sum/float(km1)
c last entry of scrtch (.+na) = - a is always ...
scrtch(n+na) = .5
c start newton iteration.
newton = 1
tol = tolrte*(tau(n) - tau(1))/float(nmk)
c start newton step
compute the 2k-1 bands of the matrix c and store in scrtch(.+nc),
c and compute the vector scrtch(.+na) = -a.
20 do 21 i=1,lenw
21 scrtch(i+nc) = 0.
do 22 i=2,n
22 scrtch(i-1+na) = 0.
sign = signst
left = kp1
do 28 j=1,nmk
xij = scrtch(j+nx)
23 if (xij .lt. scrtch(left+1))go to 25
left = left + 1
if (left .lt. kpn) go to 23
left = left - 1
25 call bsplvb(scrtch,k,1,xij,left,scrtch(1+nb))
c the tau sequence in scrtch is preceded by k additional knots
c therefore, scrtch(ll+nb) now contains b(left-2k+ll)(xij)
c which is destined for c(left-2k+ll,j), and therefore for
c w(left-k-j+ll,j)= scrtch(left-k-j+ll + (j-1)*kpkm1 + nc)
c since we store the 2k-1 bands of c in the 2k-1 r o w s of
c the work array w, and w in turn is stored in s c r t c h ,
c with w(1,1) = scrtch(1 + nc) .
c also, c being of order n-k, we would want 1 .le.
c left-2k+ll .le. n-k or
c llmin = 2k-left .le. ll .le. n-left+k = llmax .
leftmk = left-k
index = leftmk-j + (j-1)*kpkm1 + nc
llmin = max0(1,k-leftmk)
llmax = min0(k,n-leftmk)
do 26 ll=llmin,llmax
26 scrtch(ll+index) = scrtch(ll+nb)
call bsplvb(scrtch,kp1,2,xij,left,scrtch(1+nb))
id = max0(0,leftmk-kp1)
llmin = 1 - min0(0,leftmk-kp1)
do 27 ll=llmin,kp1
id = id + 1
27 scrtch(id+na) = scrtch(id+na) - sign*scrtch(ll+nb)
28 sign = -sign
call banfac(scrtch(1+nc),kpkm1,nmk,km1,km1,iflag)
go to (45,44),iflag
44 print 644
644 format(32h c in splopt is not invertible)
return
compute scrtch (.+nd) = d from scrtch (.+na) = - a .
45 i=n
46 scrtch(i-1+na) = scrtch(i-1+na) + scrtch(i+na)
i = i-1
if (i .gt. 1) go to 46
do 49 i=1,nmk
49 scrtch(i+nd) = scrtch(i+na)*(tau(i+k)-tau(i))/floatk
compute scrtch (.+nd) = x .
call banslv(scrtch(1+nc),kpkm1,nmk,km1,km1,scrtch(1+nd))
compute scrtch (.+nd) = change in xi . modify, if necessary, to
c prevent new xi from moving more than 1/3 of the way to its
c neighbors. then add to xi to obtain new xi in scrtch(.+nx).
delmax = 0.
sign = signst
do 53 i=1,nmk
del = sign*scrtch(i+nd)
delmax = amax1(delmax,abs(del))
if (del .gt. 0.) go to 51
del = amax1(del,(scrtch(i-1+nx)-scrtch(i+nx))/3.)
go to 52
51 del = amin1(del,(scrtch(i+1+nx)-scrtch(i+nx))/3.)
52 sign = -sign
53 scrtch(i+nx) = scrtch(i+nx) + del
call it a day in case change in xi was small enough or too many
c steps were taken.
if (delmax .lt. tol) go to 54
newton = newton + 1
if (newton .le. newtmx) go to 20
print 653,newtmx
653 format(33h no convergence in splopt after,i3,14h newton steps.)
54 do 55 i=1,nmk
55 t(k+i) = scrtch(i+nx)
56 do 57 i=1,k
t(i) = tau(1)
57 t(n+i) = tau(n)
return
999 iflag = 2
return
end
|
