-
Notifications
You must be signed in to change notification settings - Fork 0
/
bslib1.f
4074 lines (4067 loc) · 149 KB
/
bslib1.f
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
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c
c.... B-Spline Routines from PPPack by c. de boor
c
c.... These are avalable on NetLib.org as public domain software
c
c.... Assembled by S. Collis
c
c
c.... Be sure to compile with -r8 otherwise you'll get roundoff errors for
c.... large nx
c
subroutine splint ( tau, gtau, t, n, k, q, bcoef, iflag )
c from * a practical guide to splines * by c. de boor
c calls bsplvb, banfac/slv
c
c splint produces the b-spline coeff.s bcoef of the spline of order
c k with knots t(i), i=1,..., n + k , which takes on the value
c gtau(i) at tau(i), i=1,..., n .
c
c****** i n p u t ******
c tau.....array of length n , containing data point abscissae.
c a s s u m p t i o n . . . tau is strictly increasing
c gtau.....corresponding array of length n , containing data point or-
c dinates
c t.....knot sequence, of length n+k
c n.....number of data points and dimension of spline space s(k,t)
c k.....order of spline
c
c****** o u t p u t ******
c q.....array of size (2*k-1)*n , containing the triangular factoriz-
c ation of the coefficient matrix of the linear system for the b-
c coefficients of the spline interpolant.
c the b-coeffs for the interpolant of an additional data set
c (tau(i),htau(i)), i=1,...,n with the same data abscissae can
c be obtained without going through all the calculations in this
c routine, simply by loading htau into bcoef and then execut-
c ing the call banslv ( q, 2*k-1, n, k-1, k-1, bcoef )
c bcoef.....the b-coefficients of the interpolant, of length n
c iflag.....an integer indicating success (= 1) or failure (= 2)
c the linear system to be solved is (theoretically) invertible if
c and only if
c t(i) .lt. tau(i) .lt. t(i+k), all i.
c violation of this condition is certain to lead to iflag = 2 .
c
c****** m e t h o d ******
c the i-th equation of the linear system a*bcoef = b for the b-co-
c effs of the interpolant enforces interpolation at tau(i), i=1,...,n.
c hence, b(i) = gtau(i), all i, and a is a band matrix with 2k-1
c bands (if it is invertible).
c the matrix a is generated row by row and stored, diagonal by di-
c agonal, in the r o w s of the array q , with the main diagonal go-
c ing into row k . see comments in the program below.
c the banded system is then solved by a call to banfac (which con-
c structs the triangular factorization for a and stores it again in
c q ), followed by a call to banslv (which then obtains the solution
c bcoef by substitution).
c banfac does no pivoting, since the total positivity of the matrix
c a makes this unnecessary.
c
integer iflag,k,n, i,ilp1mx,j,jj,km1,kpkm2,left,lenq,np1
real bcoef(n),gtau(n),q(1),t(1),tau(n), taui
c dimension q(2*k-1,n), t(n+k)
c current fortran standard makes it impossible to specify precisely the
c dimension of q and t without the introduction of otherwise super-
c fluous additional arguments.
np1 = n + 1
km1 = k - 1
kpkm2 = 2*km1
left = k
c zero out all entries of q
lenq = n*(k+km1)
do 5 i=1,lenq
5 q(i) = 0.
c
c *** loop over i to construct the n interpolation equations
do 30 i=1,n
taui = tau(i)
ilp1mx = min0(i+k,np1)
c *** find left in the closed interval (i,i+k-1) such that
c t(left) .le. tau(i) .lt. t(left+1)
c matrix is singular if this is not possible
left = max0(left,i)
if (taui .lt. t(left)) go to 998
15 if (taui .lt. t(left+1)) go to 16
left = left + 1
if (left .lt. ilp1mx) go to 15
left = left - 1
if (taui .gt. t(left+1)) go to 998
c *** the i-th equation enforces interpolation at taui, hence
c a(i,j) = b(j,k,t)(taui), all j. only the k entries with j =
c left-k+1,...,left actually might be nonzero. these k numbers
c are returned, in bcoef (used for temp.storage here), by the
c following
16 call bsplvb ( t, k, 1, taui, left, bcoef )
c we therefore want bcoef(j) = b(left-k+j)(taui) to go into
c a(i,left-k+j), i.e., into q(i-(left+j)+2*k,(left+j)-k) since
c a(i+j,j) is to go into q(i+k,j), all i,j, if we consider q
c as a two-dim. array , with 2*k-1 rows (see comments in
c banfac). in the present program, we treat q as an equivalent
c one-dimensional array (because of fortran restrictions on
c dimension statements) . we therefore want bcoef(j) to go into
c entry
c i -(left+j) + 2*k + ((left+j) - k-1)*(2*k-1)
c = i-left+1 + (left -k)*(2*k-1) + (2*k-2)*j
c of q .
jj = i-left+1 + (left-k)*(k+km1)
do 30 j=1,k
jj = jj+kpkm2
30 q(jj) = bcoef(j)
c
c ***obtain factorization of a , stored again in q.
call banfac ( q, k+km1, n, km1, km1, iflag )
go to (40,999), iflag
c *** solve a*bcoef = gtau by backsubstitution
40 do 41 i=1,n
41 bcoef(i) = gtau(i)
call banslv ( q, k+km1, n, km1, km1, bcoef )
return
998 iflag = 2
999 print 699
699 format(41h linear system in splint not invertible)
return
end
subroutine bsplvb ( t, jhigh, index, x, left, biatx )
c from * a practical guide to splines * by c. de boor
c calculates the value of all possibly nonzero b-splines at x of order
c
c jout = max( jhigh , (j+1)*(index-1) )
c
c with knot sequence t .
c
c****** i n p u t ******
c t.....knot sequence, of length left + jout , assumed to be nonde-
c creasing. a s s u m p t i o n . . . .
c t(left) .lt. t(left + 1) .
c d i v i s i o n b y z e r o will result if t(left) = t(left+1)
c jhigh,
c index.....integers which determine the order jout = max(jhigh,
c (j+1)*(index-1)) of the b-splines whose values at x are to
c be returned. index is used to avoid recalculations when seve-
c ral columns of the triangular array of b-spline values are nee-
c ded (e.g., in bsplpp or in bsplvd ). precisely,
c if index = 1 ,
c the calculation starts from scratch and the entire triangular
c array of b-spline values of orders 1,2,...,jhigh is generated
c order by order , i.e., column by column .
c if index = 2 ,
c only the b-spline values of order j+1, j+2, ..., jout are ge-
c nerated, the assumption being that biatx , j , deltal , deltar
c are, on entry, as they were on exit at the previous call.
c in particular, if jhigh = 0, then jout = j+1, i.e., just
c the next column of b-spline values is generated.
c
c w a r n i n g . . . the restriction jout .le. jmax (= 20) is im-
c posed arbitrarily by the dimension statement for deltal and
c deltar below, but is n o w h e r e c h e c k e d for .
c
c x.....the point at which the b-splines are to be evaluated.
c left.....an integer chosen (usually) so that
c t(left) .le. x .le. t(left+1) .
c
c****** o u t p u t ******
c biatx.....array of length jout , with biatx(i) containing the val-
c ue at x of the polynomial of order jout which agrees with
c the b-spline b(left-jout+i,jout,t) on the interval (t(left),
c t(left+1)) .
c
c****** m e t h o d ******
c the recurrence relation
c
c x - t(i) t(i+j+1) - x
c b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x)
c t(i+j)-t(i) t(i+j+1)-t(i+1)
c
c is used (repeatedly) to generate the (j+1)-vector b(left-j,j+1)(x),
c ...,b(left,j+1)(x) from the j-vector b(left-j+1,j)(x),...,
c b(left,j)(x), storing the new values in biatx over the old. the
c facts that
c b(i,1) = 1 if t(i) .le. x .lt. t(i+1)
c and that
c b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i+j)
c are used. the particular organization of the calculations follows al-
c gorithm (8) in chapter x of the text.
c
integer index,jhigh,left, i,j,jmax,jp1
parameter (jmax = 20)
real biatx(jhigh),t(1),x, deltal(jmax),deltar(jmax),saved,term
C real biatx(jhigh),t(1),x, deltal(20),deltar(20),saved,term
c dimension biatx(jout), t(left+jout)
c current fortran standard makes it impossible to specify the length of
c t and of biatx precisely without the introduction of otherwise
c superfluous additional arguments.
data j/1/
save j,deltal,deltar
c
go to (10,20), index
10 j = 1
biatx(1) = 1.
if (j .ge. jhigh) go to 99
c
20 jp1 = j + 1
deltar(j) = t(left+j) - x
deltal(j) = x - t(left+1-j)
saved = 0.
do 26 i=1,j
term = biatx(i)/(deltar(i) + deltal(jp1-i))
biatx(i) = saved + deltar(i)*term
26 saved = deltal(jp1-i)*term
biatx(jp1) = saved
j = jp1
if (j .lt. jhigh) go to 20
c
99 return
end
subroutine banfac ( w, nroww, nrow, nbandl, nbandu, iflag )
c from * a practical guide to splines * by c. de boor
c returns in w the lu-factorization (without pivoting) of the banded
c matrix a of order nrow with (nbandl + 1 + nbandu) bands or diag-
c onals in the work array w .
c
c****** i n p u t ******
c w.....work array of size (nroww,nrow) containing the interesting
c part of a banded matrix a , with the diagonals or bands of a
c stored in the rows of w , while columns of a correspond to
c columns of w . this is the storage mode used in linpack and
c results in efficient innermost loops.
c explicitly, a has nbandl bands below the diagonal
c + 1 (main) diagonal
c + nbandu bands above the diagonal
c and thus, with middle = nbandu + 1,
c a(i+j,j) is in w(i+middle,j) for i=-nbandu,...,nbandl
c j=1,...,nrow .
c for example, the interesting entries of a (1,2)-banded matrix
c of order 9 would appear in the first 1+1+2 = 4 rows of w
c as follows.
c
c
c
c
c
c all other entries of w not identified in this way with an en-
c try of a are never referenced .
c nroww.....row dimension of the work array w .
c must be .ge. nbandl + 1 + nbandu .
c nbandl.....number of bands of a below the main diagonal
c nbandu.....number of bands of a above the main diagonal .
c
c****** o u t p u t ******
c iflag.....integer indicating success( = 1) or failure ( = 2) .
c if iflag = 1, then
c w.....contains the lu-factorization of a into a unit lower triangu-
c lar matrix l and an upper triangular matrix u (both banded)
c and stored in customary fashion over the corresponding entries
c of a . this makes it possible to solve any particular linear
c system a*x = b for x by a
c call banslv ( w, nroww, nrow, nbandl, nbandu, b )
c with the solution x contained in b on return .
c if iflag = 2, then
c one of nrow-1, nbandl,nbandu failed to be nonnegative, or else
c one of the potential pivots was found to be zero indicating
c that a does not have an lu-factorization. this implies that
c a is singular in case it is totally positive .
c
c****** m e t h o d ******
c gauss elimination w i t h o u t pivoting is used. the routine is
c intended for use with matrices a which do not require row inter-
c changes during factorization, especially for the t o t a l l y
c p o s i t i v e matrices which occur in spline calculations.
c the routine should not be used for an arbitrary banded matrix.
c
integer iflag,nbandl,nbandu,nrow,nroww, i,ipk,j,jmax,k,kmax
* ,middle,midmk,nrowm1
real w(nroww,nrow), factor,pivot
c
iflag = 1
middle = nbandu + 1
c w(middle,.) contains the main diagonal of a .
nrowm1 = nrow - 1
if (nrowm1) 999,900,1
1 if (nbandl .gt. 0) go to 10
c a is upper triangular. check that diagonal is nonzero .
do 5 i=1,nrowm1
if (w(middle,i) .eq. 0.) go to 999
5 continue
go to 900
10 if (nbandu .gt. 0) go to 20
c a is lower triangular. check that diagonal is nonzero and
c divide each column by its diagonal .
do 15 i=1,nrowm1
pivot = w(middle,i)
if(pivot .eq. 0.) go to 999
jmax = min0(nbandl, nrow - i)
do 15 j=1,jmax
15 w(middle+j,i) = w(middle+j,i)/pivot
return
c
c a is not just a triangular matrix. construct lu factorization
20 do 50 i=1,nrowm1
c w(middle,i) is pivot for i-th step .
pivot = w(middle,i)
if (pivot .eq. 0.) go to 999
c jmax is the number of (nonzero) entries in column i
c below the diagonal .
jmax = min0(nbandl,nrow - i)
c divide each entry in column i below diagonal by pivot .
do 32 j=1,jmax
32 w(middle+j,i) = w(middle+j,i)/pivot
c kmax is the number of (nonzero) entries in row i to
c the right of the diagonal .
kmax = min0(nbandu,nrow - i)
c subtract a(i,i+k)*(i-th column) from (i+k)-th column
c (below row i ) .
do 40 k=1,kmax
ipk = i + k
midmk = middle - k
factor = w(midmk,ipk)
do 40 j=1,jmax
40 w(midmk+j,ipk) = w(midmk+j,ipk) - w(middle+j,i)*factor
50 continue
c check the last diagonal entry .
900 if (w(middle,nrow) .ne. 0.) return
999 iflag = 2
return
end
subroutine banslv ( w, nroww, nrow, nbandl, nbandu, b )
c from * a practical guide to splines * by c. de boor
c companion routine to banfac . it returns the solution x of the
c linear system a*x = b in place of b , given the lu-factorization
c for a in the workarray w .
c
c****** i n p u t ******
c w, nroww,nrow,nbandl,nbandu.....describe the lu-factorization of a
c banded matrix a of order nrow as constructed in banfac .
c for details, see banfac .
c b.....right side of the system to be solved .
c
c****** o u t p u t ******
c b.....contains the solution x , of order nrow .
c
c****** m e t h o d ******
c (with a = l*u, as stored in w,) the unit lower triangular system
c l(u*x) = b is solved for y = u*x, and y stored in b . then the
c upper triangular system u*x = y is solved for x . the calcul-
c ations are so arranged that the innermost loops stay within columns.
c
integer nbandl,nbandu,nrow,nroww, i,j,jmax,middle,nrowm1
real w(nroww,nrow),b(nrow)
middle = nbandu + 1
if (nrow .eq. 1) go to 49
nrowm1 = nrow - 1
if (nbandl .eq. 0) go to 30
c forward pass
c for i=1,2,...,nrow-1, subtract right side(i)*(i-th column
c of l ) from right side (below i-th row) .
do 21 i=1,nrowm1
jmax = min0(nbandl, nrow-i)
do 21 j=1,jmax
21 b(i+j) = b(i+j) - b(i)*w(middle+j,i)
c backward pass
c for i=nrow,nrow-1,...,1, divide right side(i) by i-th diag-
c onal entry of u, then subtract right side(i)*(i-th column
c of u) from right side (above i-th row).
30 if (nbandu .gt. 0) go to 40
c a is lower triangular .
do 31 i=1,nrow
31 b(i) = b(i)/w(1,i)
return
40 i = nrow
41 b(i) = b(i)/w(middle,i)
jmax = min0(nbandu,i-1)
do 45 j=1,jmax
45 b(i-j) = b(i-j) - b(i)*w(middle-j,i)
i = i - 1
if (i .gt. 1) go to 41
49 b(1) = b(1)/w(middle,1)
return
end
real function bvalue ( t, bcoef, n, k, x, jderiv )
c from * a practical guide to splines * by c. de boor
c calls interv
c
c calculates value at x of jderiv-th derivative of spline from b-repr.
c the spline is taken to be continuous from the right, EXCEPT at the
c rightmost knot, where it is taken to be continuous from the left.
c
c****** i n p u t ******
c t, bcoef, n, k......forms the b-representation of the spline f to
c be evaluated. specifically,
c t.....knot sequence, of length n+k, assumed nondecreasing.
c bcoef.....b-coefficient sequence, of length n .
c n.....length of bcoef and dimension of spline(k,t),
c a s s u m e d positive .
c k.....order of the spline .
c
c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed
c arbitrarily by the dimension statement for aj, dl, dr below,
c but is n o w h e r e c h e c k e d for.
c
c x.....the point at which to evaluate .
c jderiv.....integer giving the order of the derivative to be evaluated
c a s s u m e d to be zero or positive.
c
c****** o u t p u t ******
c bvalue.....the value of the (jderiv)-th derivative of f at x .
c
c****** m e t h o d ******
c The nontrivial knot interval (t(i),t(i+1)) containing x is lo-
c cated with the aid of interv . The k b-coeffs of f relevant for
c this interval are then obtained from bcoef (or taken to be zero if
c not explicitly available) and are then differenced jderiv times to
c obtain the b-coeffs of (d**jderiv)f relevant for that interval.
c Precisely, with j = jderiv, we have from x.(12) of the text that
c
c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
c
c where
c / bcoef(.), , j .eq. 0
c /
c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1)
c / ----------------------------- , j .gt. 0
c / (t(.+k-j) - t(.))/(k-j)
c
c Then, we use repeatedly the fact that
c
c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) )
c with
c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
c a(.,x) = ---------------------------------------
c (x - t(.)) + (t(.+m-1) - x)
c
c to write (d**j)f(x) eventually as a linear combination of b-splines
c of order 1 , and the coefficient for b(i,1,t)(x) must then be the
c desired number (d**j)f(x). (see x.(17)-(19) of text).
c
integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmax,kmj,km1
* ,mflag,nmi,jdrvp1
parameter (kmax = 20)
C real bcoef(n),t(1),x, aj(20),dl(20),dr(20),fkmj
real bcoef(n),t(1),x, aj(kmax),dl(kmax),dr(kmax),fkmj
c dimension t(n+k)
c former fortran standard made it impossible to specify the length of t
c precisely without the introduction of otherwise superfluous addition-
c al arguments.
bvalue = 0.
if (jderiv .ge. k) go to 99
c
c *** Find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
c t(i) .le. x .lt. t(i+1) . If no such i can be found, x lies
c outside the support of the spline f , hence bvalue = 0.
c (The asymmetry in this choice of i makes f rightcontinuous, except
c at t(n+k) where it is leftcontinuous.)
call interv ( t, n+k, x, i, mflag )
if (mflag .ne. 0) go to 99
c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
km1 = k - 1
if (km1 .gt. 0) go to 1
bvalue = bcoef(i)
go to 99
c
c *** store the k b-spline coefficients relevant for the knot interval
c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
c from input to zero. set any t.s not obtainable equal to t(1) or
c to t(n+k) appropriately.
1 jcmin = 1
imk = i - k
if (imk .ge. 0) go to 8
jcmin = 1 - imk
do 5 j=1,i
5 dl(j) = x - t(i+1-j)
do 6 j=i,km1
aj(k-j) = 0.
6 dl(j) = dl(i)
go to 10
8 do 9 j=1,km1
9 dl(j) = x - t(i+1-j)
c
10 jcmax = k
nmi = n - i
if (nmi .ge. 0) go to 18
jcmax = k + nmi
do 15 j=1,jcmax
15 dr(j) = t(i+j) - x
do 16 j=jcmax,km1
aj(j+1) = 0.
16 dr(j) = dr(jcmax)
go to 20
18 do 19 j=1,km1
19 dr(j) = t(i+j) - x
c
20 do 21 jc=jcmin,jcmax
21 aj(jc) = bcoef(imk + jc)
c
c *** difference the coefficients jderiv times.
if (jderiv .eq. 0) go to 30
do 23 j=1,jderiv
kmj = k-j
fkmj = float(kmj)
ilo = kmj
do 23 jj=1,kmj
aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
23 ilo = ilo - 1
c
c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative,
c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
30 if (jderiv .eq. km1) go to 39
jdrvp1 = jderiv + 1
do 33 j=jdrvp1,km1
kmj = k-j
ilo = kmj
do 33 jj=1,kmj
aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
33 ilo = ilo - 1
39 bvalue = aj(1)
c
99 return
end
subroutine interv ( xt, lxt, x, left, mflag )
c from * a practical guide to splines * by C. de Boor
c computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) .
c
c****** i n p u t ******
c xt.....a real sequence, of length lxt , assumed to be nondecreasing
c lxt.....number of terms in the sequence xt .
c x.....the point whose location with respect to the sequence xt is
c to be determined.
c
c****** o u t p u t ******
c left, mflag.....both integers, whose value is
c
c 1 -1 if x .lt. xt(1)
c i 0 if xt(i) .le. x .lt. xt(i+1)
c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt)
c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x
c
c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0
c indicates that x lies outside the CLOSED interval
c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the
c intervals is due to the decision to make all pp functions cont-
c inuous from the right, but, by returning mflag = 0 even if
C x = xt(lxt), there is the option of having the computed pp function
c continuous from the left at xt(lxt) .
c
c****** m e t h o d ******
c The program is designed to be efficient in the common situation that
c it is called repeatedly, with x taken from an increasing or decrea-
c sing sequence. This will happen, e.g., when a pp function is to be
c graphed. The first guess for left is therefore taken to be the val-
c ue returned at the previous call and stored in the l o c a l varia-
c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec-
c essary since the present call may have nothing to do with the previ-
c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left =
c ilo and are done after just three comparisons.
c Otherwise, we repeatedly double the difference istep = ihi - ilo
c while also moving ilo and ihi in the direction of x , until
c xt(ilo) .le. x .lt. xt(ihi) ,
c after which we use bisection to get, in addition, ilo+1 = ihi .
c left = ilo is then returned.
c
integer left,lxt,mflag, ihi,ilo,istep,middle
real x,xt(lxt)
data ilo /1/
save ilo
ihi = ilo + 1
if (ihi .lt. lxt) go to 20
if (x .ge. xt(lxt)) go to 110
if (lxt .le. 1) go to 90
ilo = lxt - 1
ihi = lxt
c
20 if (x .ge. xt(ihi)) go to 40
if (x .ge. xt(ilo)) go to 100
c
c **** now x .lt. xt(ilo) . decrease ilo to capture x .
istep = 1
31 ihi = ilo
ilo = ihi - istep
if (ilo .le. 1) go to 35
if (x .ge. xt(ilo)) go to 50
istep = istep*2
go to 31
35 ilo = 1
if (x .lt. xt(1)) go to 90
go to 50
c **** now x .ge. xt(ihi) . increase ihi to capture x .
40 istep = 1
41 ilo = ihi
ihi = ilo + istep
if (ihi .ge. lxt) go to 45
if (x .lt. xt(ihi)) go to 50
istep = istep*2
go to 41
45 if (x .ge. xt(lxt)) go to 110
ihi = lxt
c
c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
50 middle = (ilo + ihi)/2
if (middle .eq. ilo) go to 100
c note. it is assumed that middle = ilo in case ihi = ilo+1 .
if (x .lt. xt(middle)) go to 53
ilo = middle
go to 50
53 ihi = middle
go to 50
c**** set output and return.
90 mflag = -1
left = 1
return
100 mflag = 0
left = ilo
return
110 mflag = 1
if (x .eq. xt(lxt)) mflag = 0
left = lxt
111 if (left .eq. 1) return
left = left - 1
if (xt(left) .lt. xt(lxt)) return
go to 111
end
*DECK BINTK
SUBROUTINE BINTK (X, Y, T, N, K, BCOEF, Q, WORK)
C***BEGIN PROLOGUE BINTK
C***PURPOSE Compute the B-representation of a spline which interpolates
C given data.
C***LIBRARY SLATEC
C***CATEGORY E1A
C***TYPE SINGLE PRECISION (BINTK-S, DBINTK-D)
C***KEYWORDS B-SPLINE, DATA FITTING, INTERPOLATION
C***AUTHOR Amos, D. E., (SNLA)
C***DESCRIPTION
C
C Written by Carl de Boor and modified by D. E. Amos
C
C Abstract
C
C BINTK is the SPLINT routine of the reference.
C
C BINTK produces the B-spline coefficients, BCOEF, of the
C B-spline of order K with knots T(I), I=1,...,N+K, which
C takes on the value Y(I) at X(I), I=1,...,N. The spline or
C any of its derivatives can be evaluated by calls to BVALU.
C The I-th equation of the linear system A*BCOEF = B for the
C coefficients of the interpolant enforces interpolation at
C X(I)), I=1,...,N. Hence, B(I) = Y(I), all I, and A is
C a band matrix with 2K-1 bands if A is invertible. The matrix
C A is generated row by row and stored, diagonal by diagonal,
C in the rows of Q, with the main diagonal going into row K.
C The banded system is then solved by a call to BNFAC (which
C constructs the triangular factorization for A and stores it
C again in Q), followed by a call to BNSLV (which then
C obtains the solution BCOEF by substitution). BNFAC does no
C pivoting, since the total positivity of the matrix A makes
C this unnecessary. The linear system to be solved is
C (theoretically) invertible if and only if
C T(I) .LT. X(I)) .LT. T(I+K), all I.
C Equality is permitted on the left for I=1 and on the right
C for I=N when K knots are used at X(1) or X(N). Otherwise,
C violation of this condition is certain to lead to an error.
C
C Description of Arguments
C Input
C X - vector of length N containing data point abscissa
C in strictly increasing order.
C Y - corresponding vector of length N containing data
C point ordinates.
C T - knot vector of length N+K
C since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K)
C .GE. X(N), this leaves only N-K knots (not nec-
C essarily X(I)) values) interior to (X(1),X(N))
C N - number of data points, N .GE. K
C K - order of the spline, K .GE. 1
C
C Output
C BCOEF - a vector of length N containing the B-spline
C coefficients
C Q - a work vector of length (2*K-1)*N, containing
C the triangular factorization of the coefficient
C matrix of the linear system being solved. The
C coefficients for the interpolant of an
C additional data set (X(I)),YY(I)), I=1,...,N
C with the same abscissa can be obtained by loading
C YY into BCOEF and then executing
C CALL BNSLV (Q,2K-1,N,K-1,K-1,BCOEF)
C WORK - work vector of length 2*K
C
C Error Conditions
C Improper input is a fatal error
C Singular system of equations is a fatal error
C
C***REFERENCES D. E. Amos, Computation with splines and B-splines,
C Report SAND78-1968, Sandia Laboratories, March 1979.
C Carl de Boor, Package for calculating with B-splines,
C SIAM Journal on Numerical Analysis 14, 3 (June 1977),
C pp. 441-472.
C Carl de Boor, A Practical Guide to Splines, Applied
C Mathematics Series 27, Springer-Verlag, New York,
C 1978.
C***ROUTINES CALLED BNFAC, BNSLV, BSPVN, XERMSG
C***REVISION HISTORY (YYMMDD)
C 800901 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE BINTK
C
INTEGER IFLAG, IWORK, K, N, I, ILP1MX, J, JJ, KM1, KPKM2, LEFT,
1 LENQ, NP1
REAL BCOEF(*), Y(*), Q(*), T(*), X(*), XI, WORK(*)
C DIMENSION Q(2*K-1,N), T(N+K)
C***FIRST EXECUTABLE STATEMENT BINTK
IF(K.LT.1) GO TO 100
IF(N.LT.K) GO TO 105
JJ = N - 1
IF(JJ.EQ.0) GO TO 6
DO 5 I=1,JJ
IF(X(I).GE.X(I+1)) GO TO 110
5 CONTINUE
6 CONTINUE
NP1 = N + 1
KM1 = K - 1
KPKM2 = 2*KM1
LEFT = K
C ZERO OUT ALL ENTRIES OF Q
LENQ = N*(K+KM1)
DO 10 I=1,LENQ
Q(I) = 0.0E0
10 CONTINUE
C
C *** LOOP OVER I TO CONSTRUCT THE N INTERPOLATION EQUATIONS
DO 50 I=1,N
XI = X(I)
ILP1MX = MIN(I+K,NP1)
C *** FIND LEFT IN THE CLOSED INTERVAL (I,I+K-1) SUCH THAT
C T(LEFT) .LE. X(I) .LT. T(LEFT+1)
C MATRIX IS SINGULAR IF THIS IS NOT POSSIBLE
LEFT = MAX(LEFT,I)
IF (XI.LT.T(LEFT)) GO TO 80
20 IF (XI.LT.T(LEFT+1)) GO TO 30
LEFT = LEFT + 1
IF (LEFT.LT.ILP1MX) GO TO 20
LEFT = LEFT - 1
IF (XI.GT.T(LEFT+1)) GO TO 80
C *** THE I-TH EQUATION ENFORCES INTERPOLATION AT XI, HENCE
C A(I,J) = B(J,K,T)(XI), ALL J. ONLY THE K ENTRIES WITH J =
C LEFT-K+1,...,LEFT ACTUALLY MIGHT BE NONZERO. THESE K NUMBERS
C ARE RETURNED, IN BCOEF (USED FOR TEMP. STORAGE HERE), BY THE
C FOLLOWING
30 CALL BSPVN(T, K, K, 1, XI, LEFT, BCOEF, WORK, IWORK)
C WE THEREFORE WANT BCOEF(J) = B(LEFT-K+J)(XI) TO GO INTO
C A(I,LEFT-K+J), I.E., INTO Q(I-(LEFT+J)+2*K,(LEFT+J)-K) SINCE
C A(I+J,J) IS TO GO INTO Q(I+K,J), ALL I,J, IF WE CONSIDER Q
C AS A TWO-DIM. ARRAY , WITH 2*K-1 ROWS (SEE COMMENTS IN
C BNFAC). IN THE PRESENT PROGRAM, WE TREAT Q AS AN EQUIVALENT
C ONE-DIMENSIONAL ARRAY (BECAUSE OF FORTRAN RESTRICTIONS ON
C DIMENSION STATEMENTS) . WE THEREFORE WANT BCOEF(J) TO GO INTO
C ENTRY
C I -(LEFT+J) + 2*K + ((LEFT+J) - K-1)*(2*K-1)
C = I-LEFT+1 + (LEFT -K)*(2*K-1) + (2*K-2)*J
C OF Q .
JJ = I - LEFT + 1 + (LEFT-K)*(K+KM1)
DO 40 J=1,K
JJ = JJ + KPKM2
Q(JJ) = BCOEF(J)
40 CONTINUE
50 CONTINUE
C
C ***OBTAIN FACTORIZATION OF A , STORED AGAIN IN Q.
CALL BNFAC(Q, K+KM1, N, KM1, KM1, IFLAG)
GO TO (60, 90), IFLAG
C *** SOLVE A*BCOEF = Y BY BACKSUBSTITUTION
60 DO 70 I=1,N
BCOEF(I) = Y(I)
70 CONTINUE
CALL BNSLV(Q, K+KM1, N, KM1, KM1, BCOEF)
RETURN
C
C
80 CONTINUE
CALL XERMSG ('SLATEC', 'BINTK',
+ 'SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ' //
+ 'BASIS FUNCTION AND THE SYSTEM IS SINGULAR.', 2, 1)
RETURN
90 CONTINUE
CALL XERMSG ('SLATEC', 'BINTK',
+ 'THE SYSTEM OF SOLVER DETECTS A SINGULAR SYSTEM ALTHOUGH ' //
+ 'THE THEORETICAL CONDITIONS FOR A SOLUTION WERE SATISFIED.',
+ 8, 1)
RETURN
100 CONTINUE
CALL XERMSG ('SLATEC', 'BINTK', 'K DOES NOT SATISFY K.GE.1', 2,
+ 1)
RETURN
105 CONTINUE
CALL XERMSG ('SLATEC', 'BINTK', 'N DOES NOT SATISFY N.GE.K', 2,
+ 1)
RETURN
110 CONTINUE
CALL XERMSG ('SLATEC', 'BINTK',
+ 'X(I) DOES NOT SATISFY X(I).LT.X(I+1) FOR SOME I', 2, 1)
RETURN
END
*DECK BNFAC
SUBROUTINE BNFAC (W, NROWW, NROW, NBANDL, NBANDU, IFLAG)
C***BEGIN PROLOGUE BNFAC
C***SUBSIDIARY
C***PURPOSE Subsidiary to BINT4 and BINTK
C***LIBRARY SLATEC
C***TYPE SINGLE PRECISION (BNFAC-S, DBNFAC-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C BNFAC is the BANFAC routine from
C * A Practical Guide to Splines * by C. de Boor
C
C Returns in W the lu-factorization (without pivoting) of the banded
C matrix A of order NROW with (NBANDL + 1 + NBANDU) bands or diag-
C onals in the work array W .
C
C ***** I N P U T ******
C W.....Work array of size (NROWW,NROW) containing the interesting
C part of a banded matrix A , with the diagonals or bands of A
C stored in the rows of W , while columns of A correspond to
C columns of W . This is the storage mode used in LINPACK and
C results in efficient innermost loops.
C Explicitly, A has NBANDL bands below the diagonal
C + 1 (main) diagonal
C + NBANDU bands above the diagonal
C and thus, with MIDDLE = NBANDU + 1,
C A(I+J,J) is in W(I+MIDDLE,J) for I=-NBANDU,...,NBANDL
C J=1,...,NROW .
C For example, the interesting entries of A (1,2)-banded matrix
C of order 9 would appear in the first 1+1+2 = 4 rows of W
C as follows.
C 13 24 35 46 57 68 79
C 12 23 34 45 56 67 78 89
C 11 22 33 44 55 66 77 88 99
C 21 32 43 54 65 76 87 98
C
C All other entries of W not identified in this way with an en-
C try of A are never referenced .
C NROWW.....Row dimension of the work array W .
C must be .GE. NBANDL + 1 + NBANDU .
C NBANDL.....Number of bands of A below the main diagonal
C NBANDU.....Number of bands of A above the main diagonal .
C
C ***** O U T P U T ******
C IFLAG.....Integer indicating success( = 1) or failure ( = 2) .
C If IFLAG = 1, then
C W.....contains the LU-factorization of A into a unit lower triangu-
C lar matrix L and an upper triangular matrix U (both banded)
C and stored in customary fashion over the corresponding entries
C of A . This makes it possible to solve any particular linear
C system A*X = B for X by A
C CALL BNSLV ( W, NROWW, NROW, NBANDL, NBANDU, B )
C with the solution X contained in B on return .
C If IFLAG = 2, then
C one of NROW-1, NBANDL,NBANDU failed to be nonnegative, or else
C one of the potential pivots was found to be zero indicating
C that A does not have an LU-factorization. This implies that
C A is singular in case it is totally positive .
C
C ***** M E T H O D ******
C Gauss elimination W I T H O U T pivoting is used. The routine is
C intended for use with matrices A which do not require row inter-
C changes during factorization, especially for the T O T A L L Y
C P O S I T I V E matrices which occur in spline calculations.
C The routine should not be used for an arbitrary banded matrix.
C
C***SEE ALSO BINT4, BINTK
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 800901 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE BNFAC
C
INTEGER IFLAG, NBANDL, NBANDU, NROW, NROWW, I, IPK, J, JMAX, K,
1 KMAX, MIDDLE, MIDMK, NROWM1
REAL W(NROWW,*), FACTOR, PIVOT
C
C***FIRST EXECUTABLE STATEMENT BNFAC
IFLAG = 1
MIDDLE = NBANDU + 1
C W(MIDDLE,.) CONTAINS THE MAIN DIAGONAL OF A .
NROWM1 = NROW - 1
IF (NROWM1) 120, 110, 10
10 IF (NBANDL.GT.0) GO TO 30
C A IS UPPER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO .
DO 20 I=1,NROWM1
IF (W(MIDDLE,I).EQ.0.0E0) GO TO 120
20 CONTINUE
GO TO 110
30 IF (NBANDU.GT.0) GO TO 60
C A IS LOWER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO AND
C DIVIDE EACH COLUMN BY ITS DIAGONAL .
DO 50 I=1,NROWM1
PIVOT = W(MIDDLE,I)
IF (PIVOT.EQ.0.0E0) GO TO 120
JMAX = MIN(NBANDL,NROW-I)
DO 40 J=1,JMAX
W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT
40 CONTINUE
50 CONTINUE
RETURN
C
C A IS NOT JUST A TRIANGULAR MATRIX. CONSTRUCT LU FACTORIZATION
60 DO 100 I=1,NROWM1
C W(MIDDLE,I) IS PIVOT FOR I-TH STEP .
PIVOT = W(MIDDLE,I)
IF (PIVOT.EQ.0.0E0) GO TO 120
C JMAX IS THE NUMBER OF (NONZERO) ENTRIES IN COLUMN I
C BELOW THE DIAGONAL .
JMAX = MIN(NBANDL,NROW-I)
C DIVIDE EACH ENTRY IN COLUMN I BELOW DIAGONAL BY PIVOT .
DO 70 J=1,JMAX
W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT
70 CONTINUE
C KMAX IS THE NUMBER OF (NONZERO) ENTRIES IN ROW I TO
C THE RIGHT OF THE DIAGONAL .
KMAX = MIN(NBANDU,NROW-I)
C SUBTRACT A(I,I+K)*(I-TH COLUMN) FROM (I+K)-TH COLUMN
C (BELOW ROW I ) .
DO 90 K=1,KMAX
IPK = I + K
MIDMK = MIDDLE - K
FACTOR = W(MIDMK,IPK)
DO 80 J=1,JMAX
W(MIDMK+J,IPK) = W(MIDMK+J,IPK) - W(MIDDLE+J,I)*FACTOR
80 CONTINUE
90 CONTINUE
100 CONTINUE
C CHECK THE LAST DIAGONAL ENTRY .
110 IF (W(MIDDLE,NROW).NE.0.0E0) RETURN
120 IFLAG = 2
RETURN
END
*DECK BNSLV
SUBROUTINE BNSLV (W, NROWW, NROW, NBANDL, NBANDU, B)
C***BEGIN PROLOGUE BNSLV
C***SUBSIDIARY
C***PURPOSE Subsidiary to BINT4 and BINTK
C***LIBRARY SLATEC
C***TYPE SINGLE PRECISION (BNSLV-S, DBNSLV-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C BNSLV is the BANSLV routine from
C * A Practical Guide to Splines * by C. de Boor
C
C Companion routine to BNFAC . It returns the solution X of the
C linear system A*X = B in place of B , given the LU-factorization
C for A in the work array W from BNFAC.
C
C ***** I N P U T ******
C W, NROWW,NROW,NBANDL,NBANDU.....Describe the LU-factorization of a
C banded matrix A of order NROW as constructed in BNFAC .
C For details, see BNFAC .
C B.....Right side of the system to be solved .
C
C ***** O U T P U T ******
C B.....Contains the solution X , of order NROW .
C
C ***** M E T H O D ******
C (With A = L*U, as stored in W,) the unit lower triangular system
C L(U*X) = B is solved for Y = U*X, and Y stored in B . Then the
C upper triangular system U*X = Y is solved for X . The calcul-
C ations are so arranged that the innermost loops stay within columns.
C
C***SEE ALSO BINT4, BINTK
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 800901 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE BNSLV
C
INTEGER NBANDL, NBANDU, NROW, NROWW, I, J, JMAX, MIDDLE, NROWM1
REAL W(NROWW,*), B(*)
C***FIRST EXECUTABLE STATEMENT BNSLV
MIDDLE = NBANDU + 1
IF (NROW.EQ.1) GO TO 80
NROWM1 = NROW - 1
IF (NBANDL.EQ.0) GO TO 30
C FORWARD PASS
C FOR I=1,2,...,NROW-1, SUBTRACT RIGHT SIDE(I)*(I-TH COLUMN