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
|
subroutine tred1(nm,n,a,d,e,e2)
c
integer i,j,k,l,n,ii,nm,jp1
double precision a(nm,n),d(n),e(n),e2(n)
double precision f,g,h,scale
c
c this subroutine is a translation of the algol procedure tred1,
c num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
c
c this subroutine reduces a real symmetric matrix
c to a symmetric tridiagonal matrix using
c orthogonal similarity transformations.
c
c on input
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement.
c
c n is the order of the matrix.
c
c a contains the real symmetric input matrix. only the
c lower triangle of the matrix need be supplied.
c
c on output
c
c a contains information about the orthogonal trans-
c formations used in the reduction in its strict lower
c triangle. the full upper triangle of a is unaltered.
c
c d contains the diagonal elements of the tridiagonal matrix.
c
c e contains the subdiagonal elements of the tridiagonal
c matrix in its last n-1 positions. e(1) is set to zero.
c
c e2 contains the squares of the corresponding elements of e.
c e2 may coincide with e if the squares are not needed.
c
c questions and comments should be directed to burton s. garbow,
c mathematics and computer science div, argonne national laboratory
c
c this version dated august 1983.
c
c ------------------------------------------------------------------
c
do 100 i = 1, n
d(i) = a(n,i)
a(n,i) = a(i,i)
100 continue
c .......... for i=n step -1 until 1 do -- ..........
do 300 ii = 1, n
i = n + 1 - ii
l = i - 1
h = 0.0d0
scale = 0.0d0
if (l .lt. 1) go to 130
c .......... scale row (algol tol then not needed) ..........
do 120 k = 1, l
120 scale = scale + dabs(d(k))
c
if (scale .ne. 0.0d0) go to 140
c
do 125 j = 1, l
d(j) = a(l,j)
a(l,j) = a(i,j)
a(i,j) = 0.0d0
125 continue
c
130 e(i) = 0.0d0
e2(i) = 0.0d0
go to 300
c
140 do 150 k = 1, l
d(k) = d(k) / scale
h = h + d(k) * d(k)
150 continue
c
e2(i) = scale * scale * h
f = d(l)
g = -dsign(dsqrt(h),f)
e(i) = scale * g
h = h - f * g
d(l) = f - g
if (l .eq. 1) go to 285
c .......... form a*u ..........
do 170 j = 1, l
170 e(j) = 0.0d0
c
do 240 j = 1, l
f = d(j)
g = e(j) + a(j,j) * f
jp1 = j + 1
if (l .lt. jp1) go to 220
c
do 200 k = jp1, l
g = g + a(k,j) * d(k)
e(k) = e(k) + a(k,j) * f
200 continue
c
220 e(j) = g
240 continue
c .......... form p ..........
f = 0.0d0
c
do 245 j = 1, l
e(j) = e(j) / h
f = f + e(j) * d(j)
245 continue
c
h = f / (h + h)
c .......... form q ..........
do 250 j = 1, l
250 e(j) = e(j) - h * d(j)
c .......... form reduced a ..........
do 280 j = 1, l
f = d(j)
g = e(j)
c
do 260 k = j, l
260 a(k,j) = a(k,j) - f * e(k) - g * d(k)
c
280 continue
c
285 do 290 j = 1, l
f = d(j)
d(j) = a(l,j)
a(l,j) = a(i,j)
a(i,j) = f * scale
290 continue
c
300 continue
c
return
end
|
