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
|
c prog5
c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12
c to appear in 'solving least squares problems', prentice-hall, 1974
c example of the use of subroutines bndacc and bndsol to solve
c sequentially the banded least squares problem which arises in
c spline curve fitting.
c
dimension x(12),y(12),b(10),g(12,5),c(12),q(4),cov(12)
c
c define functions p1 and p2 to be used in constructing
c cubic splines over uniformly spaced breakpoints.
c
p1(t)=.25*t**2*t
p2(t)=-(1.-t)**2*(1.+t)*.75+1.
zero=0.
c
write (6,230)
mdg=12
nband=4
m=12
c set ordinates of data
y(1)=2.2
y(2)=4.0
y(3)=5.0
y(4)=4.6
y(5)=2.8
y(6)=2.7
y(7)=3.8
y(8)=5.1
y(9)=6.1
y(10)=6.3
y(11)=5.0
y(12)=2.0
c set abcissas of data
do 10 i=1,m
10 x(i)=2*i
c
c begin loop thru cases using increasing nos of breakpoints.
c
do 150 nbp=5,10
nc=nbp+2
c set breakpoints
b(1)=x(1)
b(nbp)=x(m)
h=(b(nbp)-b(1))/float(nbp-1)
if (nbp.le.2) go to 30
do 20 i=3,nbp
20 b(i-1)=b(i-2)+h
30 continue
write (6,240) nbp,(b(i),i=1,nbp)
c
c initialize ir and ip before first call to bndacc.
c
ir=1
ip=1
i=1
jt=1
40 mt=0
50 continue
if (x(i).gt.b(jt+1)) go to 60
c
c set row for ith data point
c
u=(x(i)-b(jt))/h
ig=ir+mt
g(ig,1)=p1(1.-u)
g(ig,2)=p2(1.-u)
g(ig,3)=p2(u)
g(ig,4)=p1(u)
g(ig,5)=y(i)
mt=mt+1
if (i.eq.m) go to 60
i=i+1
go to 50
c
c send block of data to processor
c
60 continue
call bndacc (g,mdg,nband,ip,ir,mt,jt)
if (i.eq.m) go to 70
jt=jt+1
go to 40
c compute solution c()
70 continue
call bndsol (1,g,mdg,nband,ip,ir,c,nc,rnorm)
c write solution coefficients
write (6,180) (c(l),l=1,nc)
write (6,210) rnorm
c
c compute and print x,y,yfit,r=y-yfit
c
write (6,160)
jt=1
do 110 i=1,m
80 if (x(i).le.b(jt+1)) go to 90
jt=jt+1
go to 80
c
90 u=(x(i)-b(jt))/h
q(1)=p1(1.-u)
q(2)=p2(1.-u)
q(3)=p2(u)
q(4)=p1(u)
yfit=zero
do 100 l=1,4
ic=jt-1+l
100 yfit=yfit+c(ic)*q(l)
r=y(i)-yfit
write (6,170) i,x(i),y(i),yfit,r
110 continue
c
c compute residual vector norm.
c
if (m.le.nc) go to 150
sigsq=(rnorm**2)/float(m-nc)
sigfac=sqrt(sigsq)
write (6,220) sigfac
write (6,200)
c
c compute and print cols. of covariance.
c
do 140 j=1,nc
do 120 i=1,nc
120 cov(i)=zero
cov(j)=1.
call bndsol (2,g,mdg,nband,ip,ir,cov,nc,rdummy)
call bndsol (3,g,mdg,nband,ip,ir,cov,nc,rdummy)
c
c compute the jth col. of the covariance matrix.
do 130 i=1,nc
130 cov(i)=cov(i)*sigsq
140 write (6,190) (l,j,cov(l),l=1,nc)
150 continue
stop
160 format (4h0 i,8x,1hx,10x,1hy,6x,4hyfit,4x,8hr=y-yfit/1x)
170 format (1x,i3,4x,f6.0,4x,f6.2,4x,f6.2,4x,f8.4)
180 format (4h0c =,10f10.5/(4x,10f10.5))
190 format (4(02x,2i5,e15.7))
200 format (46h0covariance matrix of the spline coefficients.)
210 format (9h0rnorm =,e15.8)
220 format (9h0sigfac =,e15.8)
230 format (50h1prog5. execute a sequence of cubic spline fits,27h
1to a discrete set of data.)
240 format (1x////,11h0new case..,/47h0number of breakpoints, includin
1g endpoints, is,i5/14h0breakpoints =,10f10.5,/(14x,10f10.5))
end
|
