|
|
subroutine tred2(nm,n,a,d,e,z)
|
|
|
c
|
|
|
integer i,j,k,l,n,ii,nm,jp1
|
|
|
real a(nm,n),d(n),e(n),z(nm,n)
|
|
|
real f,g,h,hh,scale
|
|
|
c
|
|
|
c this subroutine is a translation of the algol procedure tred2,
|
|
|
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 to a
|
|
|
c symmetric tridiagonal matrix using and accumulating
|
|
|
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 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 z contains the orthogonal transformation matrix
|
|
|
c produced in the reduction.
|
|
|
c
|
|
|
c a and z may coincide. if distinct, a is unaltered.
|
|
|
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
|
|
|
c
|
|
|
do 80 j = i, n
|
|
|
80 z(j,i) = a(j,i)
|
|
|
c
|
|
|
d(i) = a(n,i)
|
|
|
100 continue
|
|
|
c
|
|
|
if (n .eq. 1) go to 510
|
|
|
c .......... for i=n step -1 until 2 do -- ..........
|
|
|
do 300 ii = 2, n
|
|
|
i = n + 2 - ii
|
|
|
l = i - 1
|
|
|
h = 0.0e0
|
|
|
scale = 0.0e0
|
|
|
if (l .lt. 2) 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
|
|
|
130 e(i) = d(l)
|
|
|
c
|
|
|
do 135 j = 1, l
|
|
|
d(j) = z(l,j)
|
|
|
z(i,j) = 0.0e0
|
|
|
z(j,i) = 0.0e0
|
|
|
135 continue
|
|
|
c
|
|
|
go to 290
|
|
|
c
|
|
|
140 do 150 k = 1, l
|
|
|
d(k) = d(k) / scale
|
|
|
h = h + d(k) * d(k)
|
|
|
150 continue
|
|
|
c
|
|
|
f = d(l)
|
|
|
g = -sign(sqrt(h),f)
|
|
|
e(i) = scale * g
|
|
|
h = h - f * g
|
|
|
d(l) = f - g
|
|
|
c .......... form a*u ..........
|
|
|
do 170 j = 1, l
|
|
|
170 e(j) = 0.0e0
|
|
|
c
|
|
|
do 240 j = 1, l
|
|
|
f = d(j)
|
|
|
z(j,i) = f
|
|
|
g = e(j) + z(j,j) * f
|
|
|
jp1 = j + 1
|
|
|
if (l .lt. jp1) go to 220
|
|
|
c
|
|
|
do 200 k = jp1, l
|
|
|
g = g + z(k,j) * d(k)
|
|
|
e(k) = e(k) + z(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
|
|
|
hh = f / (h + h)
|
|
|
c .......... form q ..........
|
|
|
do 250 j = 1, l
|
|
|
250 e(j) = e(j) - hh * d(j)
|
|
|
c .......... form reduced a ..........
|
|
|
do 280 j = 1, l
|
|
|
f = d(j)
|
|
|
g = e(j)
|
|
|
c
|
|
|
do 260 k = j, l
|
|
|
260 z(k,j) = z(k,j) - f * e(k) - g * d(k)
|
|
|
c
|
|
|
d(j) = z(l,j)
|
|
|
z(i,j) = 0.0e0
|
|
|
280 continue
|
|
|
c
|
|
|
290 d(i) = h
|
|
|
300 continue
|
|
|
c .......... accumulation of transformation matrices ..........
|
|
|
do 500 i = 2, n
|
|
|
l = i - 1
|
|
|
z(n,l) = z(l,l)
|
|
|
z(l,l) = 1.0e0
|
|
|
h = d(i)
|
|
|
if (h .eq. 0.0e0) go to 380
|
|
|
c
|
|
|
do 330 k = 1, l
|
|
|
330 d(k) = z(k,i) / h
|
|
|
c
|
|
|
do 360 j = 1, l
|
|
|
g = 0.0e0
|
|
|
c
|
|
|
do 340 k = 1, l
|
|
|
340 g = g + z(k,i) * z(k,j)
|
|
|
c
|
|
|
do 360 k = 1, l
|
|
|
z(k,j) = z(k,j) - g * d(k)
|
|
|
360 continue
|
|
|
c
|
|
|
380 do 400 k = 1, l
|
|
|
400 z(k,i) = 0.0e0
|
|
|
c
|
|
|
500 continue
|
|
|
c
|
|
|
510 do 520 i = 1, n
|
|
|
d(i) = z(n,i)
|
|
|
z(n,i) = 0.0e0
|
|
|
520 continue
|
|
|
c
|
|
|
z(n,n) = 1.0e0
|
|
|
e(1) = 0.0e0
|
|
|
return
|
|
|
end
|
|
|
|
