home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tsturm.f < prev   
Text File  |  1996-09-28  |  11KB  |  369 lines
  1.       subroutine tsturm(nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,
  2.      x                  ierr,rv1,rv2,rv3,rv4,rv5,rv6)
  3. c
  4.       integer i,j,k,m,n,p,q,r,s,ii,ip,jj,mm,m1,m2,nm,its,
  5.      x        ierr,group,isturm
  6.       double precision d(n),e(n),e2(n),w(mm),z(nm,mm),
  7.      x       rv1(n),rv2(n),rv3(n),rv4(n),rv5(n),rv6(n)
  8.       double precision u,v,lb,t1,t2,ub,uk,xu,x0,x1,eps1,eps2,eps3,eps4,
  9.      x       norm,tst1,tst2,epslon,pythag
  10. c
  11. c     this subroutine is a translation of the algol procedure tristurm
  12. c     by peters and wilkinson.
  13. c     handbook for auto. comp., vol.ii-linear algebra, 418-439(1971).
  14. c
  15. c     this subroutine finds those eigenvalues of a tridiagonal
  16. c     symmetric matrix which lie in a specified interval and their
  17. c     associated eigenvectors, using bisection and inverse iteration.
  18. c
  19. c     on input
  20. c
  21. c        nm must be set to the row dimension of two-dimensional
  22. c          array parameters as declared in the calling program
  23. c          dimension statement.
  24. c
  25. c        n is the order of the matrix.
  26. c
  27. c        eps1 is an absolute error tolerance for the computed
  28. c          eigenvalues.  it should be chosen commensurate with
  29. c          relative perturbations in the matrix elements of the
  30. c          order of the relative machine precision.  if the
  31. c          input eps1 is non-positive, it is reset for each
  32. c          submatrix to a default value, namely, minus the
  33. c          product of the relative machine precision and the
  34. c          1-norm of the submatrix.
  35. c
  36. c        d contains the diagonal elements of the input matrix.
  37. c
  38. c        e contains the subdiagonal elements of the input matrix
  39. c          in its last n-1 positions.  e(1) is arbitrary.
  40. c
  41. c        e2 contains the squares of the corresponding elements of e.
  42. c          e2(1) is arbitrary.
  43. c
  44. c        lb and ub define the interval to be searched for eigenvalues.
  45. c          if lb is not less than ub, no eigenvalues will be found.
  46. c
  47. c        mm should be set to an upper bound for the number of
  48. c          eigenvalues in the interval.  warning. if more than
  49. c          mm eigenvalues are determined to lie in the interval,
  50. c          an error return is made with no values or vectors found.
  51. c
  52. c     on output
  53. c
  54. c        eps1 is unaltered unless it has been reset to its
  55. c          (last) default value.
  56. c
  57. c        d and e are unaltered.
  58. c
  59. c        elements of e2, corresponding to elements of e regarded
  60. c          as negligible, have been replaced by zero causing the
  61. c          matrix to split into a direct sum of submatrices.
  62. c          e2(1) is also set to zero.
  63. c
  64. c        m is the number of eigenvalues determined to lie in (lb,ub).
  65. c
  66. c        w contains the m eigenvalues in ascending order if the matrix
  67. c          does not split.  if the matrix splits, the eigenvalues are
  68. c          in ascending order for each submatrix.  if a vector error
  69. c          exit is made, w contains those values already found.
  70. c
  71. c        z contains the associated set of orthonormal eigenvectors.
  72. c          if an error exit is made, z contains those vectors
  73. c          already found.
  74. c
  75. c        ierr is set to
  76. c          zero       for normal return,
  77. c          3*n+1      if m exceeds mm.
  78. c          4*n+r      if the eigenvector corresponding to the r-th
  79. c                     eigenvalue fails to converge in 5 iterations.
  80. c
  81. c        rv1, rv2, rv3, rv4, rv5, and rv6 are temporary storage arrays.
  82. c
  83. c     the algol procedure sturmcnt contained in tristurm
  84. c     appears in tsturm in-line.
  85. c
  86. c     calls pythag for  dsqrt(a*a + b*b) .
  87. c
  88. c     questions and comments should be directed to burton s. garbow,
  89. c     mathematics and computer science div, argonne national laboratory
  90. c
  91. c     this version dated august 1983.
  92. c
  93. c     ------------------------------------------------------------------
  94. c
  95.       ierr = 0
  96.       t1 = lb
  97.       t2 = ub
  98. c     .......... look for small sub-diagonal entries ..........
  99.       do 40 i = 1, n
  100.          if (i .eq. 1) go to 20
  101.          tst1 = dabs(d(i)) + dabs(d(i-1))
  102.          tst2 = tst1 + dabs(e(i))
  103.          if (tst2 .gt. tst1) go to 40
  104.    20    e2(i) = 0.0d0
  105.    40 continue
  106. c     .......... determine the number of eigenvalues
  107. c                in the interval ..........
  108.       p = 1
  109.       q = n
  110.       x1 = ub
  111.       isturm = 1
  112.       go to 320
  113.    60 m = s
  114.       x1 = lb
  115.       isturm = 2
  116.       go to 320
  117.    80 m = m - s
  118.       if (m .gt. mm) go to 980
  119.       q = 0
  120.       r = 0
  121. c     .......... establish and process next submatrix, refining
  122. c                interval by the gerschgorin bounds ..........
  123.   100 if (r .eq. m) go to 1001
  124.       p = q + 1
  125.       xu = d(p)
  126.       x0 = d(p)
  127.       u = 0.0d0
  128. c
  129.       do 120 q = p, n
  130.          x1 = u
  131.          u = 0.0d0
  132.          v = 0.0d0
  133.          if (q .eq. n) go to 110
  134.          u = dabs(e(q+1))
  135.          v = e2(q+1)
  136.   110    xu = dmin1(d(q)-(x1+u),xu)
  137.          x0 = dmax1(d(q)+(x1+u),x0)
  138.          if (v .eq. 0.0d0) go to 140
  139.   120 continue
  140. c
  141.   140 x1 = epslon(dmax1(dabs(xu),dabs(x0)))
  142.       if (eps1 .le. 0.0d0) eps1 = -x1
  143.       if (p .ne. q) go to 180
  144. c     .......... check for isolated root within interval ..........
  145.       if (t1 .gt. d(p) .or. d(p) .ge. t2) go to 940
  146.       r = r + 1
  147. c
  148.       do 160 i = 1, n
  149.   160 z(i,r) = 0.0d0
  150. c
  151.       w(r) = d(p)
  152.       z(p,r) = 1.0d0
  153.       go to 940
  154.   180 u = q-p+1
  155.       x1 = u * x1
  156.       lb = dmax1(t1,xu-x1)
  157.       ub = dmin1(t2,x0+x1)
  158.       x1 = lb
  159.       isturm = 3
  160.       go to 320
  161.   200 m1 = s + 1
  162.       x1 = ub
  163.       isturm = 4
  164.       go to 320
  165.   220 m2 = s
  166.       if (m1 .gt. m2) go to 940
  167. c     .......... find roots by bisection ..........
  168.       x0 = ub
  169.       isturm = 5
  170. c
  171.       do 240 i = m1, m2
  172.          rv5(i) = ub
  173.          rv4(i) = lb
  174.   240 continue
  175. c     .......... loop for k-th eigenvalue
  176. c                for k=m2 step -1 until m1 do --
  177. c                (-do- not used to legalize -computed go to-) ..........
  178.       k = m2
  179.   250    xu = lb
  180. c     .......... for i=k step -1 until m1 do -- ..........
  181.          do 260 ii = m1, k
  182.             i = m1 + k - ii
  183.             if (xu .ge. rv4(i)) go to 260
  184.             xu = rv4(i)
  185.             go to 280
  186.   260    continue
  187. c
  188.   280    if (x0 .gt. rv5(k)) x0 = rv5(k)
  189. c     .......... next bisection step ..........
  190.   300    x1 = (xu + x0) * 0.5d0
  191.          if ((x0 - xu) .le. dabs(eps1)) go to 420
  192.          tst1 = 2.0d0 * (dabs(xu) + dabs(x0))
  193.          tst2 = tst1 + (x0 - xu)
  194.          if (tst2 .eq. tst1) go to 420
  195. c     .......... in-line procedure for sturm sequence ..........
  196.   320    s = p - 1
  197.          u = 1.0d0
  198. c
  199.          do 340 i = p, q
  200.             if (u .ne. 0.0d0) go to 325
  201.             v = dabs(e(i)) / epslon(1.0d0)
  202.             if (e2(i) .eq. 0.0d0) v = 0.0d0
  203.             go to 330
  204.   325       v = e2(i) / u
  205.   330       u = d(i) - x1 - v
  206.             if (u .lt. 0.0d0) s = s + 1
  207.   340    continue
  208. c
  209.          go to (60,80,200,220,360), isturm
  210. c     .......... refine intervals ..........
  211.   360    if (s .ge. k) go to 400
  212.          xu = x1
  213.          if (s .ge. m1) go to 380
  214.          rv4(m1) = x1
  215.          go to 300
  216.   380    rv4(s+1) = x1
  217.          if (rv5(s) .gt. x1) rv5(s) = x1
  218.          go to 300
  219.   400    x0 = x1
  220.          go to 300
  221. c     .......... k-th eigenvalue found ..........
  222.   420    rv5(k) = x1
  223.       k = k - 1
  224.       if (k .ge. m1) go to 250
  225. c     .......... find vectors by inverse iteration ..........
  226.       norm = dabs(d(p))
  227.       ip = p + 1
  228. c
  229.       do 500 i = ip, q
  230.   500 norm = dmax1(norm, dabs(d(i)) + dabs(e(i)))
  231. c     .......... eps2 is the criterion for grouping,
  232. c                eps3 replaces zero pivots and equal
  233. c                roots are modified by eps3,
  234. c                eps4 is taken very small to avoid overflow ..........
  235.       eps2 = 1.0d-3 * norm
  236.       eps3 = epslon(norm)
  237.       uk = q - p + 1
  238.       eps4 = uk * eps3
  239.       uk = eps4 / dsqrt(uk)
  240.       group = 0
  241.       s = p
  242. c
  243.       do 920 k = m1, m2
  244.          r = r + 1
  245.          its = 1
  246.          w(r) = rv5(k)
  247.          x1 = rv5(k)
  248. c     .......... look for close or coincident roots ..........
  249.          if (k .eq. m1) go to 520
  250.          if (x1 - x0 .ge. eps2) group = -1
  251.          group = group + 1
  252.          if (x1 .le. x0) x1 = x0 + eps3
  253. c     .......... elimination with interchanges and
  254. c                initialization of vector ..........
  255.   520    v = 0.0d0
  256. c
  257.          do 580 i = p, q
  258.             rv6(i) = uk
  259.             if (i .eq. p) go to 560
  260.             if (dabs(e(i)) .lt. dabs(u)) go to 540
  261.             xu = u / e(i)
  262.             rv4(i) = xu
  263.             rv1(i-1) = e(i)
  264.             rv2(i-1) = d(i) - x1
  265.             rv3(i-1) = 0.0d0
  266.             if (i .ne. q) rv3(i-1) = e(i+1)
  267.             u = v - xu * rv2(i-1)
  268.             v = -xu * rv3(i-1)
  269.             go to 580
  270.   540       xu = e(i) / u
  271.             rv4(i) = xu
  272.             rv1(i-1) = u
  273.             rv2(i-1) = v
  274.             rv3(i-1) = 0.0d0
  275.   560       u = d(i) - x1 - xu * v
  276.             if (i .ne. q) v = e(i+1)
  277.   580    continue
  278. c
  279.          if (u .eq. 0.0d0) u = eps3
  280.          rv1(q) = u
  281.          rv2(q) = 0.0d0
  282.          rv3(q) = 0.0d0
  283. c     .......... back substitution
  284. c                for i=q step -1 until p do -- ..........
  285.   600    do 620 ii = p, q
  286.             i = p + q - ii
  287.             rv6(i) = (rv6(i) - u * rv2(i) - v * rv3(i)) / rv1(i)
  288.             v = u
  289.             u = rv6(i)
  290.   620    continue
  291. c     .......... orthogonalize with respect to previous
  292. c                members of group ..........
  293.          if (group .eq. 0) go to 700
  294. c
  295.          do 680 jj = 1, group
  296.             j = r - group - 1 + jj
  297.             xu = 0.0d0
  298. c
  299.             do 640 i = p, q
  300.   640       xu = xu + rv6(i) * z(i,j)
  301. c
  302.             do 660 i = p, q
  303.   660       rv6(i) = rv6(i) - xu * z(i,j)
  304. c
  305.   680    continue
  306. c
  307.   700    norm = 0.0d0
  308. c
  309.          do 720 i = p, q
  310.   720    norm = norm + dabs(rv6(i))
  311. c
  312.          if (norm .ge. 1.0d0) go to 840
  313. c     .......... forward substitution ..........
  314.          if (its .eq. 5) go to 960
  315.          if (norm .ne. 0.0d0) go to 740
  316.          rv6(s) = eps4
  317.          s = s + 1
  318.          if (s .gt. q) s = p
  319.          go to 780
  320.   740    xu = eps4 / norm
  321. c
  322.          do 760 i = p, q
  323.   760    rv6(i) = rv6(i) * xu
  324. c     .......... elimination operations on next vector
  325. c                iterate ..........
  326.   780    do 820 i = ip, q
  327.             u = rv6(i)
  328. c     .......... if rv1(i-1) .eq. e(i), a row interchange
  329. c                was performed earlier in the
  330. c                triangularization process ..........
  331.             if (rv1(i-1) .ne. e(i)) go to 800
  332.             u = rv6(i-1)
  333.             rv6(i-1) = rv6(i)
  334.   800       rv6(i) = u - rv4(i) * rv6(i-1)
  335.   820    continue
  336. c
  337.          its = its + 1
  338.          go to 600
  339. c     .......... normalize so that sum of squares is
  340. c                1 and expand to full order ..........
  341.   840    u = 0.0d0
  342. c
  343.          do 860 i = p, q
  344.   860    u = pythag(u,rv6(i))
  345. c
  346.          xu = 1.0d0 / u
  347. c
  348.          do 880 i = 1, n
  349.   880    z(i,r) = 0.0d0
  350. c
  351.          do 900 i = p, q
  352.   900    z(i,r) = rv6(i) * xu
  353. c
  354.          x0 = x1
  355.   920 continue
  356. c
  357.   940 if (q .lt. n) go to 100
  358.       go to 1001
  359. c     .......... set error -- non-converged eigenvector ..........
  360.   960 ierr = 4 * n + r
  361.       go to 1001
  362. c     .......... set error -- underestimate of number of
  363. c                eigenvalues in interval ..........
  364.   980 ierr = 3 * n + 1
  365.  1001 lb = t1
  366.       ub = t2
  367.       return
  368.       end
  369.