home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / comhes.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  119 lines
  1.       subroutine comhes(nm,n,low,igh,ar,ai,int)
  2. c
  3.       integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1
  4.       double precision ar(nm,n),ai(nm,n)
  5.       double precision xr,xi,yr,yi
  6.       integer int(igh)
  7. c
  8. c     this subroutine is a translation of the algol procedure comhes,
  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 complex 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  cbal.  if  cbal  has not been used,
  27. c          set low=1, igh=n.
  28. c
  29. c        ar and ai contain the real and imaginary parts,
  30. c          respectively, of the complex input matrix.
  31. c
  32. c     on output
  33. c
  34. c        ar and ai contain the real and imaginary parts,
  35. c          respectively, of the hessenberg matrix.  the
  36. c          multipliers which were used in the reduction
  37. c          are stored in the remaining triangles under the
  38. c          hessenberg matrix.
  39. c
  40. c        int contains information on the rows and columns
  41. c          interchanged in the reduction.
  42. c          only elements low through igh are used.
  43. c
  44. c     calls cdiv for complex division.
  45. c
  46. c     questions and comments should be directed to burton s. garbow,
  47. c     mathematics and computer science div, argonne national laboratory
  48. c
  49. c     this version dated august 1983.
  50. c
  51. c     ------------------------------------------------------------------
  52. c
  53.       la = igh - 1
  54.       kp1 = low + 1
  55.       if (la .lt. kp1) go to 200
  56. c
  57.       do 180 m = kp1, la
  58.          mm1 = m - 1
  59.          xr = 0.0d0
  60.          xi = 0.0d0
  61.          i = m
  62. c
  63.          do 100 j = m, igh
  64.             if (dabs(ar(j,mm1)) + dabs(ai(j,mm1))
  65.      x         .le. dabs(xr) + dabs(xi)) go to 100
  66.             xr = ar(j,mm1)
  67.             xi = ai(j,mm1)
  68.             i = j
  69.   100    continue
  70. c
  71.          int(m) = i
  72.          if (i .eq. m) go to 130
  73. c     .......... interchange rows and columns of ar and ai ..........
  74.          do 110 j = mm1, n
  75.             yr = ar(i,j)
  76.             ar(i,j) = ar(m,j)
  77.             ar(m,j) = yr
  78.             yi = ai(i,j)
  79.             ai(i,j) = ai(m,j)
  80.             ai(m,j) = yi
  81.   110    continue
  82. c
  83.          do 120 j = 1, igh
  84.             yr = ar(j,i)
  85.             ar(j,i) = ar(j,m)
  86.             ar(j,m) = yr
  87.             yi = ai(j,i)
  88.             ai(j,i) = ai(j,m)
  89.             ai(j,m) = yi
  90.   120    continue
  91. c     .......... end interchange ..........
  92.   130    if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 180
  93.          mp1 = m + 1
  94. c
  95.          do 160 i = mp1, igh
  96.             yr = ar(i,mm1)
  97.             yi = ai(i,mm1)
  98.             if (yr .eq. 0.0d0 .and. yi .eq. 0.0d0) go to 160
  99.             call cdiv(yr,yi,xr,xi,yr,yi)
  100.             ar(i,mm1) = yr
  101.             ai(i,mm1) = yi
  102. c
  103.             do 140 j = m, n
  104.                ar(i,j) = ar(i,j) - yr * ar(m,j) + yi * ai(m,j)
  105.                ai(i,j) = ai(i,j) - yr * ai(m,j) - yi * ar(m,j)
  106.   140       continue
  107. c
  108.             do 150 j = 1, igh
  109.                ar(j,m) = ar(j,m) + yr * ar(j,i) - yi * ai(j,i)
  110.                ai(j,m) = ai(j,m) + yr * ai(j,i) + yi * ar(j,i)
  111.   150       continue
  112. c
  113.   160    continue
  114. c
  115.   180 continue
  116. c
  117.   200 return
  118.       end
  119.