home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tred3.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  128 lines
  1.       subroutine tred3(n,nv,a,d,e,e2)
  2. c
  3.       integer i,j,k,l,n,ii,iz,jk,nv,jm1
  4.       double precision a(nv),d(n),e(n),e2(n)
  5.       double precision f,g,h,hh,scale
  6. c
  7. c     this subroutine is a translation of the algol procedure tred3,
  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, stored as
  12. c     a one-dimensional array, to a symmetric tridiagonal matrix
  13. c     using orthogonal similarity transformations.
  14. c
  15. c     on input
  16. c
  17. c        n is the order of the matrix.
  18. c
  19. c        nv must be set to the dimension of the array parameter a
  20. c          as declared in the calling program dimension statement.
  21. c
  22. c        a contains the lower triangle of the real symmetric
  23. c          input matrix, stored row-wise as a one-dimensional
  24. c          array, in its first n*(n+1)/2 positions.
  25. c
  26. c     on output
  27. c
  28. c        a contains information about the orthogonal
  29. c          transformations used in the reduction.
  30. c
  31. c        d contains the diagonal elements of the tridiagonal matrix.
  32. c
  33. c        e contains the subdiagonal elements of the tridiagonal
  34. c          matrix in its last n-1 positions.  e(1) is set to zero.
  35. c
  36. c        e2 contains the squares of the corresponding elements of e.
  37. c          e2 may coincide with e if the squares are not needed.
  38. c
  39. c     questions and comments should be directed to burton s. garbow,
  40. c     mathematics and computer science div, argonne national laboratory
  41. c
  42. c     this version dated august 1983.
  43. c
  44. c     ------------------------------------------------------------------
  45. c
  46. c     .......... for i=n step -1 until 1 do -- ..........
  47.       do 300 ii = 1, n
  48.          i = n + 1 - ii
  49.          l = i - 1
  50.          iz = (i * l) / 2
  51.          h = 0.0d0
  52.          scale = 0.0d0
  53.          if (l .lt. 1) go to 130
  54. c     .......... scale row (algol tol then not needed) ..........
  55.          do 120 k = 1, l
  56.             iz = iz + 1
  57.             d(k) = a(iz)
  58.             scale = scale + dabs(d(k))
  59.   120    continue
  60. c
  61.          if (scale .ne. 0.0d0) go to 140
  62.   130    e(i) = 0.0d0
  63.          e2(i) = 0.0d0
  64.          go to 290
  65. c
  66.   140    do 150 k = 1, l
  67.             d(k) = d(k) / scale
  68.             h = h + d(k) * d(k)
  69.   150    continue
  70. c
  71.          e2(i) = scale * scale * h
  72.          f = d(l)
  73.          g = -dsign(dsqrt(h),f)
  74.          e(i) = scale * g
  75.          h = h - f * g
  76.          d(l) = f - g
  77.          a(iz) = scale * d(l)
  78.          if (l .eq. 1) go to 290
  79.          jk = 1
  80. c
  81.          do 240 j = 1, l
  82.             f = d(j)
  83.             g = 0.0d0
  84.             jm1 = j - 1
  85.             if (jm1 .lt. 1) go to 220
  86. c
  87.             do 200 k = 1, jm1
  88.                g = g + a(jk) * d(k)
  89.                e(k) = e(k) + a(jk) * f
  90.                jk = jk + 1
  91.   200       continue
  92. c     
  93.   220       e(j) = g + a(jk) * f
  94.             jk = jk + 1
  95.   240    continue
  96. c     .......... form p ..........
  97.          f = 0.0d0
  98. c
  99.          do 245 j = 1, l
  100.             e(j) = e(j) / h
  101.             f = f + e(j) * d(j)
  102.   245    continue
  103. c
  104.          hh = f / (h + h)
  105. c     .......... form q ..........
  106.          do 250 j = 1, l
  107.   250    e(j) = e(j) - hh * d(j)
  108. c
  109.          jk = 1
  110. c     .......... form reduced a ..........
  111.          do 280 j = 1, l
  112.             f = d(j)
  113.             g = e(j)
  114. c
  115.             do 260 k = 1, j
  116.                a(jk) = a(jk) - f * e(k) - g * d(k)
  117.                jk = jk + 1
  118.   260       continue
  119. c
  120.   280    continue
  121. c
  122.   290    d(i) = a(iz+1)
  123.          a(iz+1) = scale * dsqrt(h)
  124.   300 continue
  125. c
  126.       return
  127.       end
  128.