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
|
c prog1
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 demonstrate algorithms hft and hs1 for solving least squares
c problems and algorithm cov for omputing the associated covariance
c matrices.
c
dimension a(8,8),h(8),b(8)
real gen,anoise
double precision sm
data mda/8/
c
do 180 noise=1,2
anoise=0.
if (noise.eq.2) anoise=1.e-4
write (6,230)
write (6,240) anoise
c initialize the data generation function
c ..
dummy=gen(-1.)
do 180 mn1=1,6,5
mn2=mn1+2
do 180 m=mn1,mn2
do 180 n=mn1,mn2
np1=n+1
write (6,250) m,n
c generate data
c ..
do 10 i=1,m
do 10 j=1,n
10 a(i,j)=gen(anoise)
do 20 i=1,m
20 b(i)=gen(anoise)
if(m .lt. n) go to 180
c
c ****** begin algorithm hft ******
c ..
do 30 j=1,n
30 call h12 (1,j,j+1,m,a(1,j),1,h(j),a(1,j+1),1,mda,n-j)
c ..
c the algorithm 'hft' is completed.
c
c ****** begin algorithm hs1 ******
c apply the transformations q(n)...q(1)=q to b
c replacing the previous contents of the array, b .
c ..
do 40 j=1,n
40 call h12 (2,j,j+1,m,a(1,j),1,h(j),b,1,1,1)
c solve the triangular system for the solution x.
c store x in the array, b .
c ..
do 80 k=1,n
i=np1-k
sm=0.d0
if (i.eq.n) go to 60
ip1=i+1
do 50 j=ip1,n
50 sm=sm+a(i,j)*dble(b(j))
60 if (a(i,i)) 80,70,80
70 write (6,260)
go to 180
80 b(i)=(b(i)-sm)/a(i,i)
c compute length of residual vector.
c ..
srsmsq=0.
if (n.eq.m) go to 100
mmn=m-n
do 90 j=1,mmn
npj=n+j
90 srsmsq=srsmsq+b(npj)**2
srsmsq=sqrt(srsmsq)
c ****** begin algorithm cov ******
c compute unscaled covariance matrix ((a**t)*a)**(-1)
c ..
100 do 110 j=1,n
110 a(j,j)=1./a(j,j)
if (n.eq.1) go to 140
nm1=n-1
do 130 i=1,nm1
ip1=i+1
do 130 j=ip1,n
jm1=j-1
sm=0.d0
do 120 l=i,jm1
120 sm=sm+a(i,l)*dble(a(l,j))
130 a(i,j)=-sm*a(j,j)
c ..
c the upper triangle of a has been inverted
c upon itself.
140 do 160 i=1,n
do 160 j=i,n
sm=0.d0
do 150 l=j,n
150 sm=sm+a(i,l)*dble(a(j,l))
160 a(i,j)=sm
c ..
c the upper triangular part of the
c symmetric matrix (a**t*a)**(-1) has
c replaced the upper triangular part of
c the a array.
write (6,200) (i,b(i),i=1,n)
write (6,190) srsmsq
write (6,210)
do 170 i=1,n
170 write (6,220) (i,j,a(i,j),j=i,n)
180 continue
stop
190 format (1h0,8x,17hresidual length =,e12.4)
200 format (1h0,8x,34hestimated parameters, x=a**(+)*b,,22h computed
1by 'hft,hs1'//(9x,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8,i6,e16.8))
210 format (1h0,8x,31hcovariance matrix (unscaled) of,22h estimated pa
1rameters.,19h computed by 'cov'./1x)
220 format (9x,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8,2i3,e16.8)
230 format (52h1 prog1. this program demonstrates the algorithms,19
1h hft, hs1, and cov.//40h caution.. 'prog1' does no checking for ,
252hnear rank deficient matrices. results in such cases,20h may be
3 meaningless.,/34h such cases are handled by 'prog2',11h or 'pro
4')
240 format (1h0,54hthe relative noise level of the generated data will
1 be,e16.4)
250 format (1h0////9h0 m n/1x,2i4)
260 format (1h0,8x,36h****** terminating this case due to,37h a divis
1or being exactly zero. ******)
end
|
