home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tred1.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  136 lines
  1.       subroutine tred1(nm,n,a,d,e,e2)
  2. c
  3.       integer i,j,k,l,n,ii,nm,jp1
  4.       double precision a(nm,n),d(n),e(n),e2(n)
  5.       double precision f,g,h,scale
  6. c
  7. c     this subroutine is a translation of the algol procedure tred1,
  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
  12. c     to a symmetric tridiagonal matrix using
  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        a contains information about the orthogonal trans-
  29. c          formations used in the reduction in its strict lower
  30. c          triangle.  the full upper triangle of a is unaltered.
  31. c
  32. c        d contains the diagonal elements of the tridiagonal matrix.
  33. c
  34. c        e contains the subdiagonal elements of the tridiagonal
  35. c          matrix in its last n-1 positions.  e(1) is set to zero.
  36. c
  37. c        e2 contains the squares of the corresponding elements of e.
  38. c          e2 may coincide with e if the squares are not needed.
  39. c
  40. c     questions and comments should be directed to burton s. garbow,
  41. c     mathematics and computer science div, argonne national laboratory
  42. c
  43. c     this version dated august 1983.
  44. c
  45. c     ------------------------------------------------------------------
  46. c
  47.       do 100 i = 1, n
  48.          d(i) = a(n,i)
  49.          a(n,i) = a(i,i)
  50.   100 continue
  51. c     .......... for i=n step -1 until 1 do -- ..........
  52.       do 300 ii = 1, n
  53.          i = n + 1 - ii
  54.          l = i - 1
  55.          h = 0.0d0
  56.          scale = 0.0d0
  57.          if (l .lt. 1) go to 130
  58. c     .......... scale row (algol tol then not needed) ..........
  59.          do 120 k = 1, l
  60.   120    scale = scale + dabs(d(k))
  61. c
  62.          if (scale .ne. 0.0d0) go to 140
  63. c
  64.          do 125 j = 1, l
  65.             d(j) = a(l,j)
  66.             a(l,j) = a(i,j)
  67.             a(i,j) = 0.0d0
  68.   125    continue
  69. c
  70.   130    e(i) = 0.0d0
  71.          e2(i) = 0.0d0
  72.          go to 300
  73. c
  74.   140    do 150 k = 1, l
  75.             d(k) = d(k) / scale
  76.             h = h + d(k) * d(k)
  77.   150    continue
  78. c
  79.          e2(i) = scale * scale * h
  80.          f = d(l)
  81.          g = -dsign(dsqrt(h),f)
  82.          e(i) = scale * g
  83.          h = h - f * g
  84.          d(l) = f - g
  85.          if (l .eq. 1) go to 285
  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.             g = e(j) + a(j,j) * f
  93.             jp1 = j + 1
  94.             if (l .lt. jp1) go to 220
  95. c
  96.             do 200 k = jp1, l
  97.                g = g + a(k,j) * d(k)
  98.                e(k) = e(k) + a(k,j) * f
  99.   200       continue
  100. c
  101.   220       e(j) = g
  102.   240    continue
  103. c     .......... form p ..........
  104.          f = 0.0d0
  105. c
  106.          do 245 j = 1, l
  107.             e(j) = e(j) / h
  108.             f = f + e(j) * d(j)
  109.   245    continue
  110. c
  111.          h = f / (h + h)
  112. c     .......... form q ..........
  113.          do 250 j = 1, l
  114.   250    e(j) = e(j) - h * d(j)
  115. c     .......... form reduced a ..........
  116.          do 280 j = 1, l
  117.             f = d(j)
  118.             g = e(j)
  119. c
  120.             do 260 k = j, l
  121.   260       a(k,j) = a(k,j) - f * e(k) - g * d(k)
  122. c
  123.   280    continue
  124. c
  125.   285    do 290 j = 1, l
  126.             f = d(j)
  127.             d(j) = a(l,j)
  128.             a(l,j) = a(i,j)
  129.             a(i,j) = f * scale
  130.   290    continue
  131. c
  132.   300 continue
  133. c
  134.       return
  135.       end
  136.