1
subroutine fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,
2
* ty,ny,p,c,nc,fp,fpx,fpy,mm,mynx,kx1,kx2,ky1,ky2,spx,spy,right,q,
7
integer ifsx,ifsy,ifbx,ifby,mx,my,mz,kx,ky,nx,ny,nc,mm,mynx,
10
real*8 x(mx),y(my),z(mz),tx(nx),ty(ny),c(nc),spx(mx,kx1),spy(my,ky
12
* ,right(mm),q(mynx),ax(nx,kx2),bx(nx,kx2),ay(ny,ky2),by(ny,ky2),
14
integer nrx(mx),nry(my)
16
real*8 arg,cos,fac,pinv,piv,sin,term,one,half
17
integer i,ibandx,ibandy,ic,iq,irot,it,iz,i1,i2,i3,j,k,k1,k2,l,
18
* l1,l2,ncof,nk1x,nk1y,nrold,nroldx,nroldy,number,numx,numx1,
22
c ..subroutine references..
23
c fpback,fpbspl,fpgivs,fpdisc,fprota
25
c the b-spline coefficients of the smoothing spline are calculated as
26
c the least-squares solution of the over-determined linear system of
27
c equations (ay) c (ax)' = q where
30
c (ax) = | ---------- | (ay) = | ---------- |
31
c | (1/p) (bx) | | (1/p) (by) |
37
c with c : the (ny-ky-1) x (nx-kx-1) matrix which contains the
38
c b-spline coefficients.
39
c z : the my x mx matrix which contains the function values.
40
c spx,spy: the mx x (nx-kx-1) and my x (ny-ky-1) observation
41
c matrices according to the least-squares problems in
42
c the x- and y-direction.
43
c bx,by : the (nx-2*kx-1) x (nx-kx-1) and (ny-2*ky-1) x (ny-ky-1)
44
c matrices which contain the discontinuity jumps of the
45
c derivatives of the b-splines in the x- and y-direction.
50
if(p.gt.0.) pinv = one/p
51
c it depends on the value of the flags ifsx,ifsy,ifbx and ifby and on
52
c the value of p whether the matrices (spx),(spy),(bx) and (by) still
54
if(ifsx.ne.0) go to 50
55
c calculate the non-zero elements of the matrix (spx) which is the
56
c observation matrix according to the least-squares spline approximat-
57
c ion problem in the x-direction.
63
10 if(arg.lt.tx(l1) .or. l.eq.nk1x) go to 20
68
20 call fpbspl(tx,nx,kx,arg,l,h)
75
50 if(ifsy.ne.0) go to 100
76
c calculate the non-zero elements of the matrix (spy) which is the
77
c observation matrix according to the least-squares spline approximat-
78
c ion problem in the y-direction.
84
60 if(arg.lt.ty(l1) .or. l.eq.nk1y) go to 70
89
70 call fpbspl(ty,ny,ky,arg,l,h)
96
100 if(p.le.0.) go to 120
97
c calculate the non-zero elements of the matrix (bx).
98
if(ifbx.ne.0 .or. nx.eq.2*kx1) go to 110
99
call fpdisc(tx,nx,kx2,bx,nx)
101
c calculate the non-zero elements of the matrix (by).
102
110 if(ifby.ne.0 .or. ny.eq.2*ky1) go to 120
103
call fpdisc(ty,ny,ky2,by,ny)
105
c reduce the matrix (ax) to upper triangular form (rx) using givens
106
c rotations. apply the same transformations to the rows of matrix q
107
c to obtain the my x (nx-kx-1) matrix g.
108
c store matrix (rx) into (ax) and g into q.
120
c ibandx denotes the bandwidth of the matrices (ax) and (rx).
124
150 if(nrold.eq.number) go to 180
125
if(p.le.0.) go to 260
127
c fetch a new row of matrix (bx).
132
c find the appropriate column of q.
138
c fetch a new row of matrix (spx).
143
c find the appropriate column of q.
149
c rotate the new row of matrix (ax) into triangle.
150
210 do 240 i=1,ibandx
153
if(piv.eq.0.) go to 240
154
c calculate the parameters of the givens transformation.
155
call fpgivs(piv,ax(irot,1),cos,sin)
156
c apply that transformation to the rows of matrix q.
160
call fprota(cos,sin,right(j),q(iq))
162
c apply that transformation to the columns of (ax).
163
if(i.eq.ibandx) go to 250
168
call fprota(cos,sin,h(j),ax(irot,i2))
171
250 if(nrold.eq.number) go to 270
175
c reduce the matrix (ay) to upper triangular form (ry) using givens
176
c rotations. apply the same transformations to the columns of matrix g
177
c to obtain the (ny-ky-1) x (nx-kx-1) matrix h.
178
c store matrix (ry) into (ay) and h into c.
189
c ibandy denotes the bandwidth of the matrices (ay) and (ry).
193
300 if(nrold.eq.number) go to 330
194
if(p.le.0.) go to 410
196
c fetch a new row of matrix (by).
201
c find the appropiate row of g.
207
c fetch a new row of matrix (spy)
212
c find the appropiate row of g.
219
c rotate the new row of matrix (ay) into triangle.
220
360 do 390 i=1,ibandy
223
if(piv.eq.0.) go to 390
224
c calculate the parameters of the givens transformation.
225
call fpgivs(piv,ay(irot,1),cos,sin)
226
c apply that transformation to the colums of matrix g.
229
call fprota(cos,sin,right(j),c(ic))
232
c apply that transformation to the columns of matrix (ay).
233
if(i.eq.ibandy) go to 400
238
call fprota(cos,sin,h(j),ay(irot,i2))
241
400 if(nrold.eq.number) go to 420
245
c backward substitution to obtain the b-spline coefficients as the
246
c solution of the linear system (ry) c (rx)' = h.
247
c first step: solve the system (ry) (c1) = h.
250
call fpback(ay,c(k),nk1y,ibandy,c(k),ny)
253
c second step: solve the system c (rx)' = (c1).
262
call fpback(ax,right,nk1x,ibandx,right,nx)
269
c calculate the quantities
270
c res(i,j) = (z(i,j) - s(x(i),y(j)))**2 , i=1,2,..,mx;j=1,2,..,my
271
c fp = sumi=1,mx(sumj=1,my(res(i,j)))
272
c fpx(r) = sum''i(sumj=1,my(res(i,j))) , r=1,2,...,nx-2*kx-1
273
c tx(r+kx) <= x(i) <= tx(r+kx+1)
274
c fpy(r) = sumi=1,mx(sum''j(res(i,j))) , r=1,2,...,ny-2*ky-1
275
c ty(r+ky) <= y(j) <= ty(r+ky+1)
286
c main loop for the different grid points.
295
c evaluate s(x,y) at the current grid point by making the sum of the
296
c cross products of the non-zero b-splines at (x,y), multiplied with
297
c the appropiate b-spline coefficients.
305
term = term+fac*spy(i2,l2)*c(k2)
309
c calculate the squared residual at the current grid point.
310
term = (z(iz)-term)**2
311
c adjust the different parameters.
313
fpx(numx1) = fpx(numx1)+term
314
fpy(numy1) = fpy(numy1)+term
316
if(numy.eq.nroldy) go to 530
317
fpy(numy1) = fpy(numy1)-fac
318
fpy(numy) = fpy(numy)+fac
320
if(numx.eq.nroldx) go to 540
321
fpx(numx1) = fpx(numx1)-fac
322
fpx(numx) = fpx(numx)+fac
added
removed
Lib/interpolate/fitpack/fpgrre.f
