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 / gradeq.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  178 lines
  1.       subroutine gradeq (n,ma,a,mb,b,low,igh,cperm,wk)
  2. c
  3. c     *****parameters:
  4.       integer igh,low,ma,mb,n
  5.       double precision a(ma,n),b(mb,n),cperm(n),wk(n,2)
  6. c
  7. c     *****local variables:
  8.       integer i,ighm1,im,ip1,j,jm,jp1,k
  9.       double precision cmax,rmax,suma,sumb,temp
  10. c
  11. c     *****fortran functions:
  12.       double precision dabs
  13. c
  14. c     *****subroutines called:
  15. c     none
  16. c
  17. c     ---------------------------------------------------------------
  18. c
  19. c     *****purpose:
  20. c     this subroutine grades the submatrices of a and b given by
  21. c     starting -1 low and ending -1 igh in the generalized
  22. c     eigenvalue problem a*x = (lambda)*b*x by permuting rows and
  23. c     columns such that the norm of the i-th row (column) of the
  24. c     a submatrix divided by the norm of the i-th row (column) of
  25. c     the b submatrix becomes smaller as i increases.
  26. c     ref.:  ward, r. c., balancing the generalized eigenvalue
  27. c     problem, siam j. sci. stat. comput., vol. 2, no. 2, june 1981,
  28. c     141-152.
  29. c
  30. c     *****parameter description:
  31. c
  32. c     on input:
  33. c
  34. c       ma,mb   integer
  35. c               row dimensions of the arrays containing matrices
  36. c               a and b respectively, as declared in the main calling
  37. c               program dimension statement;
  38. c
  39. c       n       integer
  40. c               order of the matrices a and b;
  41. c
  42. c       a       real(ma,n)
  43. c               contains the a matrix of the generalized eigenproblem
  44. c               defined above;
  45. c
  46. c       b       real(mb,n)
  47. c               contains the b matrix of the generalized eigenproblem
  48. c               defined above;
  49. c
  50. c       low     integer
  51. c               specifies the beginning -1 for the rows and
  52. c               columns of a and b to be graded;
  53. c
  54. c       igh     integer
  55. c               specifies the ending -1 for the rows and columns
  56. c               of a and b to be graded;
  57. c
  58. c       wk      real(n,2)
  59. c               work array that must contain at least 2*n locations.
  60. c               only locations low through igh and n+low through
  61. c               n+igh are referenced by this subroutine.
  62. c
  63. c     on output:
  64. c
  65. c       a,b     contain the permuted and graded a and b matrices;
  66. c
  67. c       cperm   real(n)
  68. c               contains in its low through igh locations the
  69. c               column permutations applied in grading the
  70. c               submatrices.  the other locations are not referenced
  71. c               by this subroutine;
  72. c
  73. c       wk      contains in its low through igh locations the row
  74. c               permutations applied in grading the submatrices.
  75. c
  76. c     *****algorithm notes:
  77. c     none.
  78. c
  79. c     *****history:
  80. c     written by r. c. ward.......
  81. c
  82. c     ---------------------------------------------------------------
  83. c
  84.       if (low .eq. igh) go to 510
  85.       ighm1 = igh-1
  86. c
  87. c     compute column norms of a / those of b
  88. c
  89.       do 420 j = low,igh
  90.          suma = 0.0d0
  91.          sumb = 0.0d0
  92.          do 410 i = low,igh
  93.             suma = suma + dabs(a(i,j))
  94.             sumb = sumb + dabs(b(i,j))
  95.   410    continue
  96.          if (sumb .eq. 0.0d0) go to 415
  97.          wk(j,2) = suma / sumb
  98.          go to 420
  99.   415    continue
  100.          wk(j,2) = 1.0d38
  101.   420 continue
  102. c
  103. c     permute columns to order them by decreasing quotients
  104. c
  105.       do 450 j = low,ighm1
  106.          cmax = wk(j,2)
  107.          jm = j
  108.          jp1 = j+1
  109.          do 430 k = jp1,igh
  110.             if (cmax .ge. wk(k,2)) go to 430
  111.             jm = k
  112.             cmax = wk(k,2)
  113.   430    continue
  114.          cperm(j) = jm
  115.          if (jm .eq. j) go to 450
  116.          temp = wk(j,2)
  117.          wk(j,2) = wk(jm,2)
  118.          wk(jm,2) = temp
  119.          do 440 i = 1,igh
  120.             temp = b(i,j)
  121.             b(i,j) = b(i,jm)
  122.             b(i,jm) = temp
  123.             temp = a(i,j)
  124.             a(i,j) = a(i,jm)
  125.             a(i,jm) = temp
  126.   440    continue
  127.   450 continue
  128.       cperm(igh) = igh
  129. c
  130. c     compute row norms of a / those of b
  131. c
  132.       do 470 i = low,igh
  133.          suma = 0.0d0
  134.          sumb = 0.0d0
  135.          do 460 j = low,igh
  136.             suma = suma + dabs(a(i,j))
  137.             sumb = sumb + dabs(b(i,j))
  138.   460    continue
  139.          if (sumb .eq. 0.0d0) go to 465
  140.          wk(i,2) = suma / sumb
  141.          go to 470
  142.   465    continue
  143.          wk(i,2) = 1.0d38
  144. c
  145. c     permute rows to order them by decreasing quotients
  146. c
  147.   470 continue
  148.       do 500 i = low,ighm1
  149.          rmax = wk(i,2)
  150.          im = i
  151.          ip1 = i+1
  152.          do 480 k = ip1,igh
  153.             if (rmax .ge. wk(k,2)) go to 480
  154.             im = k
  155.             rmax = wk(k,2)
  156.   480    continue
  157.          wk(i,1) = im
  158.          if (im .eq. i) go to 500
  159.          temp = wk(i,2)
  160.          wk(i,2) = wk(im,2)
  161.          wk(im,2) = temp
  162.          do 490 j = low,n
  163.             temp = b(i,j)
  164.             b(i,j) = b(im,j)
  165.             b(im,j) = temp
  166.             temp = a(i,j)
  167.             a(i,j) = a(im,j)
  168.             a(im,j) = temp
  169.   490    continue
  170.   500 continue
  171.       wk(igh,1) = igh
  172.   510 continue
  173.       return
  174. c
  175. c     last line of gradeq
  176. c
  177.       end
  178.