1
subroutine fppola(iopt1,iopt2,iopt3,m,u,v,z,w,rad,s,nuest,nvest,
2
* eta,tol,maxit,ib1,ib3,nc,ncc,intest,nrest,nu,tu,nv,tv,c,fp,sup,
3
* fpint,coord,f,ff,row,cs,cosi,a,q,bu,bv,spu,spv,h,index,nummer,
6
integer iopt1,iopt2,iopt3,m,nuest,nvest,maxit,ib1,ib3,nc,ncc,
7
* intest,nrest,nu,nv,lwrk,ier
8
real*8 s,eta,tol,fp,sup
10
integer index(nrest),nummer(m)
11
real*8 u(m),v(m),z(m),w(m),tu(nuest),tv(nvest),c(nc),fpint(intest)
13
* coord(intest),f(ncc),ff(nc),row(nvest),cs(nvest),cosi(5,nvest),
14
* a(ncc,ib1),q(ncc,ib3),bu(nuest,5),bv(nvest,5),spu(m,4),spv(m,4),
16
c ..user supplied function..
19
real*8 acc,arg,co,c1,c2,c3,c4,dmax,eps,fac,fac1,fac2,fpmax,fpms,
20
* f1,f2,f3,hui,huj,p,pi,pinv,piv,pi2,p1,p2,p3,r,ratio,si,sigma,
21
* sq,store,uu,u2,u3,wi,zi,rn,one,two,three,con1,con4,con9,half,ten
22
integer i,iband,iband3,iband4,ich1,ich3,ii,il,in,ipar,ipar1,irot,
23
* iter,i1,i2,i3,j,jl,jrot,j1,j2,k,l,la,lf,lh,ll,lu,lv,lwest,l1,l2,
24
* l3,l4,ncof,ncoff,nvv,nv4,nreg,nrint,nrr,nr1,nuu,nu4,num,num1,
25
* numin,nvmin,rank,iband1
28
c ..function references..
29
real*8 abs,atan,cos,fprati,sin,sqrt
31
c ..subroutine references..
32
c fporde,fpbspl,fpback,fpgivs,fprota,fprank,fpdisc,fprppo
45
ipar = iopt2*(iopt2+3)/2
48
if(iopt1.lt.0) go to 90
50
nvmin = 9+iopt2*(iopt2+1)
51
c calculation of acc, the absolute tolerance for the root of f(p)=s.
53
if(iopt1.eq.0) go to 10
55
if (nv.lt.nvmin) go to 70
58
c if iopt1 = 0 we begin by computing the weighted least-squares
59
c polymomial of the form
60
c s(u,v) = f(1)*(1-u**3)+f(2)*u**3+f(3)*(u**2-u**3)+f(4)*(u-u**3)
61
c where f(4) = 0 if iopt2> 0 , f(3) = 0 if iopt2 > 1 and
62
c f(2) = 0 if iopt3> 0.
63
c the corresponding weighted sum of squared residuals gives the upper
64
c bound sup for the smoothing factor s.
81
if(iopt3.ne.0) h(2) = 0.
82
if(iopt2.gt.1) h(3) = 0.
83
if(iopt2.gt.0) h(4) = 0.
86
if(piv.eq.0.) go to 40
87
call fpgivs(piv,a(j,1),co,si)
88
call fprota(co,si,zi,f(j))
94
call fprota(co,si,h(l),a(j,j2))
99
if(a(4,1).ne.0.) f(4) = f(4)/a(4,1)
100
if(a(3,1).ne.0.) f(3) = (f(3)-a(3,2)*f(4))/a(3,1)
101
if(a(2,1).ne.0.) f(2) = (f(2)-a(2,2)*f(3)-a(2,3)*f(4))/a(2,1)
103
* f(1) = (f(1)-a(1,2)*f(2)-a(1,3)*f(3)-a(1,4)*f(4))/a(1,1)
104
c find the b-spline representation of this least-squares polynomial
108
c3 = (f(3)+two*f(4))/three+c1
124
c test whether the least-squares polynomial is an acceptable solution
126
if(fpms.lt.acc) go to 960
127
c test whether we cannot further increase the number of knots.
128
70 if(nuest.lt.numin .or. nvest.lt.nvmin) go to 950
129
c find the initial set of interior knots of the spline in case iopt1=0.
140
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
141
c part 1 : computation of least-squares bicubic splines. c
142
c ****************************************************** c
143
c if iopt1<0 we compute the least-squares bicubic spline according c
144
c to the given set of knots. c
145
c if iopt1>=0 we compute least-squares bicubic splines with in- c
146
c creasing numbers of knots until the corresponding sum f(p=inf)<=s. c
147
c the initial set of knots then depends on the value of iopt1 c
148
c if iopt1=0 we start with one interior knot in the u-direction c
149
c (0.5) and 1+iopt2*(iopt2+1) in the v-direction. c
150
c if iopt1>0 we start with the set of knots found at the last c
151
c call of the routine. c
152
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
153
c main loop for the different sets of knots. m is a save upper bound
154
c for the number of trials.
156
c find the position of the additional knots which are needed for the
157
c b-spline representation of s(u,v).
178
c find nrint, the total number of knot intervals and nreg, the number
179
c of panels in which the approximation domain is subdivided by the
180
c intersection of knots.
187
c arrange the data points according to the panel they belong to.
188
call fporde(u,v,m,3,3,tu,nu,tv,nv,nummer,index,nreg)
189
if(iopt2.eq.0) go to 195
190
c find the b-spline coefficients cosi of the cubic spline
191
c approximations for cr(v)=rad(v)*cos(v) and sr(v) = rad(v)*sin(v)
192
c if iopt2=1, and additionally also for cr(v)**2,sr(v)**2 and
193
c 2*cr(v)*sr(v) if iopt2=2
201
c the coefficients cosi are obtained from interpolation conditions
202
c at the knots tv(i),i=4,5,...nv-4.
206
call fpbspl(tv,nv,3,arg,l2,hv)
213
row(ll) = row(ll)+hv(j)
221
if(iopt2.eq.1) go to 155
227
if(piv.eq.0.) go to 170
228
call fpgivs(piv,a(j,1),co,si)
230
call fprota(co,si,cs(l),cosi(l,j))
232
if(j.eq.nvv) go to 175
237
call fprota(co,si,row(l),a(j,j2))
245
call fpback(a,cs,nvv,nvv,cs,ncc)
250
c find ncof, the dimension of the spline and ncoff, the number
251
c of coefficients in the standard b-spline representation.
255
ncof = ipar1+nvv*(nu4-1-iopt2-iopt3)
256
c find the bandwidth of the observation matrix a.
258
if(nuu-iopt2-iopt3.le.1) iband = ncof
260
c initialize the observation matrix a.
266
c initialize the sum of squared residuals.
268
ratio = one+tu(6)/tu(5)
269
c fetch the data points in the new order. main loop for the
272
c fix certain constants for the current panel; jrot records the column
273
c number of the first non-zero element in a row of the observation
274
c matrix according to a data point of the panel.
281
if(lu.gt.iopt2) jrot = ipar1+(lu-iopt2-1)*nvv
283
c test whether there are still data points in the current panel.
285
210 if(in.eq.0) go to 380
286
c fetch a new data point.
289
c evaluate for the u-direction, the 4 non-zero b-splines at u(in)
290
call fpbspl(tu,nu,3,u(in),l1,hu)
291
c evaluate for the v-direction, the 4 non-zero b-splines at v(in)
292
call fpbspl(tv,nv,3,v(in),l2,hv)
293
c store the value of these b-splines in spu and spv resp.
298
c initialize the new row of observation matrix.
302
c calculate the non-zero elements of the new row by making the cross
303
c products of the non-zero b-splines in u- and v-direction and
304
c by taking into account the conditions of the splines.
308
c take into account the periodicity condition of the bicubic splines.
312
row(ll) = row(ll)+hv(i)
315
c take into account the other conditions of the splines.
316
if(iopt2.eq.0 .or. lu.gt.iopt2+1) go to 280
320
cs(l) = cs(l)+row(i)*cosi(l,i)
322
c fill in the non-zero elements of the new row.
327
if(jlu.gt.iopt2+2) go to 320
328
go to (290,290,300,310),jlu
338
h(2) = h(2)+huj*ratio*cs(1)
339
h(3) = h(3)+huj*ratio*cs(2)
345
320 if(jlu.gt.nu4 .and. iopt3.ne.0) go to 330
354
c rotate the row into triangle by givens transformations.
359
if(piv.eq.0.) go to 350
360
c calculate the parameters of the givens transformation.
361
call fpgivs(piv,a(irot,1),co,si)
362
c apply that transformation to the right hand side.
363
call fprota(co,si,zi,f(irot))
364
if(i.eq.iband) go to 360
365
c apply that transformation to the left hand side.
370
call fprota(co,si,h(j),a(irot,i2))
373
c add the contribution of the row to the sum of squares of residual
376
c find the number of the next data point in the panel.
380
c find dmax, the maximum value for the diagonal elements in the reduced
384
if(a(i,1).le.dmax) go to 390
387
c check whether the observation matrix is rank deficient.
390
if(a(i,1).le.sigma) go to 410
392
c backward substitution in case of full rank.
393
call fpback(a,f,ncof,iband,c,ncc)
399
c in case of rank deficiency, find the minimum norm solution.
400
410 lwest = ncof*iband+ncof+iband
401
if(lwrk.lt.lwest) go to 925
410
call fprank(q,ff,ncof,iband,ncc,sigma,c,sq,rank,wrk(la),
415
c add to the sum of squared residuals, the contribution of reducing
418
c find the coefficients in the standard b-spline representation of
420
430 call fprppo(nu,nv,iopt2,iopt3,cosi,ratio,c,ff,ncoff)
421
c test whether the least-squares spline is an acceptable solution.
423
if (fp.le.0) go to 970
427
if(abs(fpms).le.acc) then
428
if (fp.le.0) go to 970
431
c if f(p=inf) < s, accept the choice of knots.
432
if(fpms.lt.0.) go to 580
433
c test whether we cannot further increase the number of knots
434
if(m.lt.ncof) go to 935
435
c search where to add a new knot.
436
c find for each interval the sum of squared residuals fpint for the
437
c data points having the coordinate belonging to that knot interval.
438
c calculate also coord which is the same sum, weighted by the position
439
c of the data points considered.
452
460 if(in.eq.0) go to 490
460
store = store+hui*spv(in,j)*c(j1)
464
store = (w(in)*(z(in)-store))**2
465
fpint(l1) = fpint(l1)+store
466
coord(l1) = coord(l1)+store*u(in)
467
fpint(l2) = fpint(l2)+store
468
coord(l2) = coord(l2)+store*v(in)
472
c bring together the information concerning knot panels which are
473
c symmetric with respect to the origin.
477
fpint(l1) = fpint(l1)+fpint(l2)
478
coord(l1) = coord(l1)+coord(l2)-pi*fpint(l2)
480
c find the interval for which fpint is maximal on the condition that
481
c there still can be added a knot.
484
if(nuest.lt.nu+1) l1=nuu+1
485
if(nvest.lt.nv+2) l2=nuu
486
c test whether we cannot further increase the number of knots.
487
if(l1.gt.l2) go to 950
491
if(fpmax.ge.fpint(i)) go to 510
496
c calculate the position of the new knot.
497
arg = coord(l)/fpint(l)
498
c test in what direction the new knot is going to be added.
499
if(l.gt.nuu) go to 530
500
c addition in the u-direction
505
if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 500
514
c addition in the v-direction
519
if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 500
533
c restart the computations with the new set of knots.
535
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
536
c part 2: determination of the smoothing bicubic spline. c
537
c ****************************************************** c
538
c we have determined the number of knots and their position. we now c
539
c compute the coefficients of the smoothing spline sp(u,v). c
540
c the observation matrix a is extended by the rows of a matrix, expres-c
541
c sing that sp(u,v) must be a constant function in the variable c
542
c v and a cubic polynomial in the variable u. the corresponding c
543
c weights of these additional rows are set to 1/(p). iteratively c
544
c we than have to determine the value of p such that f(p) = sum((w(i)* c
545
c (z(i)-sp(u(i),v(i))))**2) be = s. c
546
c we already know that the least-squares polynomial corresponds to p=0,c
547
c and that the least-squares bicubic spline corresponds to p=infin. c
548
c the iteration process makes use of rational interpolation. since f(p)c
549
c is a convex and strictly decreasing function of p, it can be approx- c
550
c imated by a rational function of the form r(p) = (u*p+v)/(p+w). c
551
c three values of p (p1,p2,p3) with corresponding values of f(p) (f1= c
552
c f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the new value c
553
c of p such that r(p)=s. convergence is guaranteed by taking f1>0,f3<0.c
554
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
555
c evaluate the discontinuity jumps of the 3-th order derivative of
556
c the b-splines at the knots tu(l),l=5,...,nu-4.
557
580 call fpdisc(tu,nu,5,bu,nuest)
558
c evaluate the discontinuity jumps of the 3-th order derivative of
559
c the b-splines at the knots tv(l),l=5,...,nv-4.
560
call fpdisc(tv,nv,5,bv,nvest)
561
c initial value for p.
572
c find the bandwidth of the extended observation matrix.
574
if(iband4.gt.ncof) iband4 = ncof
579
c iteration process to find the root of f(p)=s.
582
c store the triangularized observation matrix into q.
591
c extend the observation matrix with the rows of a matrix, expressing
592
c that for u=constant sp(u,v) must be a constant function.
601
row(ll) = row(ll)+bv(ii,l)
605
c initialize the new row.
609
c fill in the non-zero elements of the row. jrot records the column
610
c number of the first non-zero element in the row.
611
if(j.gt.iopt2) go to 665
616
cs(k) = cs(k)+cosi(k,l)*row(l)
625
cs(k) = cs(k)+cosi(k,l)*row(l)
637
jrot = ipar1+1+(j-iopt2-1)*nvv
642
c rotate the new row into triangle by givens transformations.
643
do 710 irot=jrot,ncof
645
i2 = min0(iband1,ncof-irot)
647
if (i2.le.0) go to 720
650
c calculate the parameters of the givens transformation.
651
call fpgivs(piv,q(irot,1),co,si)
652
c apply that givens transformation to the right hand side.
653
call fprota(co,si,zi,ff(irot))
654
if(i2.eq.0) go to 720
655
c apply that givens transformation to the left hand side.
658
call fprota(co,si,h(l1),q(irot,l1))
666
c extend the observation matrix with the rows of a matrix expressing
667
c that for v=constant. sp(u,v) must be a cubic polynomial.
671
c initialize the new row
675
c fill in the non-zero elements of the row. jrot records the column
676
c number of the first non-zero element in the row.
680
if(il.eq.nu4 .and. iopt3.ne.0) go to 760
681
if(il.gt.iopt2+1) go to 750
682
go to (735,740,745),il
686
740 h(1) = h(1)+bu(ii,l)
687
h(2) = bu(ii,l)*cosi(1,j)
688
h(3) = bu(ii,l)*cosi(2,j)
691
745 h(1) = h(1)+bu(ii,l)
692
h(2) = bu(ii,l)*cosi(1,j)*ratio
693
h(3) = bu(ii,l)*cosi(2,j)*ratio
694
h(4) = bu(ii,l)*cosi(3,j)
695
h(5) = bu(ii,l)*cosi(4,j)
696
h(6) = bu(ii,l)*cosi(5,j)
707
if(ii.gt.iopt2+1) jrot = ipar1+(ii-iopt2-2)*nvv+j
708
c rotate the new row into triangle by givens transformations.
709
do 800 irot=jrot,ncof
711
i2 = min0(iband3,ncof-irot)
713
if (i2.le.0) go to 810
716
c calculate the parameters of the givens transformation.
717
call fpgivs(piv,q(irot,1),co,si)
718
c apply that givens transformation to the right hand side.
719
call fprota(co,si,zi,ff(irot))
720
if(i2.eq.0) go to 810
721
c apply that givens transformation to the left hand side.
724
call fprota(co,si,h(l1),q(irot,l1))
732
c find dmax, the maximum value for the diagonal elements in the
736
if(q(i,1).le.dmax) go to 820
739
c check whether the matrix is rank deficient.
742
if(q(i,1).le.sigma) go to 840
744
c backward substitution in case of full rank.
745
call fpback(q,ff,ncof,iband4,c,ncc)
748
c in case of rank deficiency, find the minimum norm solution.
749
840 lwest = ncof*iband4+ncof+iband4
750
if(lwrk.lt.lwest) go to 925
754
call fprank(q,ff,ncof,iband4,ncc,sigma,c,sq,rank,wrk(la),
759
c find the coefficients in the standard b-spline representation of
761
call fprppo(nu,nv,iopt2,iopt3,cosi,ratio,c,ff,ncoff)
770
860 if(in.eq.0) go to 890
778
store = store+hui*spv(in,j)*c(j1)
782
fp = fp+(w(in)*(z(in)-store))**2
786
c test whether the approximation sp(u,v) is an acceptable solution
788
if(abs(fpms).le.acc) go to 980
789
c test whether the maximum allowable number of iterations has been
791
if(iter.eq.maxit) go to 940
792
c carry out one more step of the iteration process.
795
if(ich3.ne.0) go to 900
796
if((f2-f3).gt.acc) go to 895
797
c our initial choice of p is too large.
801
if(p.le.p1) p = p1*con9 + p2*con1
803
895 if(f2.lt.0.) ich3 = 1
804
900 if(ich1.ne.0) go to 910
805
if((f1-f2).gt.acc) go to 905
806
c our initial choice of p is too small
810
if(p3.lt.0.) go to 920
811
if(p.ge.p3) p = p2*con1 +p3*con9
813
905 if(f2.gt.0.) ich1 = 1
814
c test whether the iteration process proceeds as theoretically
816
910 if(f2.ge.f1 .or. f2.le.f3) go to 945
817
c find the new value of p.
818
p = fprati(p1,f1,p2,f2,p3,f3)
820
c error codes and messages.
837
980 if(ncof.ne.rank) ier = -rank
added
removed
Lib/sandbox/spline/fitpack/fppola.f
