1
subroutine fpopsp(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,r,mr,r0,r1,dr,
2
* iopt,ider,tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpu,fpv,
4
c given the set of function values r(i,j) defined on the rectangular
5
c grid (u(i),v(j)),i=1,2,...,mu;j=1,2,...,mv, fpopsp determines a
6
c smooth bicubic spline approximation with given knots tu(i),i=1,..,nu
7
c in the u-direction and tv(j),j=1,2,...,nv in the v-direction. this
8
c spline sp(u,v) will be periodic in the variable v and will satisfy
9
c the following constraints
11
c s(tu(1),v) = dr(1) , tv(4) <=v<= tv(nv-3)
13
c s(tu(nu),v) = dr(4) , tv(4) <=v<= tv(nv-3)
15
c and (if iopt(2) = 1)
18
c ------------ = dr(2)*cos(v)+dr(3)*sin(v) , tv(4) <=v<= tv(nv-3)
21
c and (if iopt(3) = 1)
24
c ------------- = dr(5)*cos(v)+dr(6)*sin(v) , tv(4) <=v<= tv(nv-3)
27
c where the parameters dr(i) correspond to the derivative values at the
28
c poles as defined in subroutine spgrid.
30
c the b-spline coefficients of sp(u,v) are determined as the least-
31
c squares solution of an overdetermined linear system which depends
32
c on the value of p and on the values dr(i),i=1,...,6. the correspond-
33
c ing sum of squared residuals sq is a simple quadratic function in
34
c the variables dr(i). these may or may not be provided. the values
35
c dr(i) which are not given will be determined so as to minimize the
36
c resulting sum of squared residuals sq. in that case the user must
37
c provide some initial guess dr(i) and some estimate (dr(i)-step,
38
c dr(i)+step) of the range of possible values for these latter.
40
c sp(u,v) also depends on the parameter p (p>0) in such a way that
41
c - if p tends to infinity, sp(u,v) becomes the least-squares spline
42
c with given knots, satisfying the constraints.
43
c - if p tends to zero, sp(u,v) becomes the least-squares polynomial,
44
c satisfying the constraints.
45
c - the function f(p)=sumi=1,mu(sumj=1,mv((r(i,j)-sp(u(i),v(j)))**2)
46
c is continuous and strictly decreasing for p>0.
48
c ..scalar arguments..
49
integer ifsu,ifsv,ifbu,ifbv,mu,mv,mr,nu,nv,nuest,nvest,
53
integer ider(4),nru(mu),nrv(mv),iopt(3)
54
real*8 u(mu),v(mv),r(mr),dr(6),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv)
58
real*8 res,sq,sqq,sq0,sq1,step1,step2,three
59
integer i,id0,iop0,iop1,i1,j,l,lau,lav1,lav2,la0,la1,lbu,lbv,lb0,
60
* lb1,lc0,lc1,lcs,lq,lri,lsu,lsv,l1,l2,mm,mvnu,number
63
real*8 delta(6),drr(6),sum(6),a(6,6),g(6)
64
c ..function references..
66
c ..subroutine references..
71
c we partition the working space
75
mm = max0(nuest,mv+nvest)
77
mvnu = nuest*(mv+nvest-8)
90
c we calculate the smoothing spline sp(u,v) according to the input
91
c values dr(i),i=1,...,6.
96
call fpgrsp(ifsu,ifsv,ifbu,ifbv,0,u,mu,v,mv,r,mr,dr,
97
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,
98
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
99
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(la0),wrk(la1),wrk(lb0),
100
* wrk(lb1),wrk(lc0),wrk(lc1),wrk(lcs),nru,nrv)
103
if(id0.eq.0) sq0 = (r0-dr(1))**2
104
if(id1.eq.0) sq1 = (r1-dr(4))**2
106
c in case all derivative values dr(i) are given (step<=0) or in case
107
c we have spline interpolation, we accept this spline as a solution.
109
if(step(1).le.0. .and. step(2).le.0.) return
113
c number denotes the number of derivative values dr(i) that still must
114
c be optimized. let us denote these parameters by g(j),j=1,...,number.
116
if(id0.gt.0) go to 20
120
20 if(iop0.eq.0) go to 30
121
if(ider(2).ne.0) go to 30
122
step2 = step(1)*three/(tu(5)-tu(4))
125
delta(number+1) = step2
126
delta(number+2) = step2
128
30 if(id1.gt.0) go to 40
131
delta(number) = step(2)
132
40 if(iop1.eq.0) go to 50
133
if(ider(4).ne.0) go to 50
134
step2 = step(2)*three/(tu(nu)-tu(nu-4))
137
delta(number+1) = step2
138
delta(number+2) = step2
140
50 if(number.eq.0) return
141
c the sum of squared residulas sq is a quadratic polynomial in the
142
c parameters g(j). we determine the unknown coefficients of this
143
c polymomial by calculating (number+1)*(number+2)/2 different splines
144
c according to specific values for g(j).
149
call fpgrsp(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,r,mr,drr,
150
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sum(i),fp,fpu,fpv,mm,mvnu,
151
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
152
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(la0),wrk(la1),wrk(lb0),
153
* wrk(lb1),wrk(lc0),wrk(lc1),wrk(lcs),nru,nrv)
154
if(id0.eq.0) sq0 = (r0-drr(1))**2
155
if(id1.eq.0) sq1 = (r1-drr(4))**2
156
sum(i) = sum(i)+sq0+sq1
158
call fpgrsp(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,r,mr,drr,
159
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sqq,fp,fpu,fpv,mm,mvnu,
160
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
161
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(la0),wrk(la1),wrk(lb0),
162
* wrk(lb1),wrk(lc0),wrk(lc1),wrk(lcs),nru,nrv)
163
if(id0.eq.0) sq0 = (r0-drr(1))**2
164
if(id1.eq.0) sq1 = (r1-drr(4))**2
167
a(i,i) = (sum(i)+sqq-sq-sq)/step1**2
168
if(a(i,i).le.0.) go to 110
169
g(i) = (sqq-sum(i))/(step1+step1)
171
if(number.eq.1) go to 90
175
drr(l1) = dr(l1)+step1
180
drr(l2) = dr(l2)+step2
181
call fpgrsp(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,r,mr,drr,
182
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sqq,fp,fpu,fpv,mm,mvnu,
183
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
184
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(la0),wrk(la1),wrk(lb0),
185
* wrk(lb1),wrk(lc0),wrk(lc1),wrk(lcs),nru,nrv)
186
if(id0.eq.0) sq0 = (r0-drr(1))**2
187
if(id1.eq.0) sq1 = (r1-drr(4))**2
189
a(i,j) = (sq+sqq-sum(i)-sum(j))/(step1*step2)
194
c the optimal values g(j) are found as the solution of the system
195
c d (sq) / d (g(j)) = 0 , j=1,...,number.
196
90 call fpsysy(a,number,g)
201
c we determine the spline sp(u,v) according to the optimal values g(j).
202
110 call fpgrsp(ifsu,ifsv,ifbu,ifbv,0,u,mu,v,mv,r,mr,dr,
203
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,
204
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
205
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(la0),wrk(la1),wrk(lb0),
206
* wrk(lb1),wrk(lc0),wrk(lc1),wrk(lcs),nru,nrv)
207
if(id0.eq.0) sq0 = (r0-dr(1))**2
208
if(id1.eq.0) sq1 = (r1-dr(4))**2
added
removed
Lib/interpolate/fitpack/fpopsp.f
