home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tred2.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  165 lines
  1.       subroutine tred2(nm,n,a,d,e,z)
  2. c
  3.       integer i,j,k,l,n,ii,nm,jp1
  4.       double precision a(nm,n),d(n),e(n),z(nm,n)
  5.       double precision f,g,h,hh,scale
  6. c
  7. c     this subroutine is a translation of the algol procedure tred2,
  8. c     num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
  9. c     handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
  10. c
  11. c     this subroutine reduces a real symmetric matrix to a
  12. c     symmetric tridiagonal matrix using and accumulating
  13. c     orthogonal similarity transformations.
  14. c
  15. c     on input
  16. c
  17. c        nm must be set to the row dimension of two-dimensional
  18. c          array parameters as declared in the calling program
  19. c          dimension statement.
  20. c
  21. c        n is the order of the matrix.
  22. c
  23. c        a contains the real symmetric input matrix.  only the
  24. c          lower triangle of the matrix need be supplied.
  25. c
  26. c     on output
  27. c
  28. c        d contains the diagonal elements of the tridiagonal matrix.
  29. c
  30. c        e contains the subdiagonal elements of the tridiagonal
  31. c          matrix in its last n-1 positions.  e(1) is set to zero.
  32. c
  33. c        z contains the orthogonal transformation matrix
  34. c          produced in the reduction.
  35. c
  36. c        a and z may coincide.  if distinct, a is unaltered.
  37. c
  38. c     questions and comments should be directed to burton s. garbow,
  39. c     mathematics and computer science div, argonne national laboratory
  40. c
  41. c     this version dated august 1983.
  42. c
  43. c     ------------------------------------------------------------------
  44. c
  45.       do 100 i = 1, n
  46. c
  47.          do 80 j = i, n
  48.    80    z(j,i) = a(j,i)
  49. c
  50.          d(i) = a(n,i)
  51.   100 continue
  52. c
  53.       if (n .eq. 1) go to 510
  54. c     .......... for i=n step -1 until 2 do -- ..........
  55.       do 300 ii = 2, n
  56.          i = n + 2 - ii
  57.          l = i - 1
  58.          h = 0.0d0
  59.          scale = 0.0d0
  60.          if (l .lt. 2) go to 130
  61. c     .......... scale row (algol tol then not needed) ..........
  62.          do 120 k = 1, l
  63.   120    scale = scale + dabs(d(k))
  64. c
  65.          if (scale .ne. 0.0d0) go to 140
  66.   130    e(i) = d(l)
  67. c
  68.          do 135 j = 1, l
  69.             d(j) = z(l,j)
  70.             z(i,j) = 0.0d0
  71.             z(j,i) = 0.0d0
  72.   135    continue
  73. c
  74.          go to 290
  75. c
  76.   140    do 150 k = 1, l
  77.             d(k) = d(k) / scale
  78.             h = h + d(k) * d(k)
  79.   150    continue
  80. c
  81.          f = d(l)
  82.          g = -dsign(dsqrt(h),f)
  83.          e(i) = scale * g
  84.          h = h - f * g
  85.          d(l) = f - g
  86. c     .......... form a*u ..........
  87.          do 170 j = 1, l
  88.   170    e(j) = 0.0d0
  89. c
  90.          do 240 j = 1, l
  91.             f = d(j)
  92.             z(j,i) = f
  93.             g = e(j) + z(j,j) * f
  94.             jp1 = j + 1
  95.             if (l .lt. jp1) go to 220
  96. c
  97.             do 200 k = jp1, l
  98.                g = g + z(k,j) * d(k)
  99.                e(k) = e(k) + z(k,j) * f
  100.   200       continue
  101. c
  102.   220       e(j) = g
  103.   240    continue
  104. c     .......... form p ..........
  105.          f = 0.0d0
  106. c
  107.          do 245 j = 1, l
  108.             e(j) = e(j) / h
  109.             f = f + e(j) * d(j)
  110.   245    continue
  111. c
  112.          hh = f / (h + h)
  113. c     .......... form q ..........
  114.          do 250 j = 1, l
  115.   250    e(j) = e(j) - hh * d(j)
  116. c     .......... form reduced a ..........
  117.          do 280 j = 1, l
  118.             f = d(j)
  119.             g = e(j)
  120. c
  121.             do 260 k = j, l
  122.   260       z(k,j) = z(k,j) - f * e(k) - g * d(k)
  123. c
  124.             d(j) = z(l,j)
  125.             z(i,j) = 0.0d0
  126.   280    continue
  127. c
  128.   290    d(i) = h
  129.   300 continue
  130. c     .......... accumulation of transformation matrices ..........
  131.       do 500 i = 2, n
  132.          l = i - 1
  133.          z(n,l) = z(l,l)
  134.          z(l,l) = 1.0d0
  135.          h = d(i)
  136.          if (h .eq. 0.0d0) go to 380
  137. c
  138.          do 330 k = 1, l
  139.   330    d(k) = z(k,i) / h
  140. c
  141.          do 360 j = 1, l
  142.             g = 0.0d0
  143. c
  144.             do 340 k = 1, l
  145.   340       g = g + z(k,i) * z(k,j)
  146. c
  147.             do 360 k = 1, l
  148.                z(k,j) = z(k,j) - g * d(k)
  149.   360    continue
  150. c
  151.   380    do 400 k = 1, l
  152.   400    z(k,i) = 0.0d0
  153. c
  154.   500 continue
  155. c
  156.   510 do 520 i = 1, n
  157.          d(i) = z(n,i)
  158.          z(n,i) = 0.0d0
  159.   520 continue
  160. c
  161.       z(n,n) = 1.0d0
  162.       e(1) = 0.0d0
  163.       return
  164.       end
  165.