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
|
subroutine fpcoco(iopt,m,x,y,w,v,s,nest,maxtr,maxbin,n,t,c,sq,sx,
* bind,e,wrk,lwrk,iwrk,kwrk,ier)
c ..scalar arguments..
real*8 s,sq
integer iopt,m,nest,maxtr,maxbin,n,lwrk,kwrk,ier
c ..array arguments..
integer iwrk(kwrk)
real*8 x(m),y(m),w(m),v(m),t(nest),c(nest),sx(m),e(nest),wrk(lwrk)
*
logical bind(nest)
c ..local scalars..
integer i,ia,ib,ic,iq,iu,iz,izz,i1,j,k,l,l1,m1,nmax,nr,n4,n6,n8,
* ji,jib,jjb,jl,jr,ju,mb,nm
real*8 sql,sqmax,term,tj,xi,half
c ..subroutine references..
c fpcosp,fpbspl,fpadno,fpdeno,fpseno,fpfrno
c ..
c set constant
half = 0.5e0
c determine the maximal admissible number of knots.
nmax = m+4
c the initial choice of knots depends on the value of iopt.
c if iopt=0 the program starts with the minimal number of knots
c so that can be guarantied that the concavity/convexity constraints
c will be satisfied.
c if iopt = 1 the program will continue from the point on where she
c left at the foregoing call.
if(iopt.gt.0) go to 80
c find the minimal number of knots.
c a knot is located at the data point x(i), i=2,3,...m-1 if
c 1) v(i) ^= 0 and
c 2) v(i)*v(i-1) <= 0 or v(i)*v(i+1) <= 0.
m1 = m-1
n = 4
do 20 i=2,m1
if(v(i).eq.0. .or. (v(i)*v(i-1).gt.0. .and.
* v(i)*v(i+1).gt.0.)) go to 20
n = n+1
c test whether the required storage space exceeds the available one.
if(n+4.gt.nest) go to 200
t(n) = x(i)
20 continue
c find the position of the knots t(1),...t(4) and t(n-3),...t(n) which
c are needed for the b-spline representation of s(x).
do 30 i=1,4
t(i) = x(1)
n = n+1
t(n) = x(m)
30 continue
c test whether the minimum number of knots exceeds the maximum number.
if(n.gt.nmax) go to 210
c main loop for the different sets of knots.
c find corresponding values e(j) to the knots t(j+3),j=1,2,...n-6
c e(j) will take the value -1,1, or 0 according to the requirement
c that s(x) must be locally convex or concave at t(j+3) or that the
c sign of s''(x) is unrestricted at that point.
40 i= 1
xi = x(1)
j = 4
tj = t(4)
n6 = n-6
do 70 l=1,n6
50 if(xi.eq.tj) go to 60
i = i+1
xi = x(i)
go to 50
60 e(l) = v(i)
j = j+1
tj = t(j)
70 continue
c we partition the working space
nm = n+maxbin
mb = maxbin+1
ia = 1
ib = ia+4*n
ic = ib+nm*maxbin
iz = ic+n
izz = iz+n
iu = izz+n
iq = iu+maxbin
ji = 1
ju = ji+maxtr
jl = ju+maxtr
jr = jl+maxtr
jjb = jr+maxtr
jib = jjb+mb
c given the set of knots t(j),j=1,2,...n, find the least-squares cubic
c spline which satisfies the imposed concavity/convexity constraints.
call fpcosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,nm,mb,wrk(ia),
*
* wrk(ib),wrk(ic),wrk(iz),wrk(izz),wrk(iu),wrk(iq),iwrk(ji),
* iwrk(ju),iwrk(jl),iwrk(jr),iwrk(jjb),iwrk(jib),ier)
c if sq <= s or in case of abnormal exit from fpcosp, control is
c repassed to the driver program.
if(sq.le.s .or. ier.gt.0) go to 300
c calculate for each knot interval t(l-1) <= xi <= t(l) the
c sum((wi*(yi-s(xi)))**2).
c find the interval t(k-1) <= x <= t(k) for which this sum is maximal
c on the condition that this interval contains at least one interior
c data point x(nr) and that s(x) is not given there by a straight line.
80 sqmax = 0.
sql = 0.
l = 5
nr = 0
i1 = 1
n4 = n-4
do 110 i=1,m
term = (w(i)*(sx(i)-y(i)))**2
if(x(i).lt.t(l) .or. l.gt.n4) go to 100
term = term*half
sql = sql+term
if(i-i1.le.1 .or. (bind(l-4).and.bind(l-3))) go to 90
if(sql.le.sqmax) go to 90
k = l
sqmax = sql
nr = i1+(i-i1)/2
90 l = l+1
i1 = i
sql = 0.
100 sql = sql+term
110 continue
if(m-i1.le.1 .or. (bind(l-4).and.bind(l-3))) go to 120
if(sql.le.sqmax) go to 120
k = l
nr = i1+(m-i1)/2
c if no such interval is found, control is repassed to the driver
c program (ier = -1).
120 if(nr.eq.0) go to 190
c if s(x) is given by the same straight line in two succeeding knot
c intervals t(l-1) <= x <= t(l) and t(l) <= x <= t(l+1),delete t(l)
n8 = n-8
l1 = 0
if(n8.le.0) go to 150
do 140 i=1,n8
if(.not. (bind(i).and.bind(i+1).and.bind(i+2))) go to 140
l = i+4-l1
if(k.gt.l) k = k-1
n = n-1
l1 = l1+1
do 130 j=l,n
t(j) = t(j+1)
130 continue
140 continue
c test whether we cannot further increase the number of knots.
150 if(n.eq.nmax) go to 180
if(n.eq.nest) go to 170
c locate an additional knot at the point x(nr).
j = n
do 160 i=k,n
t(j+1) = t(j)
j = j-1
160 continue
t(k) = x(nr)
n = n+1
c restart the computations with the new set of knots.
go to 40
c error codes and messages.
170 ier = -3
go to 300
180 ier = -2
go to 300
190 ier = -1
go to 300
200 ier = 4
go to 300
210 ier = 5
300 return
end
|
