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