home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / orthes.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  108 lines
  1.       subroutine orthes(nm,n,low,igh,a,ort)
  2. c
  3.       integer i,j,m,n,ii,jj,la,mp,nm,igh,kp1,low
  4.       double precision a(nm,n),ort(igh)
  5.       double precision f,g,h,scale
  6. c
  7. c     this subroutine is a translation of the algol procedure orthes,
  8. c     num. math. 12, 349-368(1968) by martin and wilkinson.
  9. c     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
  10. c
  11. c     given a real general matrix, this subroutine
  12. c     reduces a submatrix situated in rows and columns
  13. c     low through igh to upper hessenberg form by
  14. c     orthogonal 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        low and igh are integers determined by the balancing
  25. c          subroutine  balanc.  if  balanc  has not been used,
  26. c          set low=1, igh=n.
  27. c
  28. c        a contains the input matrix.
  29. c
  30. c     on output
  31. c
  32. c        a contains the hessenberg matrix.  information about
  33. c          the orthogonal transformations used in the reduction
  34. c          is stored in the remaining triangle under the
  35. c          hessenberg matrix.
  36. c
  37. c        ort contains further information about the transformations.
  38. c          only elements low through igh are used.
  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.       la = igh - 1
  48.       kp1 = low + 1
  49.       if (la .lt. kp1) go to 200
  50. c
  51.       do 180 m = kp1, la
  52.          h = 0.0d0
  53.          ort(m) = 0.0d0
  54.          scale = 0.0d0
  55. c     .......... scale column (algol tol then not needed) ..........
  56.          do 90 i = m, igh
  57.    90    scale = scale + dabs(a(i,m-1))
  58. c
  59.          if (scale .eq. 0.0d0) go to 180
  60.          mp = m + igh
  61. c     .......... for i=igh step -1 until m do -- ..........
  62.          do 100 ii = m, igh
  63.             i = mp - ii
  64.             ort(i) = a(i,m-1) / scale
  65.             h = h + ort(i) * ort(i)
  66.   100    continue
  67. c
  68.          g = -dsign(dsqrt(h),ort(m))
  69.          h = h - ort(m) * g
  70.          ort(m) = ort(m) - g
  71. c     .......... form (i-(u*ut)/h) * a ..........
  72.          do 130 j = m, n
  73.             f = 0.0d0
  74. c     .......... for i=igh step -1 until m do -- ..........
  75.             do 110 ii = m, igh
  76.                i = mp - ii
  77.                f = f + ort(i) * a(i,j)
  78.   110       continue
  79. c
  80.             f = f / h
  81. c
  82.             do 120 i = m, igh
  83.   120       a(i,j) = a(i,j) - f * ort(i)
  84. c
  85.   130    continue
  86. c     .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) ..........
  87.          do 160 i = 1, igh
  88.             f = 0.0d0
  89. c     .......... for j=igh step -1 until m do -- ..........
  90.             do 140 jj = m, igh
  91.                j = mp - jj
  92.                f = f + ort(j) * a(i,j)
  93.   140       continue
  94. c
  95.             f = f / h
  96. c
  97.             do 150 j = m, igh
  98.   150       a(i,j) = a(i,j) - f * ort(j)
  99. c
  100.   160    continue
  101. c
  102.          ort(m) = scale * ort(m)
  103.          a(m,m-1) = scale * g
  104.   180 continue
  105. c
  106.   200 return
  107.       end
  108.