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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
|
subroutine fpclos(iopt,idim,m,u,mx,x,w,k,s,nest,tol,maxit,k1,k2,
* n,t,nc,c,fp,fpint,z,a1,a2,b,g1,g2,q,nrdata,ier)
c ..
c ..scalar arguments..
real*8 s,tol,fp
integer iopt,idim,m,mx,k,nest,maxit,k1,k2,n,nc,ier
c ..array arguments..
real*8 u(m),x(mx),w(m),t(nest),c(nc),fpint(nest),z(nc),a1(nest,k1)
*,
* a2(nest,k),b(nest,k2),g1(nest,k2),g2(nest,k1),q(m,k1)
integer nrdata(nest)
c ..local scalars..
real*8 acc,cos,d1,fac,fpart,fpms,fpold,fp0,f1,f2,f3,p,per,pinv,piv
*,
* p1,p2,p3,sin,store,term,ui,wi,rn,one,con1,con4,con9,half
integer i,ich1,ich3,ij,ik,it,iter,i1,i2,i3,j,jj,jk,jper,j1,j2,kk,
* kk1,k3,l,l0,l1,l5,mm,m1,new,nk1,nk2,nmax,nmin,nplus,npl1,
* nrint,n10,n11,n7,n8
c ..local arrays..
real*8 h(6),h1(7),h2(6),xi(10)
c ..function references..
real*8 abs,fprati
integer max0,min0
c ..subroutine references..
c fpbacp,fpbspl,fpgivs,fpdisc,fpknot,fprota
c ..
c set constants
one = 0.1e+01
con1 = 0.1e0
con9 = 0.9e0
con4 = 0.4e-01
half = 0.5e0
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c part 1: determination of the number of knots and their position c
c ************************************************************** c
c given a set of knots we compute the least-squares closed curve c
c sinf(u). if the sum f(p=inf) <= s we accept the choice of knots. c
c if iopt=-1 sinf(u) is the requested curve c
c if iopt=0 or iopt=1 we check whether we can accept the knots: c
c if fp <=s we will continue with the current set of knots. c
c if fp > s we will increase the number of knots and compute the c
c corresponding least-squares curve until finally fp<=s. c
c the initial choice of knots depends on the value of s and iopt. c
c if s=0 we have spline interpolation; in that case the number of c
c knots equals nmax = m+2*k. c
c if s > 0 and c
c iopt=0 we first compute the least-squares polynomial curve of c
c degree k; n = nmin = 2*k+2. since s(u) must be periodic we c
c find that s(u) reduces to a fixed point. c
c iopt=1 we start with the set of knots found at the last c
c call of the routine, except for the case that s > fp0; then c
c we compute directly the least-squares polynomial curve. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
m1 = m-1
kk = k
kk1 = k1
k3 = 3*k+1
nmin = 2*k1
c determine the length of the period of the splines.
per = u(m)-u(1)
if(iopt.lt.0) go to 50
c calculation of acc, the absolute tolerance for the root of f(p)=s.
acc = tol*s
c determine nmax, the number of knots for periodic spline interpolation
nmax = m+2*k
if(s.gt.0. .or. nmax.eq.nmin) go to 30
c if s=0, s(u) is an interpolating curve.
n = nmax
c test whether the required storage space exceeds the available one.
if(n.gt.nest) go to 620
c find the position of the interior knots in case of interpolation.
5 if((k/2)*2 .eq.k) go to 20
do 10 i=2,m1
j = i+k
t(j) = u(i)
10 continue
if(s.gt.0.) go to 50
kk = k-1
kk1 = k
if(kk.gt.0) go to 50
t(1) = t(m)-per
t(2) = u(1)
t(m+1) = u(m)
t(m+2) = t(3)+per
jj = 0
do 15 i=1,m1
j = i
do 12 j1=1,idim
jj = jj+1
c(j) = x(jj)
j = j+n
12 continue
15 continue
jj = 1
j = m
do 17 j1=1,idim
c(j) = c(jj)
j = j+n
jj = jj+n
17 continue
fp = 0.
fpint(n) = fp0
fpint(n-1) = 0.
nrdata(n) = 0
go to 630
20 do 25 i=2,m1
j = i+k
t(j) = (u(i)+u(i-1))*half
25 continue
go to 50
c if s > 0 our initial choice depends on the value of iopt.
c if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
c polynomial curve. (i.e. a constant point).
c if iopt=1 and fp0>s we start computing the least-squares closed
c curve according the set of knots found at the last call of the
c routine.
30 if(iopt.eq.0) go to 35
if(n.eq.nmin) go to 35
fp0 = fpint(n)
fpold = fpint(n-1)
nplus = nrdata(n)
if(fp0.gt.s) go to 50
c the case that s(u) is a fixed point is treated separetely.
c fp0 denotes the corresponding sum of squared residuals.
35 fp0 = 0.
d1 = 0.
do 37 j=1,idim
z(j) = 0.
37 continue
jj = 0
do 45 it=1,m1
wi = w(it)
call fpgivs(wi,d1,cos,sin)
do 40 j=1,idim
jj = jj+1
fac = wi*x(jj)
call fprota(cos,sin,fac,z(j))
fp0 = fp0+fac**2
40 continue
45 continue
do 47 j=1,idim
z(j) = z(j)/d1
47 continue
c test whether that fixed point is a solution of our problem.
fpms = fp0-s
if(fpms.lt.acc .or. nmax.eq.nmin) go to 640
fpold = fp0
c test whether the required storage space exceeds the available one.
if(n.ge.nest) go to 620
c start computing the least-squares closed curve with one
c interior knot.
nplus = 1
n = nmin+1
mm = (m+1)/2
t(k2) = u(mm)
nrdata(1) = mm-2
nrdata(2) = m1-mm
c main loop for the different sets of knots. m is a save upper
c bound for the number of trials.
50 do 340 iter=1,m
c find nrint, the number of knot intervals.
nrint = n-nmin+1
c find the position of the additional knots which are needed for
c the b-spline representation of s(u). if we take
c t(k+1) = u(1), t(n-k) = u(m)
c t(k+1-j) = t(n-k-j) - per, j=1,2,...k
c t(n-k+j) = t(k+1+j) + per, j=1,2,...k
c then s(u) will be a smooth closed curve if the b-spline
c coefficients satisfy the following conditions
c c((i-1)*n+n7+j) = c((i-1)*n+j), j=1,...k,i=1,2,...,idim (**)
c with n7=n-2*k-1.
t(k1) = u(1)
nk1 = n-k1
nk2 = nk1+1
t(nk2) = u(m)
do 60 j=1,k
i1 = nk2+j
i2 = nk2-j
j1 = k1+j
j2 = k1-j
t(i1) = t(j1)+per
t(j2) = t(i2)-per
60 continue
c compute the b-spline coefficients of the least-squares closed curve
c sinf(u). the observation matrix a is built up row by row while
c taking into account condition (**) and is reduced to triangular
c form by givens transformations .
c at the same time fp=f(p=inf) is computed.
c the n7 x n7 triangularised upper matrix a has the form
c ! a1 ' !
c a = ! ' a2 !
c ! 0 ' !
c with a2 a n7 x k matrix and a1 a n10 x n10 upper triangular
c matrix of bandwith k+1 ( n10 = n7-k).
c initialization.
do 65 i=1,nc
z(i) = 0.
65 continue
do 70 i=1,nk1
do 70 j=1,kk1
a1(i,j) = 0.
70 continue
n7 = nk1-k
n10 = n7-kk
jper = 0
fp = 0.
l = k1
jj = 0
do 290 it=1,m1
c fetch the current data point u(it),x(it)
ui = u(it)
wi = w(it)
do 75 j=1,idim
jj = jj+1
xi(j) = x(jj)*wi
75 continue
c search for knot interval t(l) <= ui < t(l+1).
80 if(ui.lt.t(l+1)) go to 85
l = l+1
go to 80
c evaluate the (k+1) non-zero b-splines at ui and store them in q.
85 call fpbspl(t,n,k,ui,l,h)
do 90 i=1,k1
q(it,i) = h(i)
h(i) = h(i)*wi
90 continue
l5 = l-k1
c test whether the b-splines nj,k+1(u),j=1+n7,...nk1 are all zero at ui
if(l5.lt.n10) go to 285
if(jper.ne.0) go to 160
c initialize the matrix a2.
do 95 i=1,n7
do 95 j=1,kk
a2(i,j) = 0.
95 continue
jk = n10+1
do 110 i=1,kk
ik = jk
do 100 j=1,kk1
if(ik.le.0) go to 105
a2(ik,i) = a1(ik,j)
ik = ik-1
100 continue
105 jk = jk+1
110 continue
jper = 1
c if one of the b-splines nj,k+1(u),j=n7+1,...nk1 is not zero at ui
c we take account of condition (**) for setting up the new row
c of the observation matrix a. this row is stored in the arrays h1
c (the part with respect to a1) and h2 (the part with
c respect to a2).
160 do 170 i=1,kk
h1(i) = 0.
h2(i) = 0.
170 continue
h1(kk1) = 0.
j = l5-n10
do 210 i=1,kk1
j = j+1
l0 = j
180 l1 = l0-kk
if(l1.le.0) go to 200
if(l1.le.n10) go to 190
l0 = l1-n10
go to 180
190 h1(l1) = h(i)
go to 210
200 h2(l0) = h2(l0)+h(i)
210 continue
c rotate the new row of the observation matrix into triangle
c by givens transformations.
if(n10.le.0) go to 250
c rotation with the rows 1,2,...n10 of matrix a.
do 240 j=1,n10
piv = h1(1)
if(piv.ne.0.) go to 214
do 212 i=1,kk
h1(i) = h1(i+1)
212 continue
h1(kk1) = 0.
go to 240
c calculate the parameters of the givens transformation.
214 call fpgivs(piv,a1(j,1),cos,sin)
c transformation to the right hand side.
j1 = j
do 217 j2=1,idim
call fprota(cos,sin,xi(j2),z(j1))
j1 = j1+n
217 continue
c transformations to the left hand side with respect to a2.
do 220 i=1,kk
call fprota(cos,sin,h2(i),a2(j,i))
220 continue
if(j.eq.n10) go to 250
i2 = min0(n10-j,kk)
c transformations to the left hand side with respect to a1.
do 230 i=1,i2
i1 = i+1
call fprota(cos,sin,h1(i1),a1(j,i1))
h1(i) = h1(i1)
230 continue
h1(i1) = 0.
240 continue
c rotation with the rows n10+1,...n7 of matrix a.
250 do 270 j=1,kk
ij = n10+j
if(ij.le.0) go to 270
piv = h2(j)
if(piv.eq.0.) go to 270
c calculate the parameters of the givens transformation.
call fpgivs(piv,a2(ij,j),cos,sin)
c transformations to right hand side.
j1 = ij
do 255 j2=1,idim
call fprota(cos,sin,xi(j2),z(j1))
j1 = j1+n
255 continue
if(j.eq.kk) go to 280
j1 = j+1
c transformations to left hand side.
do 260 i=j1,kk
call fprota(cos,sin,h2(i),a2(ij,i))
260 continue
270 continue
c add contribution of this row to the sum of squares of residual
c right hand sides.
280 do 282 j2=1,idim
fp = fp+xi(j2)**2
282 continue
go to 290
c rotation of the new row of the observation matrix into
c triangle in case the b-splines nj,k+1(u),j=n7+1,...n-k-1 are all zero
c at ui.
285 j = l5
do 140 i=1,kk1
j = j+1
piv = h(i)
if(piv.eq.0.) go to 140
c calculate the parameters of the givens transformation.
call fpgivs(piv,a1(j,1),cos,sin)
c transformations to right hand side.
j1 = j
do 125 j2=1,idim
call fprota(cos,sin,xi(j2),z(j1))
j1 = j1+n
125 continue
if(i.eq.kk1) go to 150
i2 = 1
i3 = i+1
c transformations to left hand side.
do 130 i1=i3,kk1
i2 = i2+1
call fprota(cos,sin,h(i1),a1(j,i2))
130 continue
140 continue
c add contribution of this row to the sum of squares of residual
c right hand sides.
150 do 155 j2=1,idim
fp = fp+xi(j2)**2
155 continue
290 continue
fpint(n) = fp0
fpint(n-1) = fpold
nrdata(n) = nplus
c backward substitution to obtain the b-spline coefficients .
j1 = 1
do 292 j2=1,idim
call fpbacp(a1,a2,z(j1),n7,kk,c(j1),kk1,nest)
j1 = j1+n
292 continue
c calculate from condition (**) the remaining coefficients.
do 297 i=1,k
j1 = i
do 295 j=1,idim
j2 = j1+n7
c(j2) = c(j1)
j1 = j1+n
295 continue
297 continue
if(iopt.lt.0) go to 660
c test whether the approximation sinf(u) is an acceptable solution.
fpms = fp-s
if(abs(fpms).lt.acc) go to 660
c if f(p=inf) < s accept the choice of knots.
if(fpms.lt.0.) go to 350
c if n=nmax, sinf(u) is an interpolating curve.
if(n.eq.nmax) go to 630
c increase the number of knots.
c if n=nest we cannot increase the number of knots because of the
c storage capacity limitation.
if(n.eq.nest) go to 620
c determine the number of knots nplus we are going to add.
npl1 = nplus*2
rn = nplus
if(fpold-fp.gt.acc) npl1 = rn*fpms/(fpold-fp)
nplus = min0(nplus*2,max0(npl1,nplus/2,1))
fpold = fp
c compute the sum of squared residuals for each knot interval
c t(j+k) <= ui <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
fpart = 0.
i = 1
l = k1
jj = 0
do 320 it=1,m1
if(u(it).lt.t(l)) go to 300
new = 1
l = l+1
300 term = 0.
l0 = l-k2
do 310 j2=1,idim
fac = 0.
j1 = l0
do 305 j=1,k1
j1 = j1+1
fac = fac+c(j1)*q(it,j)
305 continue
jj = jj+1
term = term+(w(it)*(fac-x(jj)))**2
l0 = l0+n
310 continue
fpart = fpart+term
if(new.eq.0) go to 320
if(l.gt.k2) go to 315
fpint(nrint) = term
new = 0
go to 320
315 store = term*half
fpint(i) = fpart-store
i = i+1
fpart = store
new = 0
320 continue
fpint(nrint) = fpint(nrint)+fpart
do 330 l=1,nplus
c add a new knot
call fpknot(u,m,t,n,fpint,nrdata,nrint,nest,1)
c if n=nmax we locate the knots as for interpolation
if(n.eq.nmax) go to 5
c test whether we cannot further increase the number of knots.
if(n.eq.nest) go to 340
330 continue
c restart the computations with the new set of knots.
340 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c part 2: determination of the smoothing closed curve sp(u). c
c ********************************************************** c
c we have determined the number of knots and their position. c
c we now compute the b-spline coefficients of the smoothing curve c
c sp(u). the observation matrix a is extended by the rows of matrix c
c b expressing that the kth derivative discontinuities of sp(u) at c
c the interior knots t(k+2),...t(n-k-1) must be zero. the corres- c
c ponding weights of these additional rows are set to 1/p. c
c iteratively we then have to determine the value of p such that f(p),c
c the sum of squared residuals be = s. we already know that the least-c
c squares polynomial curve corresponds to p=0, and that the least- c
c squares periodic spline curve corresponds to p=infinity. the c
c iteration process which is proposed here, makes use of rational c
c interpolation. since f(p) is a convex and strictly decreasing c
c function of p, it can be approximated by a rational function c
c r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond- c
c ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used c
c to calculate the new value of p such that r(p)=s. convergence is c
c guaranteed by taking f1>0 and f3<0. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c evaluate the discontinuity jump of the kth derivative of the
c b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
350 call fpdisc(t,n,k2,b,nest)
c initial value for p.
p1 = 0.
f1 = fp0-s
p3 = -one
f3 = fpms
n11 = n10-1
n8 = n7-1
p = 0.
l = n7
do 352 i=1,k
j = k+1-i
p = p+a2(l,j)
l = l-1
if(l.eq.0) go to 356
352 continue
do 354 i=1,n10
p = p+a1(i,1)
354 continue
356 rn = n7
p = rn/p
ich1 = 0
ich3 = 0
c iteration process to find the root of f(p) = s.
do 595 iter=1,maxit
c form the matrix g as the matrix a extended by the rows of matrix b.
c the rows of matrix b with weight 1/p are rotated into
c the triangularised observation matrix a.
c after triangularisation our n7 x n7 matrix g takes the form
c ! g1 ' !
c g = ! ' g2 !
c ! 0 ' !
c with g2 a n7 x (k+1) matrix and g1 a n11 x n11 upper triangular
c matrix of bandwidth k+2. ( n11 = n7-k-1)
pinv = one/p
c store matrix a into g
do 358 i=1,nc
c(i) = z(i)
358 continue
do 360 i=1,n7
g1(i,k1) = a1(i,k1)
g1(i,k2) = 0.
g2(i,1) = 0.
do 360 j=1,k
g1(i,j) = a1(i,j)
g2(i,j+1) = a2(i,j)
360 continue
l = n10
do 370 j=1,k1
if(l.le.0) go to 375
g2(l,1) = a1(l,j)
l = l-1
370 continue
375 do 540 it=1,n8
c fetch a new row of matrix b and store it in the arrays h1 (the part
c with respect to g1) and h2 (the part with respect to g2).
do 380 j=1,idim
xi(j) = 0.
380 continue
do 385 i=1,k1
h1(i) = 0.
h2(i) = 0.
385 continue
h1(k2) = 0.
if(it.gt.n11) go to 420
l = it
l0 = it
do 390 j=1,k2
if(l0.eq.n10) go to 400
h1(j) = b(it,j)*pinv
l0 = l0+1
390 continue
go to 470
400 l0 = 1
do 410 l1=j,k2
h2(l0) = b(it,l1)*pinv
l0 = l0+1
410 continue
go to 470
420 l = 1
i = it-n10
do 460 j=1,k2
i = i+1
l0 = i
430 l1 = l0-k1
if(l1.le.0) go to 450
if(l1.le.n11) go to 440
l0 = l1-n11
go to 430
440 h1(l1) = b(it,j)*pinv
go to 460
450 h2(l0) = h2(l0)+b(it,j)*pinv
460 continue
if(n11.le.0) go to 510
c rotate this row into triangle by givens transformations
c rotation with the rows l,l+1,...n11.
470 do 500 j=l,n11
piv = h1(1)
c calculate the parameters of the givens transformation.
call fpgivs(piv,g1(j,1),cos,sin)
c transformation to right hand side.
j1 = j
do 475 j2=1,idim
call fprota(cos,sin,xi(j2),c(j1))
j1 = j1+n
475 continue
c transformation to the left hand side with respect to g2.
do 480 i=1,k1
call fprota(cos,sin,h2(i),g2(j,i))
480 continue
if(j.eq.n11) go to 510
i2 = min0(n11-j,k1)
c transformation to the left hand side with respect to g1.
do 490 i=1,i2
i1 = i+1
call fprota(cos,sin,h1(i1),g1(j,i1))
h1(i) = h1(i1)
490 continue
h1(i1) = 0.
500 continue
c rotation with the rows n11+1,...n7
510 do 530 j=1,k1
ij = n11+j
if(ij.le.0) go to 530
piv = h2(j)
c calculate the parameters of the givens transformation
call fpgivs(piv,g2(ij,j),cos,sin)
c transformation to the right hand side.
j1 = ij
do 515 j2=1,idim
call fprota(cos,sin,xi(j2),c(j1))
j1 = j1+n
515 continue
if(j.eq.k1) go to 540
j1 = j+1
c transformation to the left hand side.
do 520 i=j1,k1
call fprota(cos,sin,h2(i),g2(ij,i))
520 continue
530 continue
540 continue
c backward substitution to obtain the b-spline coefficients
j1 = 1
do 542 j2=1,idim
call fpbacp(g1,g2,c(j1),n7,k1,c(j1),k2,nest)
j1 = j1+n
542 continue
c calculate from condition (**) the remaining b-spline coefficients.
do 547 i=1,k
j1 = i
do 545 j=1,idim
j2 = j1+n7
c(j2) = c(j1)
j1 = j1+n
545 continue
547 continue
c computation of f(p).
fp = 0.
l = k1
jj = 0
do 570 it=1,m1
if(u(it).lt.t(l)) go to 550
l = l+1
550 l0 = l-k2
term = 0.
do 565 j2=1,idim
fac = 0.
j1 = l0
do 560 j=1,k1
j1 = j1+1
fac = fac+c(j1)*q(it,j)
560 continue
jj = jj+1
term = term+(fac-x(jj))**2
l0 = l0+n
565 continue
fp = fp+term*w(it)**2
570 continue
c test whether the approximation sp(u) is an acceptable solution.
fpms = fp-s
if(abs(fpms).lt.acc) go to 660
c test whether the maximal number of iterations is reached.
if(iter.eq.maxit) go to 600
c carry out one more step of the iteration process.
p2 = p
f2 = fpms
if(ich3.ne.0) go to 580
if((f2-f3) .gt. acc) go to 575
c our initial choice of p is too large.
p3 = p2
f3 = f2
p = p*con4
if(p.le.p1) p = p1*con9 +p2*con1
go to 595
575 if(f2.lt.0.) ich3 = 1
580 if(ich1.ne.0) go to 590
if((f1-f2) .gt. acc) go to 585
c our initial choice of p is too small
p1 = p2
f1 = f2
p = p/con4
if(p3.lt.0.) go to 595
if(p.ge.p3) p = p2*con1 +p3*con9
go to 595
585 if(f2.gt.0.) ich1 = 1
c test whether the iteration process proceeds as theoretically
c expected.
590 if(f2.ge.f1 .or. f2.le.f3) go to 610
c find the new value for p.
p = fprati(p1,f1,p2,f2,p3,f3)
595 continue
c error codes and messages.
600 ier = 3
go to 660
610 ier = 2
go to 660
620 ier = 1
go to 660
630 ier = -1
go to 660
640 ier = -2
c the point (z(1),z(2),...,z(idim)) is a solution of our problem.
c a constant function is a spline of degree k with all b-spline
c coefficients equal to that constant.
do 650 i=1,k1
rn = k1-i
t(i) = u(1)-rn*per
j = i+k1
rn = i-1
t(j) = u(m)+rn*per
650 continue
n = nmin
j1 = 0
do 658 j=1,idim
fac = z(j)
j2 = j1
do 654 i=1,k1
j2 = j2+1
c(j2) = fac
654 continue
j1 = j1+n
658 continue
fp = fp0
fpint(n) = fp0
fpint(n-1) = 0.
nrdata(n) = 0
660 return
end
|
