home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / bandr.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  265 lines
  1.       subroutine bandr(nm,n,mb,a,d,e,e2,matz,z)
  2. C  REFORMULATED S2 IN LOOP 500 TO AVOID OVERFLOW. (9/29/89 BSG)
  3. c
  4.       integer j,k,l,n,r,i1,i2,j1,j2,kr,mb,mr,m1,nm,n2,r1,ugl,maxl,maxr
  5.       double precision a(nm,mb),d(n),e(n),e2(n),z(nm,n)
  6.       double precision g,u,b1,b2,c2,f1,f2,s2,dmin,dminrt
  7.       logical matz
  8. c
  9. c     this subroutine is a translation of the algol procedure bandrd,
  10. c     num. math. 12, 231-241(1968) by schwarz.
  11. c     handbook for auto. comp., vol.ii-linear algebra, 273-283(1971).
  12. c
  13. c     this subroutine reduces a real symmetric band matrix
  14. c     to a symmetric tridiagonal matrix using and optionally
  15. c     accumulating orthogonal 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        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
  38. c        matz should be set to .true. if the transformation matrix is
  39. c          to be accumulated, and to .false. otherwise.
  40. c
  41. c     on output
  42. c
  43. c        a has been destroyed, except for its last two columns which
  44. c          contain a copy of the tridiagonal matrix.
  45. c
  46. c        d contains the diagonal elements of the tridiagonal matrix.
  47. c
  48. c        e contains the subdiagonal elements of the tridiagonal
  49. c          matrix in its last n-1 positions.  e(1) is set to zero.
  50. c
  51. c        e2 contains the squares of the corresponding elements of e.
  52. c          e2 may coincide with e if the squares are not needed.
  53. c
  54. c        z contains the orthogonal transformation matrix produced in
  55. c          the reduction if matz has been set to .true.  otherwise, z
  56. c          is not referenced.
  57. c
  58. c     questions and comments should be directed to burton s. garbow,
  59. c     mathematics and computer science div, argonne national laboratory
  60. c
  61. c     this version dated september 1989.
  62. c
  63. c     ------------------------------------------------------------------
  64. c
  65.       dmin = 2.0d0**(-64)
  66.       dminrt = 2.0d0**(-32)
  67. c     .......... initialize diagonal scaling matrix ..........
  68.       do 30 j = 1, n
  69.    30 d(j) = 1.0d0
  70. c
  71.       if (.not. matz) go to 60
  72. c
  73.       do 50 j = 1, n
  74. c
  75.          do 40 k = 1, n
  76.    40    z(j,k) = 0.0d0
  77. c
  78.          z(j,j) = 1.0d0
  79.    50 continue
  80. c
  81.    60 m1 = mb - 1
  82.       if (m1 - 1) 900, 800, 70
  83.    70 n2 = n - 2
  84. c
  85.       do 700 k = 1, n2
  86.          maxr = min0(m1,n-k)
  87. c     .......... for r=maxr step -1 until 2 do -- ..........
  88.          do 600 r1 = 2, maxr
  89.             r = maxr + 2 - r1
  90.             kr = k + r
  91.             mr = mb - r
  92.             g = a(kr,mr)
  93.             a(kr-1,1) = a(kr-1,mr+1)
  94.             ugl = k
  95. c
  96.             do 500 j = kr, n, m1
  97.                j1 = j - 1
  98.                j2 = j1 - 1
  99.                if (g .eq. 0.0d0) go to 600
  100.                b1 = a(j1,1) / g
  101.                b2 = b1 * d(j1) / d(j)
  102.                IF (ABS(B1) .GT. 1.0D0) THEN
  103.                   U = 1.0D0 / B1
  104.                   S2 = U / (U + B2)
  105.                ELSE 
  106.                   S2 = 1.0D0 / (1.0D0 + B1 * B2)
  107.                ENDIF
  108. c
  109.                if (s2 .ge. 0.5d0 ) go to 450
  110.                b1 = g / a(j1,1)
  111.                b2 = b1 * d(j) / d(j1)
  112.                c2 = 1.0d0 - s2
  113.                d(j1) = c2 * d(j1)
  114.                d(j) = c2 * d(j)
  115.                f1 = 2.0d0 * a(j,m1)
  116.                f2 = b1 * a(j1,mb)
  117.                a(j,m1) = -b2 * (b1 * a(j,m1) - a(j,mb)) - f2 + a(j,m1)
  118.                a(j1,mb) = b2 * (b2 * a(j,mb) + f1) + a(j1,mb)
  119.                a(j,mb) = b1 * (f2 - f1) + a(j,mb)
  120. c
  121.                do 200 l = ugl, j2
  122.                   i2 = mb - j + l
  123.                   u = a(j1,i2+1) + b2 * a(j,i2)
  124.                   a(j,i2) = -b1 * a(j1,i2+1) + a(j,i2)
  125.                   a(j1,i2+1) = u
  126.   200          continue
  127. c
  128.                ugl = j
  129.                a(j1,1) = a(j1,1) + b2 * g
  130.                if (j .eq. n) go to 350
  131.                maxl = min0(m1,n-j1)
  132. c
  133.                do 300 l = 2, maxl
  134.                   i1 = j1 + l
  135.                   i2 = mb - l
  136.                   u = a(i1,i2) + b2 * a(i1,i2+1)
  137.                   a(i1,i2+1) = -b1 * a(i1,i2) + a(i1,i2+1)
  138.                   a(i1,i2) = u
  139.   300          continue
  140. c
  141.                i1 = j + m1
  142.                if (i1 .gt. n) go to 350
  143.                g = b2 * a(i1,1)
  144.   350          if (.not. matz) go to 500
  145. c
  146.                do 400 l = 1, n
  147.                   u = z(l,j1) + b2 * z(l,j)
  148.                   z(l,j) = -b1 * z(l,j1) + z(l,j)
  149.                   z(l,j1) = u
  150.   400          continue
  151. c
  152.                go to 500
  153. c
  154.   450          u = d(j1)
  155.                d(j1) = s2 * d(j)
  156.                d(j) = s2 * u
  157.                f1 = 2.0d0 * a(j,m1)
  158.                f2 = b1 * a(j,mb)
  159.                u = b1 * (f2 - f1) + a(j1,mb)
  160.                a(j,m1) = b2 * (b1 * a(j,m1) - a(j1,mb)) + f2 - a(j,m1)
  161.                a(j1,mb) = b2 * (b2 * a(j1,mb) + f1) + a(j,mb)
  162.                a(j,mb) = u
  163. c
  164.                do 460 l = ugl, j2
  165.                   i2 = mb - j + l
  166.                   u = b2 * a(j1,i2+1) + a(j,i2)
  167.                   a(j,i2) = -a(j1,i2+1) + b1 * a(j,i2)
  168.                   a(j1,i2+1) = u
  169.   460          continue
  170. c
  171.                ugl = j
  172.                a(j1,1) = b2 * a(j1,1) + g
  173.                if (j .eq. n) go to 480
  174.                maxl = min0(m1,n-j1)
  175. c
  176.                do 470 l = 2, maxl
  177.                   i1 = j1 + l
  178.                   i2 = mb - l
  179.                   u = b2 * a(i1,i2) + a(i1,i2+1)
  180.                   a(i1,i2+1) = -a(i1,i2) + b1 * a(i1,i2+1)
  181.                   a(i1,i2) = u
  182.   470          continue
  183. c
  184.                i1 = j + m1
  185.                if (i1 .gt. n) go to 480
  186.                g = a(i1,1)
  187.                a(i1,1) = b1 * a(i1,1)
  188.   480          if (.not. matz) go to 500
  189. c
  190.                do 490 l = 1, n
  191.                   u = b2 * z(l,j1) + z(l,j)
  192.                   z(l,j) = -z(l,j1) + b1 * z(l,j)
  193.                   z(l,j1) = u
  194.   490          continue
  195. c
  196.   500       continue
  197. c
  198.   600    continue
  199. c
  200.          if (mod(k,64) .ne. 0) go to 700
  201. c     .......... rescale to avoid underflow or overflow ..........
  202.          do 650 j = k, n
  203.             if (d(j) .ge. dmin) go to 650
  204.             maxl = max0(1,mb+1-j)
  205. c
  206.             do 610 l = maxl, m1
  207.   610       a(j,l) = dminrt * a(j,l)
  208. c
  209.             if (j .eq. n) go to 630
  210.             maxl = min0(m1,n-j)
  211. c
  212.             do 620 l = 1, maxl
  213.                i1 = j + l
  214.                i2 = mb - l
  215.                a(i1,i2) = dminrt * a(i1,i2)
  216.   620       continue
  217. c
  218.   630       if (.not. matz) go to 645
  219. c
  220.             do 640 l = 1, n
  221.   640       z(l,j) = dminrt * z(l,j)
  222. c
  223.   645       a(j,mb) = dmin * a(j,mb)
  224.             d(j) = d(j) / dmin
  225.   650    continue
  226. c
  227.   700 continue
  228. c     .......... form square root of scaling matrix ..........
  229.   800 do 810 j = 2, n
  230.   810 e(j) = dsqrt(d(j))
  231. c
  232.       if (.not. matz) go to 840
  233. c
  234.       do 830 j = 1, n
  235. c
  236.          do 820 k = 2, n
  237.   820    z(j,k) = e(k) * z(j,k)
  238. c
  239.   830 continue
  240. c
  241.   840 u = 1.0d0
  242. c
  243.       do 850 j = 2, n
  244.          a(j,m1) = u * e(j) * a(j,m1)
  245.          u = e(j)
  246.          e2(j) = a(j,m1) ** 2
  247.          a(j,mb) = d(j) * a(j,mb)
  248.          d(j) = a(j,mb)
  249.          e(j) = a(j,m1)
  250.   850 continue
  251. c
  252.       d(1) = a(1,mb)
  253.       e(1) = 0.0d0
  254.       e2(1) = 0.0d0
  255.       go to 1001
  256. c
  257.   900 do 950 j = 1, n
  258.          d(j) = a(j,mb)
  259.          e(j) = 0.0d0
  260.          e2(j) = 0.0d0
  261.   950 continue
  262. c
  263.  1001 return
  264.       end
  265.