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
|
c prog2
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 algorithm hfti for solving least squares problems
c and algorithm cov for computing the associated unscaled
c covariance matrix.
c
dimension a(8,8),h(8),b(8),g(8)
real gen,anoise
integer ip(8)
double precision sm
data mda /8/
c
do 180 noise=1,2
anorm=500.
anoise=0.
tau=.5
if (noise.eq.1) go to 10
anoise=1.e-4
tau=anorm*anoise*10.
10 continue
c initialize the data generation function
c ..
dummy=gen(-1.)
write (6,230)
write (6,240) anoise,anorm,tau
c
do 180 mn1=1,6,5
mn2=mn1+2
do 180 m=mn1,mn2
do 180 n=mn1,mn2
write (6,250) m,n
c generate data
c ..
do 20 i=1,m
do 20 j=1,n
20 a(i,j)=gen(anoise)
do 30 i=1,m
30 b(i)=gen(anoise)
c
c ****** call hfti ******
c
call hfti(a,mda,m,n,b,1,1,tau,krank,srsmsq,h,g,ip)
c
c
write (6,260) krank
write (6,200) (i,b(i),i=1,n)
write (6,190) srsmsq
if (krank.lt.n) go to 180
c ****** algorithm cov bigins here ******
c ..
do 40 j=1,n
40 a(j,j)=1./a(j,j)
if (n.eq.1) go to 70
nm1=n-1
do 60 i=1,nm1
ip1=i+1
do 60 j=ip1,n
jm1=j-1
sm=0.d0
do 50 l=i,jm1
50 sm=sm+a(i,l)*dble(a(l,j))
60 a(i,j)=-sm*a(j,j)
c ..
c the upper triangle of a has been inverted
c upon itself.
70 do 90 i=1,n
do 90 j=i,n
sm=0.d0
do 80 l=j,n
80 sm=sm+a(i,l)*dble(a(j,l))
90 a(i,j)=sm
if (n.lt.2) go to 160
do 150 ii=2,n
i=n+1-ii
if (ip(i).eq.i) go to 150
k=ip(i)
tmp=a(i,i)
a(i,i)=a(k,k)
a(k,k)=tmp
if (i.eq.1) go to 110
do 100 l=2,i
tmp=a(l-1,i)
a(l-1,i)=a(l-1,k)
100 a(l-1,k)=tmp
110 ip1=i+1
km1=k-1
if (ip1.gt.km1) go to 130
do 120 l=ip1,km1
tmp=a(i,l)
a(i,l)=a(l,k)
120 a(l,k)=tmp
130 if (k.eq.n) go to 150
kp1=k+1
do 140 l=kp1,n
tmp=a(i,l)
a(i,l)=a(k,l)
140 a(k,l)=tmp
150 continue
160 continue
c ..
c covariance has been computed and repermuted.
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,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 'hfti' //(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 prog2. this program demonstates the algorithms,16
1h hfti and cov.)
240 format (1h0,54hthe relative noise level of the generated data will
1 be,e16.4/33h0the matrix norm is approximately,e12.4/43h0the absol
2ute pseudorank tolerance, tau, is,e12.4)
250 format (1h0////9h0 m n/1x,2i4)
260 format (1h0,8x,12hpseudorank =,i4)
end
|
