home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / bisect.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  243 lines
  1.       subroutine bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5)
  2. c
  3.       integer i,j,k,l,m,n,p,q,r,s,ii,mm,m1,m2,tag,ierr,isturm
  4.       double precision d(n),e(n),e2(n),w(mm),rv4(n),rv5(n)
  5.       double precision u,v,lb,t1,t2,ub,xu,x0,x1,eps1,tst1,tst2,epslon
  6.       integer ind(mm)
  7. c
  8. c     this subroutine is a translation of the bisection technique
  9. c     in the algol procedure tristurm by peters and wilkinson.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 418-439(1971).
  11. c
  12. c     this subroutine finds those eigenvalues of a tridiagonal
  13. c     symmetric matrix which lie in a specified interval,
  14. c     using bisection.
  15. c
  16. c     on input
  17. c
  18. c        n is the order of the matrix.
  19. c
  20. c        eps1 is an absolute error tolerance for the computed
  21. c          eigenvalues.  if the input eps1 is non-positive,
  22. c          it is reset for each submatrix to a default value,
  23. c          namely, minus the product of the relative machine
  24. c          precision and the 1-norm of the submatrix.
  25. c
  26. c        d contains the diagonal elements of the input matrix.
  27. c
  28. c        e contains the subdiagonal elements of the input matrix
  29. c          in its last n-1 positions.  e(1) is arbitrary.
  30. c
  31. c        e2 contains the squares of the corresponding elements of e.
  32. c          e2(1) is arbitrary.
  33. c
  34. c        lb and ub define the interval to be searched for eigenvalues.
  35. c          if lb is not less than ub, no eigenvalues will be found.
  36. c
  37. c        mm should be set to an upper bound for the number of
  38. c          eigenvalues in the interval.  warning. if more than
  39. c          mm eigenvalues are determined to lie in the interval,
  40. c          an error return is made with no eigenvalues found.
  41. c
  42. c     on output
  43. c
  44. c        eps1 is unaltered unless it has been reset to its
  45. c          (last) default value.
  46. c
  47. c        d and e are unaltered.
  48. c
  49. c        elements of e2, corresponding to elements of e regarded
  50. c          as negligible, have been replaced by zero causing the
  51. c          matrix to split into a direct sum of submatrices.
  52. c          e2(1) is also set to zero.
  53. c
  54. c        m is the number of eigenvalues determined to lie in (lb,ub).
  55. c
  56. c        w contains the m eigenvalues in ascending order.
  57. c
  58. c        ind contains in its first m positions the submatrix indices
  59. c          associated with the corresponding eigenvalues in w --
  60. c          1 for eigenvalues belonging to the first submatrix from
  61. c          the top, 2 for those belonging to the second submatrix, etc..
  62. c
  63. c        ierr is set to
  64. c          zero       for normal return,
  65. c          3*n+1      if m exceeds mm.
  66. c
  67. c        rv4 and rv5 are temporary storage arrays.
  68. c
  69. c     the algol procedure sturmcnt contained in tristurm
  70. c     appears in bisect in-line.
  71. c
  72. c     note that subroutine tql1 or imtql1 is generally faster than
  73. c     bisect, if more than n/4 eigenvalues are to be found.
  74. c
  75. c     questions and comments should be directed to burton s. garbow,
  76. c     mathematics and computer science div, argonne national laboratory
  77. c
  78. c     this version dated august 1983.
  79. c
  80. c     ------------------------------------------------------------------
  81. c
  82.       ierr = 0
  83.       tag = 0
  84.       t1 = lb
  85.       t2 = ub
  86. c     .......... look for small sub-diagonal entries ..........
  87.       do 40 i = 1, n
  88.          if (i .eq. 1) go to 20
  89.          tst1 = dabs(d(i)) + dabs(d(i-1))
  90.          tst2 = tst1 + dabs(e(i))
  91.          if (tst2 .gt. tst1) go to 40
  92.    20    e2(i) = 0.0d0
  93.    40 continue
  94. c     .......... determine the number of eigenvalues
  95. c                in the interval ..........
  96.       p = 1
  97.       q = n
  98.       x1 = ub
  99.       isturm = 1
  100.       go to 320
  101.    60 m = s
  102.       x1 = lb
  103.       isturm = 2
  104.       go to 320
  105.    80 m = m - s
  106.       if (m .gt. mm) go to 980
  107.       q = 0
  108.       r = 0
  109. c     .......... establish and process next submatrix, refining
  110. c                interval by the gerschgorin bounds ..........
  111.   100 if (r .eq. m) go to 1001
  112.       tag = tag + 1
  113.       p = q + 1
  114.       xu = d(p)
  115.       x0 = d(p)
  116.       u = 0.0d0
  117. c
  118.       do 120 q = p, n
  119.          x1 = u
  120.          u = 0.0d0
  121.          v = 0.0d0
  122.          if (q .eq. n) go to 110
  123.          u = dabs(e(q+1))
  124.          v = e2(q+1)
  125.   110    xu = dmin1(d(q)-(x1+u),xu)
  126.          x0 = dmax1(d(q)+(x1+u),x0)
  127.          if (v .eq. 0.0d0) go to 140
  128.   120 continue
  129. c
  130.   140 x1 = epslon(dmax1(dabs(xu),dabs(x0)))
  131.       if (eps1 .le. 0.0d0) eps1 = -x1
  132.       if (p .ne. q) go to 180
  133. c     .......... check for isolated root within interval ..........
  134.       if (t1 .gt. d(p) .or. d(p) .ge. t2) go to 940
  135.       m1 = p
  136.       m2 = p
  137.       rv5(p) = d(p)
  138.       go to 900
  139.   180 x1 = x1 * (q - p + 1)
  140.       lb = dmax1(t1,xu-x1)
  141.       ub = dmin1(t2,x0+x1)
  142.       x1 = lb
  143.       isturm = 3
  144.       go to 320
  145.   200 m1 = s + 1
  146.       x1 = ub
  147.       isturm = 4
  148.       go to 320
  149.   220 m2 = s
  150.       if (m1 .gt. m2) go to 940
  151. c     .......... find roots by bisection ..........
  152.       x0 = ub
  153.       isturm = 5
  154. c
  155.       do 240 i = m1, m2
  156.          rv5(i) = ub
  157.          rv4(i) = lb
  158.   240 continue
  159. c     .......... loop for k-th eigenvalue
  160. c                for k=m2 step -1 until m1 do --
  161. c                (-do- not used to legalize -computed go to-) ..........
  162.       k = m2
  163.   250    xu = lb
  164. c     .......... for i=k step -1 until m1 do -- ..........
  165.          do 260 ii = m1, k
  166.             i = m1 + k - ii
  167.             if (xu .ge. rv4(i)) go to 260
  168.             xu = rv4(i)
  169.             go to 280
  170.   260    continue
  171. c
  172.   280    if (x0 .gt. rv5(k)) x0 = rv5(k)
  173. c     .......... next bisection step ..........
  174.   300    x1 = (xu + x0) * 0.5d0
  175.          if ((x0 - xu) .le. dabs(eps1)) go to 420
  176.          tst1 = 2.0d0 * (dabs(xu) + dabs(x0))
  177.          tst2 = tst1 + (x0 - xu)
  178.          if (tst2 .eq. tst1) go to 420
  179. c     .......... in-line procedure for sturm sequence ..........
  180.   320    s = p - 1
  181.          u = 1.0d0
  182. c
  183.          do 340 i = p, q
  184.             if (u .ne. 0.0d0) go to 325
  185.             v = dabs(e(i)) / epslon(1.0d0)
  186.             if (e2(i) .eq. 0.0d0) v = 0.0d0
  187.             go to 330
  188.   325       v = e2(i) / u
  189.   330       u = d(i) - x1 - v
  190.             if (u .lt. 0.0d0) s = s + 1
  191.   340    continue
  192. c
  193.          go to (60,80,200,220,360), isturm
  194. c     .......... refine intervals ..........
  195.   360    if (s .ge. k) go to 400
  196.          xu = x1
  197.          if (s .ge. m1) go to 380
  198.          rv4(m1) = x1
  199.          go to 300
  200.   380    rv4(s+1) = x1
  201.          if (rv5(s) .gt. x1) rv5(s) = x1
  202.          go to 300
  203.   400    x0 = x1
  204.          go to 300
  205. c     .......... k-th eigenvalue found ..........
  206.   420    rv5(k) = x1
  207.       k = k - 1
  208.       if (k .ge. m1) go to 250
  209. c     .......... order eigenvalues tagged with their
  210. c                submatrix associations ..........
  211.   900 s = r
  212.       r = r + m2 - m1 + 1
  213.       j = 1
  214.       k = m1
  215. c
  216.       do 920 l = 1, r
  217.          if (j .gt. s) go to 910
  218.          if (k .gt. m2) go to 940
  219.          if (rv5(k) .ge. w(l)) go to 915
  220. c
  221.          do 905 ii = j, s
  222.             i = l + s - ii
  223.             w(i+1) = w(i)
  224.             ind(i+1) = ind(i)
  225.   905    continue
  226. c
  227.   910    w(l) = rv5(k)
  228.          ind(l) = tag
  229.          k = k + 1
  230.          go to 920
  231.   915    j = j + 1
  232.   920 continue
  233. c
  234.   940 if (q .lt. n) go to 100
  235.       go to 1001
  236. c     .......... set error -- underestimate of number of
  237. c                eigenvalues in interval ..........
  238.   980 ierr = 3 * n + 1
  239.  1001 lb = t1
  240.       ub = t2
  241.       return
  242.       end
  243.