1
subroutine fpopdi(ifsu,ifsv,ifbu,ifbv,u,mu,v,mv,z,mz,z0,dz,
2
* iopt,ider,tu,nu,tv,nv,nuest,nvest,p,step,c,nc,fp,fpu,fpv,
4
c given the set of function values z(i,j) defined on the rectangular
5
c grid (u(i),v(j)),i=1,2,...,mu;j=1,2,...,mv, fpopdi 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) = dz(1) , tv(4) <=v<= tv(nv-3)
13
c and (if iopt(2) = 1)
16
c ------------ = dz(2)*cos(v)+dz(3)*sin(v) , tv(4) <=v<= tv(nv-3)
19
c and (if iopt(3) = 1)
21
c s(tu(nu),v) = 0 tv(4) <=v<= tv(nv-3)
23
c where the parameters dz(i) correspond to the derivative values g(i,j)
24
c as defined in subroutine pogrid.
26
c the b-spline coefficients of sp(u,v) are determined as the least-
27
c squares solution of an overdetermined linear system which depends
28
c on the value of p and on the values dz(i),i=1,2,3. the correspond-
29
c ing sum of squared residuals sq is a simple quadratic function in
30
c the variables dz(i). these may or may not be provided. the values
31
c dz(i) which are not given will be determined so as to minimize the
32
c resulting sum of squared residuals sq. in that case the user must
33
c provide some initial guess dz(i) and some estimate (dz(i)-step,
34
c dz(i)+step) of the range of possible values for these latter.
36
c sp(u,v) also depends on the parameter p (p>0) in such a way that
37
c - if p tends to infinity, sp(u,v) becomes the least-squares spline
38
c with given knots, satisfying the constraints.
39
c - if p tends to zero, sp(u,v) becomes the least-squares polynomial,
40
c satisfying the constraints.
41
c - the function f(p)=sumi=1,mu(sumj=1,mv((z(i,j)-sp(u(i),v(j)))**2)
42
c is continuous and strictly decreasing for p>0.
44
c ..scalar arguments..
45
integer ifsu,ifsv,ifbu,ifbv,mu,mv,mz,nu,nv,nuest,nvest,
49
integer ider(2),nru(mu),nrv(mv),iopt(3)
50
real*8 u(mu),v(mv),z(mz),dz(3),tu(nu),tv(nv),c(nc),fpu(nu),fpv(nv)
54
real*8 res,sq,sqq,step1,step2,three
55
integer i,id0,iop0,iop1,i1,j,l,laa,lau,lav1,lav2,lbb,lbu,lbv,
56
* lcc,lcs,lq,lri,lsu,lsv,l1,l2,mm,mvnu,number
59
real*8 delta(3),dzz(3),sum(3),a(6,6),g(6)
60
c ..function references..
62
c ..subroutine references..
67
c we partition the working space
71
mm = max0(nuest,mv+nvest)
73
mvnu = nuest*(mv+nvest-8)
83
c we calculate the smoothing spline sp(u,v) according to the input
84
c values dz(i),i=1,2,3.
87
call fpgrdi(ifsu,ifsv,ifbu,ifbv,0,u,mu,v,mv,z,mz,dz,
88
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,
89
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
90
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(laa),wrk(lbb),
91
* wrk(lcc),wrk(lcs),nru,nrv)
97
c in case all derivative values dz(i) are given (step<=0) or in case
98
c we have spline interpolation, we accept this spline as a solution.
99
5 if(step.le.0. .or. sq.le.0.) return
103
c number denotes the number of derivative values dz(i) that still must
104
c be optimized. let us denote these parameters by g(j),j=1,...,number.
106
if(id0.gt.0) go to 10
110
10 if(iop0.eq.0) go to 20
111
if(ider(2).ne.0) go to 20
112
step2 = step*three/tu(5)
115
delta(number+1) = step2
116
delta(number+2) = step2
118
20 if(number.eq.0) return
119
c the sum of squared residuals sq is a quadratic polynomial in the
120
c parameters g(j). we determine the unknown coefficients of this
121
c polymomial by calculating (number+1)*(number+2)/2 different splines
122
c according to specific values for g(j).
127
call fpgrdi(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,z,mz,dzz,
128
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sum(i),fp,fpu,fpv,mm,mvnu,
129
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
130
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(laa),wrk(lbb),
131
* wrk(lcc),wrk(lcs),nru,nrv)
132
if(id0.eq.0) sum(i) = sum(i)+(z0-dzz(1))**2
134
call fpgrdi(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,z,mz,dzz,
135
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sqq,fp,fpu,fpv,mm,mvnu,
136
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
137
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(laa),wrk(lbb),
138
* wrk(lcc),wrk(lcs),nru,nrv)
139
if(id0.eq.0) sqq = sqq+(z0-dzz(1))**2
140
a(i,i) = (sum(i)+sqq-sq-sq)/step1**2
141
if(a(i,i).le.0.) go to 80
142
g(i) = (sqq-sum(i))/(step1+step1)
145
if(number.eq.1) go to 60
149
dzz(l1) = dz(l1)+step1
154
dzz(l2) = dz(l2)+step2
155
call fpgrdi(ifsu,ifsv,ifbu,ifbv,1,u,mu,v,mv,z,mz,dzz,
156
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sqq,fp,fpu,fpv,mm,mvnu,
157
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
158
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(laa),wrk(lbb),
159
* wrk(lcc),wrk(lcs),nru,nrv)
160
if(id0.eq.0) sqq = sqq+(z0-dzz(1))**2
161
a(i,j) = (sq+sqq-sum(i)-sum(j))/(step1*step2)
166
c the optimal values g(j) are found as the solution of the system
167
c d (sq) / d (g(j)) = 0 , j=1,...,number.
168
60 call fpsysy(a,number,g)
173
c we determine the spline sp(u,v) according to the optimal values g(j).
174
80 call fpgrdi(ifsu,ifsv,ifbu,ifbv,0,u,mu,v,mv,z,mz,dz,
175
* iop0,iop1,tu,nu,tv,nv,p,c,nc,sq,fp,fpu,fpv,mm,mvnu,
176
* wrk(lsu),wrk(lsv),wrk(lri),wrk(lq),wrk(lau),wrk(lav1),
177
* wrk(lav2),wrk(lbu),wrk(lbv),wrk(laa),wrk(lbb),
178
* wrk(lcc),wrk(lcs),nru,nrv)
179
if(id0.eq.0) fp = fp+(z0-dz(1))**2
added
removed
Lib/interpolate/fitpack/fpopdi.f
