home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tridib.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  269 lines
  1.       subroutine tridib(n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5)
  2. c
  3.       integer i,j,k,l,m,n,p,q,r,s,ii,m1,m2,m11,m22,tag,ierr,isturm
  4.       double precision d(n),e(n),e2(n),w(m),rv4(n),rv5(n)
  5.       double precision u,v,lb,t1,t2,ub,xu,x0,x1,eps1,tst1,tst2,epslon
  6.       integer ind(m)
  7. c
  8. c     this subroutine is a translation of the algol procedure bisect,
  9. c     num. math. 9, 386-393(1967) by barth, martin, and wilkinson.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 249-256(1971).
  11. c
  12. c     this subroutine finds those eigenvalues of a tridiagonal
  13. c     symmetric matrix between specified boundary indices,
  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        m11 specifies the lower boundary index for the desired
  35. c          eigenvalues.
  36. c
  37. c        m specifies the number of eigenvalues desired.  the upper
  38. c          boundary index m22 is then obtained as m22=m11+m-1.
  39. c
  40. c     on output
  41. c
  42. c        eps1 is unaltered unless it has been reset to its
  43. c          (last) default value.
  44. c
  45. c        d and e are unaltered.
  46. c
  47. c        elements of e2, corresponding to elements of e regarded
  48. c          as negligible, have been replaced by zero causing the
  49. c          matrix to split into a direct sum of submatrices.
  50. c          e2(1) is also set to zero.
  51. c
  52. c        lb and ub define an interval containing exactly the desired
  53. c          eigenvalues.
  54. c
  55. c        w contains, in its first m positions, the eigenvalues
  56. c          between indices m11 and m22 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 multiple eigenvalues at index m11 make
  66. c                     unique selection impossible,
  67. c          3*n+2      if multiple eigenvalues at index m22 make
  68. c                     unique selection impossible.
  69. c
  70. c        rv4 and rv5 are temporary storage arrays.
  71. c
  72. c     note that subroutine tql1, imtql1, or tqlrat is generally faster
  73. c     than tridib, 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.       xu = d(1)
  85.       x0 = d(1)
  86.       u = 0.0d0
  87. c     .......... look for small sub-diagonal entries and determine an
  88. c                interval containing all the eigenvalues ..........
  89.       do 40 i = 1, n
  90.          x1 = u
  91.          u = 0.0d0
  92.          if (i .ne. n) u = dabs(e(i+1))
  93.          xu = dmin1(d(i)-(x1+u),xu)
  94.          x0 = dmax1(d(i)+(x1+u),x0)
  95.          if (i .eq. 1) go to 20
  96.          tst1 = dabs(d(i)) + dabs(d(i-1))
  97.          tst2 = tst1 + dabs(e(i))
  98.          if (tst2 .gt. tst1) go to 40
  99.    20    e2(i) = 0.0d0
  100.    40 continue
  101. c
  102.       x1 = n
  103.       x1 = x1 * epslon(dmax1(dabs(xu),dabs(x0)))
  104.       xu = xu - x1
  105.       t1 = xu
  106.       x0 = x0 + x1
  107.       t2 = x0
  108. c     .......... determine an interval containing exactly
  109. c                the desired eigenvalues ..........
  110.       p = 1
  111.       q = n
  112.       m1 = m11 - 1
  113.       if (m1 .eq. 0) go to 75
  114.       isturm = 1
  115.    50 v = x1
  116.       x1 = xu + (x0 - xu) * 0.5d0
  117.       if (x1 .eq. v) go to 980
  118.       go to 320
  119.    60 if (s - m1) 65, 73, 70
  120.    65 xu = x1
  121.       go to 50
  122.    70 x0 = x1
  123.       go to 50
  124.    73 xu = x1
  125.       t1 = x1
  126.    75 m22 = m1 + m
  127.       if (m22 .eq. n) go to 90
  128.       x0 = t2
  129.       isturm = 2
  130.       go to 50
  131.    80 if (s - m22) 65, 85, 70
  132.    85 t2 = x1
  133.    90 q = 0
  134.       r = 0
  135. c     .......... establish and process next submatrix, refining
  136. c                interval by the gerschgorin bounds ..........
  137.   100 if (r .eq. m) go to 1001
  138.       tag = tag + 1
  139.       p = q + 1
  140.       xu = d(p)
  141.       x0 = d(p)
  142.       u = 0.0d0
  143. c
  144.       do 120 q = p, n
  145.          x1 = u
  146.          u = 0.0d0
  147.          v = 0.0d0
  148.          if (q .eq. n) go to 110
  149.          u = dabs(e(q+1))
  150.          v = e2(q+1)
  151.   110    xu = dmin1(d(q)-(x1+u),xu)
  152.          x0 = dmax1(d(q)+(x1+u),x0)
  153.          if (v .eq. 0.0d0) go to 140
  154.   120 continue
  155. c
  156.   140 x1 = epslon(dmax1(dabs(xu),dabs(x0)))
  157.       if (eps1 .le. 0.0d0) eps1 = -x1
  158.       if (p .ne. q) go to 180
  159. c     .......... check for isolated root within interval ..........
  160.       if (t1 .gt. d(p) .or. d(p) .ge. t2) go to 940
  161.       m1 = p
  162.       m2 = p
  163.       rv5(p) = d(p)
  164.       go to 900
  165.   180 x1 = x1 * (q - p + 1)
  166.       lb = dmax1(t1,xu-x1)
  167.       ub = dmin1(t2,x0+x1)
  168.       x1 = lb
  169.       isturm = 3
  170.       go to 320
  171.   200 m1 = s + 1
  172.       x1 = ub
  173.       isturm = 4
  174.       go to 320
  175.   220 m2 = s
  176.       if (m1 .gt. m2) go to 940
  177. c     .......... find roots by bisection ..........
  178.       x0 = ub
  179.       isturm = 5
  180. c
  181.       do 240 i = m1, m2
  182.          rv5(i) = ub
  183.          rv4(i) = lb
  184.   240 continue
  185. c     .......... loop for k-th eigenvalue
  186. c                for k=m2 step -1 until m1 do --
  187. c                (-do- not used to legalize -computed go to-) ..........
  188.       k = m2
  189.   250    xu = lb
  190. c     .......... for i=k step -1 until m1 do -- ..........
  191.          do 260 ii = m1, k
  192.             i = m1 + k - ii
  193.             if (xu .ge. rv4(i)) go to 260
  194.             xu = rv4(i)
  195.             go to 280
  196.   260    continue
  197. c
  198.   280    if (x0 .gt. rv5(k)) x0 = rv5(k)
  199. c     .......... next bisection step ..........
  200.   300    x1 = (xu + x0) * 0.5d0
  201.          if ((x0 - xu) .le. dabs(eps1)) go to 420
  202.          tst1 = 2.0d0 * (dabs(xu) + dabs(x0))
  203.          tst2 = tst1 + (x0 - xu)
  204.          if (tst2 .eq. tst1) go to 420
  205. c     .......... in-line procedure for sturm sequence ..........
  206.   320    s = p - 1
  207.          u = 1.0d0
  208. c
  209.          do 340 i = p, q
  210.             if (u .ne. 0.0d0) go to 325
  211.             v = dabs(e(i)) / epslon(1.0d0)
  212.             if (e2(i) .eq. 0.0d0) v = 0.0d0
  213.             go to 330
  214.   325       v = e2(i) / u
  215.   330       u = d(i) - x1 - v
  216.             if (u .lt. 0.0d0) s = s + 1
  217.   340    continue
  218. c
  219.          go to (60,80,200,220,360), isturm
  220. c     .......... refine intervals ..........
  221.   360    if (s .ge. k) go to 400
  222.          xu = x1
  223.          if (s .ge. m1) go to 380
  224.          rv4(m1) = x1
  225.          go to 300
  226.   380    rv4(s+1) = x1
  227.          if (rv5(s) .gt. x1) rv5(s) = x1
  228.          go to 300
  229.   400    x0 = x1
  230.          go to 300
  231. c     .......... k-th eigenvalue found ..........
  232.   420    rv5(k) = x1
  233.       k = k - 1
  234.       if (k .ge. m1) go to 250
  235. c     .......... order eigenvalues tagged with their
  236. c                submatrix associations ..........
  237.   900 s = r
  238.       r = r + m2 - m1 + 1
  239.       j = 1
  240.       k = m1
  241. c
  242.       do 920 l = 1, r
  243.          if (j .gt. s) go to 910
  244.          if (k .gt. m2) go to 940
  245.          if (rv5(k) .ge. w(l)) go to 915
  246. c
  247.          do 905 ii = j, s
  248.             i = l + s - ii
  249.             w(i+1) = w(i)
  250.             ind(i+1) = ind(i)
  251.   905    continue
  252. c
  253.   910    w(l) = rv5(k)
  254.          ind(l) = tag
  255.          k = k + 1
  256.          go to 920
  257.   915    j = j + 1
  258.   920 continue
  259. c
  260.   940 if (q .lt. n) go to 100
  261.       go to 1001
  262. c     .......... set error -- interval cannot be found containing
  263. c                exactly the desired eigenvalues ..........
  264.   980 ierr = 3 * n + isturm
  265.  1001 lb = t1
  266.       ub = t2
  267.       return
  268.       end
  269.