|
|
subroutine tred1(nm,n,a,d,e,e2)
|
|
|
c
|
|
|
integer i,j,k,l,n,ii,nm,jp1
|
|
|
real a(nm,n),d(n),e(n),e2(n)
|
|
|
real 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.0e0
|
|
|
scale = 0.0e0
|
|
|
if (l .lt. 1) go to 130
|
|
|
c .......... scale row (algol tol then not needed) ..........
|
|
|
do 120 k = 1, l
|
|
|
120 scale = scale + abs(d(k))
|
|
|
c
|
|
|
if (scale .ne. 0.0e0) go to 140
|
|
|
c
|
|
|
do 125 j = 1, l
|
|
|
d(j) = a(l,j)
|
|
|
a(l,j) = a(i,j)
|
|
|
a(i,j) = 0.0e0
|
|
|
125 continue
|
|
|
c
|
|
|
130 e(i) = 0.0e0
|
|
|
e2(i) = 0.0e0
|
|
|
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 = -sign(sqrt(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.0e0
|
|
|
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.0e0
|
|
|
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
|
|
|
|
