home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / balanc.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  167 lines
  1.       subroutine balanc(nm,n,a,low,igh,scale)
  2. c
  3.       integer i,j,k,l,m,n,jj,nm,igh,low,iexc
  4.       double precision a(nm,n),scale(n)
  5.       double precision c,f,g,r,s,b2,radix
  6.       logical noconv
  7. c
  8. c     this subroutine is a translation of the algol procedure balance,
  9. c     num. math. 13, 293-304(1969) by parlett and reinsch.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971).
  11. c
  12. c     this subroutine balances a real matrix and isolates
  13. c     eigenvalues whenever possible.
  14. c
  15. c     on input
  16. c
  17. c        nm must be set to the row dimension of two-dimensional
  18. c          array parameters as declared in the calling program
  19. c          dimension statement.
  20. c
  21. c        n is the order of the matrix.
  22. c
  23. c        a contains the input matrix to be balanced.
  24. c
  25. c     on output
  26. c
  27. c        a contains the balanced matrix.
  28. c
  29. c        low and igh are two integers such that a(i,j)
  30. c          is equal to zero if
  31. c           (1) i is greater than j and
  32. c           (2) j=1,...,low-1 or i=igh+1,...,n.
  33. c
  34. c        scale contains information determining the
  35. c           permutations and scaling factors used.
  36. c
  37. c     suppose that the principal submatrix in rows low through igh
  38. c     has been balanced, that p(j) denotes the index interchanged
  39. c     with j during the permutation step, and that the elements
  40. c     of the diagonal matrix used are denoted by d(i,j).  then
  41. c        scale(j) = p(j),    for j = 1,...,low-1
  42. c                 = d(j,j),      j = low,...,igh
  43. c                 = p(j)         j = igh+1,...,n.
  44. c     the order in which the interchanges are made is n to igh+1,
  45. c     then 1 to low-1.
  46. c
  47. c     note that 1 is returned for igh if igh is zero formally.
  48. c
  49. c     the algol procedure exc contained in balance appears in
  50. c     balanc  in line.  (note that the algol roles of identifiers
  51. c     k,l have been reversed.)
  52. c
  53. c     questions and comments should be directed to burton s. garbow,
  54. c     mathematics and computer science div, argonne national laboratory
  55. c
  56. c     this version dated august 1983.
  57. c
  58. c     ------------------------------------------------------------------
  59. c
  60.       radix = 16.0d0
  61. c
  62.       b2 = radix * radix
  63.       k = 1
  64.       l = n
  65.       go to 100
  66. c     .......... in-line procedure for row and
  67. c                column exchange ..........
  68.    20 scale(m) = j
  69.       if (j .eq. m) go to 50
  70. c
  71.       do 30 i = 1, l
  72.          f = a(i,j)
  73.          a(i,j) = a(i,m)
  74.          a(i,m) = f
  75.    30 continue
  76. c
  77.       do 40 i = k, n
  78.          f = a(j,i)
  79.          a(j,i) = a(m,i)
  80.          a(m,i) = f
  81.    40 continue
  82. c
  83.    50 go to (80,130), iexc
  84. c     .......... search for rows isolating an eigenvalue
  85. c                and push them down ..........
  86.    80 if (l .eq. 1) go to 280
  87.       l = l - 1
  88. c     .......... for j=l step -1 until 1 do -- ..........
  89.   100 do 120 jj = 1, l
  90.          j = l + 1 - jj
  91. c
  92.          do 110 i = 1, l
  93.             if (i .eq. j) go to 110
  94.             if (a(j,i) .ne. 0.0d0) go to 120
  95.   110    continue
  96. c
  97.          m = l
  98.          iexc = 1
  99.          go to 20
  100.   120 continue
  101. c
  102.       go to 140
  103. c     .......... search for columns isolating an eigenvalue
  104. c                and push them left ..........
  105.   130 k = k + 1
  106. c
  107.   140 do 170 j = k, l
  108. c
  109.          do 150 i = k, l
  110.             if (i .eq. j) go to 150
  111.             if (a(i,j) .ne. 0.0d0) go to 170
  112.   150    continue
  113. c
  114.          m = k
  115.          iexc = 2
  116.          go to 20
  117.   170 continue
  118. c     .......... now balance the submatrix in rows k to l ..........
  119.       do 180 i = k, l
  120.   180 scale(i) = 1.0d0
  121. c     .......... iterative loop for norm reduction ..........
  122.   190 noconv = .false.
  123. c
  124.       do 270 i = k, l
  125.          c = 0.0d0
  126.          r = 0.0d0
  127. c
  128.          do 200 j = k, l
  129.             if (j .eq. i) go to 200
  130.             c = c + dabs(a(j,i))
  131.             r = r + dabs(a(i,j))
  132.   200    continue
  133. c     .......... guard against zero c or r due to underflow ..........
  134.          if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270
  135.          g = r / radix
  136.          f = 1.0d0
  137.          s = c + r
  138.   210    if (c .ge. g) go to 220
  139.          f = f * radix
  140.          c = c * b2
  141.          go to 210
  142.   220    g = r * radix
  143.   230    if (c .lt. g) go to 240
  144.          f = f / radix
  145.          c = c / b2
  146.          go to 230
  147. c     .......... now balance ..........
  148.   240    if ((c + r) / f .ge. 0.95d0 * s) go to 270
  149.          g = 1.0d0 / f
  150.          scale(i) = scale(i) * f
  151.          noconv = .true.
  152. c
  153.          do 250 j = k, n
  154.   250    a(i,j) = a(i,j) * g
  155. c
  156.          do 260 j = 1, l
  157.   260    a(j,i) = a(j,i) * f
  158. c
  159.   270 continue
  160. c
  161.       if (noconv) go to 190
  162. c
  163.   280 low = k
  164.       igh = l
  165.       return
  166.       end
  167.