subroutine bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5) c integer i,j,k,l,m,n,p,q,r,s,ii,mm,m1,m2,tag,ierr,isturm double precision d(n),e(n),e2(n),w(mm),rv4(n),rv5(n) double precision u,v,lb,t1,t2,ub,xu,x0,x1,eps1,tst1,tst2,epslon integer ind(mm) c c this subroutine is a translation of the bisection technique c in the algol procedure tristurm by peters and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 418-439(1971). c c this subroutine finds those eigenvalues of a tridiagonal c symmetric matrix which lie in a specified interval, c using bisection. c c on input c c n is the order of the matrix. c c eps1 is an absolute error tolerance for the computed c eigenvalues. if the input eps1 is non-positive, c it is reset for each submatrix to a default value, c namely, minus the product of the relative machine c precision and the 1-norm of the submatrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c e2 contains the squares of the corresponding elements of e. c e2(1) is arbitrary. c c lb and ub define the interval to be searched for eigenvalues. c if lb is not less than ub, no eigenvalues will be found. c c mm should be set to an upper bound for the number of c eigenvalues in the interval. warning. if more than c mm eigenvalues are determined to lie in the interval, c an error return is made with no eigenvalues found. c c on output c c eps1 is unaltered unless it has been reset to its c (last) default value. c c d and e are unaltered. c c elements of e2, corresponding to elements of e regarded c as negligible, have been replaced by zero causing the c matrix to split into a direct sum of submatrices. c e2(1) is also set to zero. c c m is the number of eigenvalues determined to lie in (lb,ub). c c w contains the m eigenvalues in ascending order. c c ind contains in its first m positions the submatrix indices c associated with the corresponding eigenvalues in w -- c 1 for eigenvalues belonging to the first submatrix from c the top, 2 for those belonging to the second submatrix, etc.. c c ierr is set to c zero for normal return, c 3*n+1 if m exceeds mm. c c rv4 and rv5 are temporary storage arrays. c c the algol procedure sturmcnt contained in tristurm c appears in bisect in-line. c c note that subroutine tql1 or imtql1 is generally faster than c bisect, if more than n/4 eigenvalues are to be found. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 tag = 0 t1 = lb t2 = ub c .......... look for small sub-diagonal entries .......... do 40 i = 1, n if (i .eq. 1) go to 20 tst1 = dabs(d(i)) + dabs(d(i-1)) tst2 = tst1 + dabs(e(i)) if (tst2 .gt. tst1) go to 40 20 e2(i) = 0.0d0 40 continue c .......... determine the number of eigenvalues c in the interval .......... p = 1 q = n x1 = ub isturm = 1 go to 320 60 m = s x1 = lb isturm = 2 go to 320 80 m = m - s if (m .gt. mm) go to 980 q = 0 r = 0 c .......... establish and process next submatrix, refining c interval by the gerschgorin bounds .......... 100 if (r .eq. m) go to 1001 tag = tag + 1 p = q + 1 xu = d(p) x0 = d(p) u = 0.0d0 c do 120 q = p, n x1 = u u = 0.0d0 v = 0.0d0 if (q .eq. n) go to 110 u = dabs(e(q+1)) v = e2(q+1) 110 xu = dmin1(d(q)-(x1+u),xu) x0 = dmax1(d(q)+(x1+u),x0) if (v .eq. 0.0d0) go to 140 120 continue c 140 x1 = epslon(dmax1(dabs(xu),dabs(x0))) if (eps1 .le. 0.0d0) eps1 = -x1 if (p .ne. q) go to 180 c .......... check for isolated root within interval .......... if (t1 .gt. d(p) .or. d(p) .ge. t2) go to 940 m1 = p m2 = p rv5(p) = d(p) go to 900 180 x1 = x1 * (q - p + 1) lb = dmax1(t1,xu-x1) ub = dmin1(t2,x0+x1) x1 = lb isturm = 3 go to 320 200 m1 = s + 1 x1 = ub isturm = 4 go to 320 220 m2 = s if (m1 .gt. m2) go to 940 c .......... find roots by bisection .......... x0 = ub isturm = 5 c do 240 i = m1, m2 rv5(i) = ub rv4(i) = lb 240 continue c .......... loop for k-th eigenvalue c for k=m2 step -1 until m1 do -- c (-do- not used to legalize -computed go to-) .......... k = m2 250 xu = lb c .......... for i=k step -1 until m1 do -- .......... do 260 ii = m1, k i = m1 + k - ii if (xu .ge. rv4(i)) go to 260 xu = rv4(i) go to 280 260 continue c 280 if (x0 .gt. rv5(k)) x0 = rv5(k) c .......... next bisection step .......... 300 x1 = (xu + x0) * 0.5d0 if ((x0 - xu) .le. dabs(eps1)) go to 420 tst1 = 2.0d0 * (dabs(xu) + dabs(x0)) tst2 = tst1 + (x0 - xu) if (tst2 .eq. tst1) go to 420 c .......... in-line procedure for sturm sequence .......... 320 s = p - 1 u = 1.0d0 c do 340 i = p, q if (u .ne. 0.0d0) go to 325 v = dabs(e(i)) / epslon(1.0d0) if (e2(i) .eq. 0.0d0) v = 0.0d0 go to 330 325 v = e2(i) / u 330 u = d(i) - x1 - v if (u .lt. 0.0d0) s = s + 1 340 continue c go to (60,80,200,220,360), isturm c .......... refine intervals .......... 360 if (s .ge. k) go to 400 xu = x1 if (s .ge. m1) go to 380 rv4(m1) = x1 go to 300 380 rv4(s+1) = x1 if (rv5(s) .gt. x1) rv5(s) = x1 go to 300 400 x0 = x1 go to 300 c .......... k-th eigenvalue found .......... 420 rv5(k) = x1 k = k - 1 if (k .ge. m1) go to 250 c .......... order eigenvalues tagged with their c submatrix associations .......... 900 s = r r = r + m2 - m1 + 1 j = 1 k = m1 c do 920 l = 1, r if (j .gt. s) go to 910 if (k .gt. m2) go to 940 if (rv5(k) .ge. w(l)) go to 915 c do 905 ii = j, s i = l + s - ii w(i+1) = w(i) ind(i+1) = ind(i) 905 continue c 910 w(l) = rv5(k) ind(l) = tag k = k + 1 go to 920 915 j = j + 1 920 continue c 940 if (q .lt. n) go to 100 go to 1001 c .......... set error -- underestimate of number of c eigenvalues in interval .......... 980 ierr = 3 * n + 1 1001 lb = t1 ub = t2 return end