home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / balgen / reduce.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  182 lines
  1.       subroutine reduce (n,ma,a,mb,b,low,igh,cscale,wk)
  2. c
  3. c     *****parameters:
  4.       integer igh,low,ma,mb,n
  5.       double precision a(ma,n),b(mb,n),cscale(n),wk(n)
  6. c
  7. c     *****local variables:
  8.       integer i,iflow,ii,ip1,is,j,jp1,k,l,lm1,m
  9.       double precision f
  10. c
  11. c     *****functions:
  12. c     none
  13. c
  14. c     *****subroutines called:
  15. c     none
  16. c
  17. c     ---------------------------------------------------------------
  18. c
  19. c     *****purpose:
  20. c     this subroutine reduces, if possible, the order of the
  21. c     generalized eigenvalue problem a*x = (lambda)*b*x by permuting
  22. c     the rows and columns of a and b so that they each have the
  23. c     form
  24. c                       u  x  y
  25. c                       0  c  z
  26. c                       0  0  r
  27. c
  28. c     where u and r are upper triangular and c, x, y, and z are
  29. c     arbitrary.  thus, the isolated eigenvalues corresponding to
  30. c     the triangular matrices are obtained by a division, leaving
  31. c     only eigenvalues corresponding to the center matrices to be
  32. c     computed.
  33. c     ref.:  ward, r. c., balancing the generalized eigenvalue
  34. c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
  35. c     141-152.
  36. c
  37. c     *****parameter description:
  38. c
  39. c     on input:
  40. c
  41. c       ma,mb   integer
  42. c               row dimensions of the arrays containing matrices
  43. c               a and b respectively, as declared in the main calling
  44. c               program dimension statement;
  45. c
  46. c       n       integer
  47. c               order of the matrices a and b;
  48. c
  49. c       a       real(ma,n)
  50. c               contains the a matrix of the generalized eigenproblem
  51. c               defined above;
  52. c
  53. c       b       real(mb,n)
  54. c               contains the b matrix of the generalized eigenproblem
  55. c               defined above.
  56. c
  57. c     on output:
  58. c
  59. c       a,b     contain the permuted a and b matrices;
  60. c
  61. c       low     integer
  62. c               beginning -1 of the submatrices of a and b
  63. c               containing the non-isolated eigenvalues;
  64. c
  65. c       igh     integer
  66. c               ending -1 of the submatrices of a and b
  67. c               containing the non-isolated eigenvalues.  if
  68. c               igh = 1 (low = 1 also), the permuted a and b
  69. c               matrices are upper triangular;
  70. c
  71. c       cscale  real(n)
  72. c               contains the required column permutations in its
  73. c               first low-1 and its igh+1 through n locations;
  74. c
  75. c       wk      real(n)
  76. c               contains the required row permutations in its first
  77. c               low-1 and its igh+1 through n locations.
  78. c
  79. c     *****algorithm notes:
  80. c     none
  81. c
  82. c     *****history:
  83. c     written by r. c. ward.......
  84. c
  85. c     ---------------------------------------------------------------
  86. c
  87.       k = 1
  88.       l = n
  89.       go to 20
  90. c
  91. c     find row with one nonzero in columns 1 through l
  92. c
  93.    10 continue
  94.       l = lm1
  95.       if (l .ne. 1) go to 20
  96.       wk(1) = 1
  97.       cscale(1) = 1
  98.       go to 200
  99.    20 continue
  100.       lm1 = l-1
  101.       do 70 ii = 1,l
  102.          i = l+1-ii
  103.          do 30 j = 1,lm1
  104.             jp1 = j+1
  105.             if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 40
  106.    30    continue
  107.          j = l
  108.          go to 60
  109.    40    continue
  110.          do 50 j = jp1,l
  111.             if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 70
  112.    50    continue
  113.          j = jp1-1
  114.    60    continue
  115.          m = l
  116.          iflow = 1
  117.          go to 150
  118.    70 continue
  119.       go to 90
  120. c
  121. c     find column with one nonzero in rows k through n
  122. c
  123.    80 continue
  124.       k = k+1
  125.    90 continue
  126.       do 140 j = k,l
  127.          do 100 i = k,lm1
  128.             ip1 = i+1
  129.             if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 110
  130.   100    continue
  131.          i = l
  132.          go to 130
  133.   110    continue
  134.          do 120 i = ip1,l
  135.             if (a(i,j) .ne. 0.0d0 .or. b(i,j) .ne. 0.0d0) go to 140
  136.   120    continue
  137.          i = ip1-1
  138.   130    continue
  139.          m = k
  140.          iflow = 2
  141.          go to 150
  142.   140 continue
  143.       go to 200
  144. c
  145. c     permute rows m and i
  146. c
  147.   150 continue
  148.       wk(m) = i
  149.       if (i .eq. m) go to 170
  150.       do 160 is = k,n
  151.          f = a(i,is)
  152.          a(i,is) = a(m,is)
  153.          a(m,is) = f
  154.          f = b(i,is)
  155.          b(i,is) = b(m,is)
  156.          b(m,is) = f
  157.   160 continue
  158. c
  159. c     permute columns m and j
  160. c
  161.   170 continue
  162.       cscale(m) = j
  163.       if (j .eq. m) go to 190
  164.       do 180 is = 1,l
  165.          f = a(is,j)
  166.          a(is,j) = a(is,m)
  167.          a(is,m) = f
  168.          f = b(is,j)
  169.          b(is,j) = b(is,m)
  170.          b(is,m) = f
  171.   180 continue
  172.   190 continue
  173.       go to (10,80), iflow
  174.   200 continue
  175.       low = k
  176.       igh = l
  177.       return
  178. c
  179. c     last line of reduce
  180. c
  181.       end
  182.