home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / elmhes.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  99 lines
  1.       subroutine elmhes(nm,n,low,igh,a,int)
  2. c
  3.       integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1
  4.       double precision a(nm,n)
  5.       double precision x,y
  6.       integer int(igh)
  7. c
  8. c     this subroutine is a translation of the algol procedure elmhes,
  9. c     num. math. 12, 349-368(1968) by martin and wilkinson.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
  11. c
  12. c     given a real general matrix, this subroutine
  13. c     reduces a submatrix situated in rows and columns
  14. c     low through igh to upper hessenberg form by
  15. c     stabilized elementary similarity transformations.
  16. c
  17. c     on input
  18. c
  19. c        nm must be set to the row dimension of two-dimensional
  20. c          array parameters as declared in the calling program
  21. c          dimension statement.
  22. c
  23. c        n is the order of the matrix.
  24. c
  25. c        low and igh are integers determined by the balancing
  26. c          subroutine  balanc.  if  balanc  has not been used,
  27. c          set low=1, igh=n.
  28. c
  29. c        a contains the input matrix.
  30. c
  31. c     on output
  32. c
  33. c        a contains the hessenberg matrix.  the multipliers
  34. c          which were used in the reduction are stored in the
  35. c          remaining triangle under the hessenberg matrix.
  36. c
  37. c        int contains information on the rows and columns
  38. c          interchanged in the reduction.
  39. c          only elements low through igh are used.
  40. c
  41. c     questions and comments should be directed to burton s. garbow,
  42. c     mathematics and computer science div, argonne national laboratory
  43. c
  44. c     this version dated august 1983.
  45. c
  46. c     ------------------------------------------------------------------
  47. c
  48.       la = igh - 1
  49.       kp1 = low + 1
  50.       if (la .lt. kp1) go to 200
  51. c
  52.       do 180 m = kp1, la
  53.          mm1 = m - 1
  54.          x = 0.0d0
  55.          i = m
  56. c
  57.          do 100 j = m, igh
  58.             if (dabs(a(j,mm1)) .le. dabs(x)) go to 100
  59.             x = a(j,mm1)
  60.             i = j
  61.   100    continue
  62. c
  63.          int(m) = i
  64.          if (i .eq. m) go to 130
  65. c     .......... interchange rows and columns of a ..........
  66.          do 110 j = mm1, n
  67.             y = a(i,j)
  68.             a(i,j) = a(m,j)
  69.             a(m,j) = y
  70.   110    continue
  71. c
  72.          do 120 j = 1, igh
  73.             y = a(j,i)
  74.             a(j,i) = a(j,m)
  75.             a(j,m) = y
  76.   120    continue
  77. c     .......... end interchange ..........
  78.   130    if (x .eq. 0.0d0) go to 180
  79.          mp1 = m + 1
  80. c
  81.          do 160 i = mp1, igh
  82.             y = a(i,mm1)
  83.             if (y .eq. 0.0d0) go to 160
  84.             y = y / x
  85.             a(i,mm1) = y
  86. c
  87.             do 140 j = m, n
  88.   140       a(i,j) = a(i,j) - y * a(m,j)
  89. c
  90.             do 150 j = 1, igh
  91.   150       a(j,m) = a(j,m) + y * a(j,i)
  92. c
  93.   160    continue
  94. c
  95.   180 continue
  96. c
  97.   200 return
  98.       end
  99.