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
|
c prog3
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 the use of subroutine svdrs to compute the
c singular value decomposition of a matrix, a, and solve a least
c squares problem, a*x=b.
c
c the s.v.d. a= u*(s,0)**t*v**t is
c computed so that..
c (1) u**t*b replaces the m+1 st. col. of b.
c
c (2) u**t replaces the m by m identity in
c the first m cols. of b.
c
c (3) v replaces the first n rows and cols.
c of a.
c
c (4) the diagonal entries of the s matrix
c replace the first n entries of the array s.
c
c the array s( ) must be dimensioned at least 3*n .
c
dimension a(8,8),b(8,9),s(24),x(8),aa(8,8)
real gen,anoise
double precision sm
data mda,mdb/8,8/
c
do 150 noise=1,2
anoise = 0.
rho = 1.e-3
if(noise .eq. 1) go to 5
anoise = 1.e-4
rho = 10. * anoise
5 continue
write(6,230)
write(6,240) anoise,rho
c initialize data generation function
c ..
dummy= gen(-1.)
c
do 150 mn1=1,6,5
mn2=mn1+2
do 150 m=mn1,mn2
do 150 n=mn1,mn2
write (6,160) m,n
do 20 i=1,m
do 10 j=1,m
10 b(i,j)=0.
b(i,i)=1.
do 20 j=1,n
a(i,j)= gen(anoise)
20 aa(i,j)=a(i,j)
do 30 i=1,m
30 b(i,m+1)= gen(anoise)
c
c the arrays are now filled in..
c compute the s.v.d.
c ******************************************************
call svdrs (a,mda,m,n,b,mdb,m+1,s)
c ******************************************************
c
write (6,170)
write (6,220) (i,s(i),i=1,n)
write (6,180)
write (6,220) (i,b(i,m+1),i=1,m)
c
c test for disparity of ratio of singular values.
c ..
krank=n
tau=rho*s(1)
do 40 i=1,n
if (s(i).le.tau) go to 50
40 continue
go to 55
50 krank=i-1
55 write(6,250) tau, krank
c compute solution vector assuming pseudorank is krank
c ..
60 do 70 i=1,krank
70 b(i,m+1)=b(i,m+1)/s(i)
do 90 i=1,n
sm=0.d0
do 80 j=1,krank
80 sm=sm+a(i,j)*dble(b(j,m+1))
90 x(i)=sm
c compute predicted norm of residual vector.
c ..
srsmsq=0.
if (krank.eq.m) go to 110
kp1=krank+1
do 100 i=kp1,m
100 srsmsq=srsmsq+b(i,m+1)**2
srsmsq=sqrt(srsmsq)
c
110 continue
write (6,190)
write (6,220) (i,x(i),i=1,n)
write (6,200) srsmsq
c compute the frobenius norm of a**t- v*(s,0)*u**t.
c
c compute v*s first.
c
minmn=min0(m,n)
do 120 j=1,minmn
do 120 i=1,n
120 a(i,j)=a(i,j)*s(j)
dn=0.
do 140 j=1,m
do 140 i=1,n
sm=0.d0
do 130 l=1,minmn
130 sm=sm+a(i,l)*dble(b(l,j))
c computed difference of (i,j) th entry
c of a**t-v*(s,0)*u**t.
c ..
t=aa(j,i)-sm
140 dn=dn+t**2
dn=sqrt(dn)/(sqrt(float(n))*s(1))
write (6,210) dn
150 continue
stop
160 format (1h0////9h0 m n/1x,2i4)
170 format (1h0,8x,25hsingular values of matrix)
180 format (1h0,8x,30htransformed right side, u**t*b)
190 format (1h0,8x,33hestimated parameters, x=a**(+)*b,21h computed b
1y 'svdrs' )
200 format (1h0,8x,24hresidual vector length =,e12.4)
210 format (1h0,8x,43hfrobenius norm (a-u*(s,0)**t*v**t)/(sqrt(n),
*22h*spectral norm of a) =,e12.4)
220 format (9x,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8,i5,e16.8)
230 format(51h1prog3. this program demonstrates the algorithm,
* 9h, svdrs. )
240 format(55h0the relative noise level of the generated data will be,
*e16.4/44h0the relative tolerance, rho, for pseudorank,
*17h determination is,e16.4)
250 format(1h0,8x,36habsolute pseudorank tolerance, tau =,
*e12.4,10x,12hpseudorank =,i5)
end
|
