home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / bqr.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  279 lines
  1.       subroutine bqr(nm,n,mb,a,t,r,ierr,nv,rv)
  2. c
  3.       integer i,j,k,l,m,n,ii,ik,jk,jm,kj,kk,km,ll,mb,mk,mn,mz,
  4.      x        m1,m2,m3,m4,ni,nm,nv,its,kj1,m21,m31,ierr,imult
  5.       double precision a(nm,mb),rv(nv)
  6.       double precision f,g,q,r,s,t,tst1,tst2,scale,pythag
  7. c
  8. c     this subroutine is a translation of the algol procedure bqr,
  9. c     num. math. 16, 85-92(1970) by martin, reinsch, and wilkinson.
  10. c     handbook for auto. comp., vol ii-linear algebra, 266-272(1971).
  11. c
  12. c     this subroutine finds the eigenvalue of smallest (usually)
  13. c     magnitude of a real symmetric band matrix using the
  14. c     qr algorithm with shifts of origin.  consecutive calls
  15. c     can be made to find further eigenvalues.
  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        mb is the (half) band width of the matrix, defined as the
  26. c          number of adjacent diagonals, including the principal
  27. c          diagonal, required to specify the non-zero portion of the
  28. c          lower triangle of the matrix.
  29. c
  30. c        a contains the lower triangle of the symmetric band input
  31. c          matrix stored as an n by mb array.  its lowest subdiagonal
  32. c          is stored in the last n+1-mb positions of the first column,
  33. c          its next subdiagonal in the last n+2-mb positions of the
  34. c          second column, further subdiagonals similarly, and finally
  35. c          its principal diagonal in the n positions of the last column.
  36. c          contents of storages not part of the matrix are arbitrary.
  37. c          on a subsequent call, its output contents from the previous
  38. c          call should be passed.
  39. c
  40. c        t specifies the shift (of eigenvalues) applied to the diagonal
  41. c          of a in forming the input matrix. what is actually determined
  42. c          is the eigenvalue of a+ti (i is the identity matrix) nearest
  43. c          to t.  on a subsequent call, the output value of t from the
  44. c          previous call should be passed if the next nearest eigenvalue
  45. c          is sought.
  46. c
  47. c        r should be specified as zero on the first call, and as its
  48. c          output value from the previous call on a subsequent call.
  49. c          it is used to determine when the last row and column of
  50. c          the transformed band matrix can be regarded as negligible.
  51. c
  52. c        nv must be set to the dimension of the array parameter rv
  53. c          as declared in the calling program dimension statement.
  54. c
  55. c     on output
  56. c
  57. c        a contains the transformed band matrix.  the matrix a+ti
  58. c          derived from the output parameters is similar to the
  59. c          input a+ti to within rounding errors.  its last row and
  60. c          column are null (if ierr is zero).
  61. c
  62. c        t contains the computed eigenvalue of a+ti (if ierr is zero).
  63. c
  64. c        r contains the maximum of its input value and the norm of the
  65. c          last column of the input matrix a.
  66. c
  67. c        ierr is set to
  68. c          zero       for normal return,
  69. c          n          if the eigenvalue has not been
  70. c                     determined after 30 iterations.
  71. c
  72. c        rv is a temporary storage array of dimension at least
  73. c          (2*mb**2+4*mb-3).  the first (3*mb-2) locations correspond
  74. c          to the algol array b, the next (2*mb-1) locations correspond
  75. c          to the algol array h, and the final (2*mb**2-mb) locations
  76. c          correspond to the mb by (2*mb-1) algol array u.
  77. c
  78. c     note. for a subsequent call, n should be replaced by n-1, but
  79. c     mb should not be altered even when it exceeds the current n.
  80. c
  81. c     calls pythag for  dsqrt(a*a + b*b) .
  82. c
  83. c     questions and comments should be directed to burton s. garbow,
  84. c     mathematics and computer science div, argonne national laboratory
  85. c
  86. c     this version dated august 1983.
  87. c
  88. c     ------------------------------------------------------------------
  89. c
  90.       ierr = 0
  91.       m1 = min0(mb,n)
  92.       m = m1 - 1
  93.       m2 = m + m
  94.       m21 = m2 + 1
  95.       m3 = m21 + m
  96.       m31 = m3 + 1
  97.       m4 = m31 + m2
  98.       mn = m + n
  99.       mz = mb - m1
  100.       its = 0
  101. c     .......... test for convergence ..........
  102.    40 g = a(n,mb)
  103.       if (m .eq. 0) go to 360
  104.       f = 0.0d0
  105. c
  106.       do 50 k = 1, m
  107.          mk = k + mz
  108.          f = f + dabs(a(n,mk))
  109.    50 continue
  110. c
  111.       if (its .eq. 0 .and. f .gt. r) r = f
  112.       tst1 = r
  113.       tst2 = tst1 + f
  114.       if (tst2 .le. tst1) go to 360
  115.       if (its .eq. 30) go to 1000
  116.       its = its + 1
  117. c     .......... form shift from bottom 2 by 2 minor ..........
  118.       if (f .gt. 0.25d0 * r .and. its .lt. 5) go to 90
  119.       f = a(n,mb-1)
  120.       if (f .eq. 0.0d0) go to 70
  121.       q = (a(n-1,mb) - g) / (2.0d0 * f)
  122.       s = pythag(q,1.0d0)
  123.       g = g - f / (q + dsign(s,q))
  124.    70 t = t + g
  125. c
  126.       do 80 i = 1, n
  127.    80 a(i,mb) = a(i,mb) - g
  128. c
  129.    90 do 100 k = m31, m4
  130.   100 rv(k) = 0.0d0
  131. c
  132.       do 350 ii = 1, mn
  133.          i = ii - m
  134.          ni = n - ii
  135.          if (ni .lt. 0) go to 230
  136. c     .......... form column of shifted matrix a-g*i ..........
  137.          l = max0(1,2-i)
  138. c
  139.          do 110 k = 1, m3
  140.   110    rv(k) = 0.0d0
  141. c
  142.          do 120 k = l, m1
  143.             km = k + m
  144.             mk = k + mz
  145.             rv(km) = a(ii,mk)
  146.   120    continue
  147. c
  148.          ll = min0(m,ni)
  149.          if (ll .eq. 0) go to 135
  150. c
  151.          do 130 k = 1, ll
  152.             km = k + m21
  153.             ik = ii + k
  154.             mk = mb - k
  155.             rv(km) = a(ik,mk)
  156.   130    continue
  157. c     .......... pre-multiply with householder reflections ..........
  158.   135    ll = m2
  159.          imult = 0
  160. c     .......... multiplication procedure ..........
  161.   140    kj = m4 - m1
  162. c
  163.          do 170 j = 1, ll
  164.             kj = kj + m1
  165.             jm = j + m3
  166.             if (rv(jm) .eq. 0.0d0) go to 170
  167.             f = 0.0d0
  168. c
  169.             do 150 k = 1, m1
  170.                kj = kj + 1
  171.                jk = j + k - 1
  172.                f = f + rv(kj) * rv(jk)
  173.   150       continue
  174. c
  175.             f = f / rv(jm)
  176.             kj = kj - m1
  177. c
  178.             do 160 k = 1, m1
  179.                kj = kj + 1
  180.                jk = j + k - 1
  181.                rv(jk) = rv(jk) - rv(kj) * f
  182.   160       continue
  183. c
  184.             kj = kj - m1
  185.   170    continue
  186. c
  187.          if (imult .ne. 0) go to 280
  188. c     .......... householder reflection ..........
  189.          f = rv(m21)
  190.          s = 0.0d0
  191.          rv(m4) = 0.0d0
  192.          scale = 0.0d0
  193. c
  194.          do 180 k = m21, m3
  195.   180    scale = scale + dabs(rv(k))
  196. c
  197.          if (scale .eq. 0.0d0) go to 210
  198. c
  199.          do 190 k = m21, m3
  200.   190    s = s + (rv(k)/scale)**2
  201. c
  202.          s = scale * scale * s
  203.          g = -dsign(dsqrt(s),f)
  204.          rv(m21) = g
  205.          rv(m4) = s - f * g
  206.          kj = m4 + m2 * m1 + 1
  207.          rv(kj) = f - g
  208. c
  209.          do 200 k = 2, m1
  210.             kj = kj + 1
  211.             km = k + m2
  212.             rv(kj) = rv(km)
  213.   200    continue
  214. c     .......... save column of triangular factor r ..........
  215.   210    do 220 k = l, m1
  216.             km = k + m
  217.             mk = k + mz
  218.             a(ii,mk) = rv(km)
  219.   220    continue
  220. c
  221.   230    l = max0(1,m1+1-i)
  222.          if (i .le. 0) go to 300
  223. c     .......... perform additional steps ..........
  224.          do 240 k = 1, m21
  225.   240    rv(k) = 0.0d0
  226. c
  227.          ll = min0(m1,ni+m1)
  228. c     .......... get row of triangular factor r ..........
  229.          do 250 kk = 1, ll
  230.             k = kk - 1
  231.             km = k + m1
  232.             ik = i + k
  233.             mk = mb - k
  234.             rv(km) = a(ik,mk)
  235.   250    continue
  236. c     .......... post-multiply with householder reflections ..........
  237.          ll = m1
  238.          imult = 1
  239.          go to 140
  240. c     .......... store column of new a matrix ..........
  241.   280    do 290 k = l, m1
  242.             mk = k + mz
  243.             a(i,mk) = rv(k)
  244.   290    continue
  245. c     .......... update householder reflections ..........
  246.   300    if (l .gt. 1) l = l - 1
  247.          kj1 = m4 + l * m1
  248. c
  249.          do 320 j = l, m2
  250.             jm = j + m3
  251.             rv(jm) = rv(jm+1)
  252. c
  253.             do 320 k = 1, m1
  254.                kj1 = kj1 + 1
  255.                kj = kj1 - m1
  256.                rv(kj) = rv(kj1)
  257.   320    continue
  258. c
  259.   350 continue
  260. c
  261.       go to 40
  262. c     .......... convergence ..........
  263.   360 t = t + g
  264. c
  265.       do 380 i = 1, n
  266.   380 a(i,mb) = a(i,mb) - g
  267. c
  268.       do 400 k = 1, m1
  269.          mk = k + mz
  270.          a(n,mk) = 0.0d0
  271.   400 continue
  272. c
  273.       go to 1001
  274. c     .......... set error -- no convergence to
  275. c                eigenvalue after 30 iterations ..........
  276.  1000 ierr = n
  277.  1001 return
  278.       end
  279.