tred2.f
164 lines
| 4.0 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | 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 | ||||
