home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / htridi.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  155 lines
  1.       subroutine htridi(nm,n,ar,ai,d,e,e2,tau)
  2. c
  3.       integer i,j,k,l,n,ii,nm,jp1
  4.       double precision ar(nm,n),ai(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 tred1, 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
  13. c     to a real symmetric tridiagonal matrix using
  14. c     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        ar and ai contain the real and imaginary parts,
  25. c          respectively, of the complex hermitian input matrix.
  26. c          only the lower triangle of the matrix need be supplied.
  27. c
  28. c     on output
  29. c
  30. c        ar and ai contain information about the unitary trans-
  31. c          formations used in the reduction in their full lower
  32. c          triangles.  their strict upper triangles and the
  33. c          diagonal of ar are unaltered.
  34. c
  35. c        d contains the diagonal elements of the the tridiagonal matrix.
  36. c
  37. c        e contains the subdiagonal elements of the tridiagonal
  38. c          matrix in its last n-1 positions.  e(1) is set to zero.
  39. c
  40. c        e2 contains the squares of the corresponding elements of e.
  41. c          e2 may coincide with e if the squares are not needed.
  42. c
  43. c        tau contains further information about the transformations.
  44. c
  45. c     calls pythag for  dsqrt(a*a + b*b) .
  46. c
  47. c     questions and comments should be directed to burton s. garbow,
  48. c     mathematics and computer science div, argonne national laboratory
  49. c
  50. c     this version dated august 1983.
  51. c
  52. c     ------------------------------------------------------------------
  53. c
  54.       tau(1,n) = 1.0d0
  55.       tau(2,n) = 0.0d0
  56. c
  57.       do 100 i = 1, n
  58.   100 d(i) = ar(i,i)
  59. c     .......... for i=n step -1 until 1 do -- ..........
  60.       do 300 ii = 1, n
  61.          i = n + 1 - ii
  62.          l = i - 1
  63.          h = 0.0d0
  64.          scale = 0.0d0
  65.          if (l .lt. 1) go to 130
  66. c     .......... scale row (algol tol then not needed) ..........
  67.          do 120 k = 1, l
  68.   120    scale = scale + dabs(ar(i,k)) + dabs(ai(i,k))
  69. c
  70.          if (scale .ne. 0.0d0) go to 140
  71.          tau(1,l) = 1.0d0
  72.          tau(2,l) = 0.0d0
  73.   130    e(i) = 0.0d0
  74.          e2(i) = 0.0d0
  75.          go to 290
  76. c
  77.   140    do 150 k = 1, l
  78.             ar(i,k) = ar(i,k) / scale
  79.             ai(i,k) = ai(i,k) / scale
  80.             h = h + ar(i,k) * ar(i,k) + ai(i,k) * ai(i,k)
  81.   150    continue
  82. c
  83.          e2(i) = scale * scale * h
  84.          g = dsqrt(h)
  85.          e(i) = scale * g
  86.          f = pythag(ar(i,l),ai(i,l))
  87. c     .......... form next diagonal element of matrix t ..........
  88.          if (f .eq. 0.0d0) go to 160
  89.          tau(1,l) = (ai(i,l) * tau(2,i) - ar(i,l) * tau(1,i)) / f
  90.          si = (ar(i,l) * tau(2,i) + ai(i,l) * tau(1,i)) / f
  91.          h = h + f * g
  92.          g = 1.0d0 + g / f
  93.          ar(i,l) = g * ar(i,l)
  94.          ai(i,l) = g * ai(i,l)
  95.          if (l .eq. 1) go to 270
  96.          go to 170
  97.   160    tau(1,l) = -tau(1,i)
  98.          si = tau(2,i)
  99.          ar(i,l) = g
  100.   170    f = 0.0d0
  101. c
  102.          do 240 j = 1, l
  103.             g = 0.0d0
  104.             gi = 0.0d0
  105. c     .......... form element of a*u ..........
  106.             do 180 k = 1, j
  107.                g = g + ar(j,k) * ar(i,k) + ai(j,k) * ai(i,k)
  108.                gi = gi - ar(j,k) * ai(i,k) + ai(j,k) * ar(i,k)
  109.   180       continue
  110. c
  111.             jp1 = j + 1
  112.             if (l .lt. jp1) go to 220
  113. c
  114.             do 200 k = jp1, l
  115.                g = g + ar(k,j) * ar(i,k) - ai(k,j) * ai(i,k)
  116.                gi = gi - ar(k,j) * ai(i,k) - ai(k,j) * ar(i,k)
  117.   200       continue
  118. c     .......... form element of p ..........
  119.   220       e(j) = g / h
  120.             tau(2,j) = gi / h
  121.             f = f + e(j) * ar(i,j) - tau(2,j) * ai(i,j)
  122.   240    continue
  123. c
  124.          hh = f / (h + h)
  125. c     .......... form reduced a ..........
  126.          do 260 j = 1, l
  127.             f = ar(i,j)
  128.             g = e(j) - hh * f
  129.             e(j) = g
  130.             fi = -ai(i,j)
  131.             gi = tau(2,j) - hh * fi
  132.             tau(2,j) = -gi
  133. c
  134.             do 260 k = 1, j
  135.                ar(j,k) = ar(j,k) - f * e(k) - g * ar(i,k)
  136.      x                           + fi * tau(2,k) + gi * ai(i,k)
  137.                ai(j,k) = ai(j,k) - f * tau(2,k) - g * ai(i,k)
  138.      x                           - fi * e(k) - gi * ar(i,k)
  139.   260    continue
  140. c
  141.   270    do 280 k = 1, l
  142.             ar(i,k) = scale * ar(i,k)
  143.             ai(i,k) = scale * ai(i,k)
  144.   280    continue
  145. c
  146.          tau(2,l) = -si
  147.   290    hh = d(i)
  148.          d(i) = ar(i,i)
  149.          ar(i,i) = hh
  150.          ai(i,i) = scale * dsqrt(h)
  151.   300 continue
  152. c
  153.       return
  154.       end
  155.