1 subroutine tridib(n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5) 2c 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) 7c 8c this subroutine is a translation of the algol procedure bisect, 9c num. math. 9, 386-393(1967) by barth, martin, and wilkinson. 10c handbook for auto. comp., vol.ii-linear algebra, 249-256(1971). 11c 12c this subroutine finds those eigenvalues of a tridiagonal 13c symmetric matrix between specified boundary indices, 14c using bisection. 15c 16c on input 17c 18c n is the order of the matrix. 19c 20c eps1 is an absolute error tolerance for the computed 21c eigenvalues. if the input eps1 is non-positive, 22c it is reset for each submatrix to a default value, 23c namely, minus the product of the relative machine 24c precision and the 1-norm of the submatrix. 25c 26c d contains the diagonal elements of the input matrix. 27c 28c e contains the subdiagonal elements of the input matrix 29c in its last n-1 positions. e(1) is arbitrary. 30c 31c e2 contains the squares of the corresponding elements of e. 32c e2(1) is arbitrary. 33c 34c m11 specifies the lower boundary index for the desired 35c eigenvalues. 36c 37c m specifies the number of eigenvalues desired. the upper 38c boundary index m22 is then obtained as m22=m11+m-1. 39c 40c on output 41c 42c eps1 is unaltered unless it has been reset to its 43c (last) default value. 44c 45c d and e are unaltered. 46c 47c elements of e2, corresponding to elements of e regarded 48c as negligible, have been replaced by zero causing the 49c matrix to split into a direct sum of submatrices. 50c e2(1) is also set to zero. 51c 52c lb and ub define an interval containing exactly the desired 53c eigenvalues. 54c 55c w contains, in its first m positions, the eigenvalues 56c between indices m11 and m22 in ascending order. 57c 58c ind contains in its first m positions the submatrix indices 59c associated with the corresponding eigenvalues in w -- 60c 1 for eigenvalues belonging to the first submatrix from 61c the top, 2 for those belonging to the second submatrix, etc.. 62c 63c ierr is set to 64c zero for normal return, 65c 3*n+1 if multiple eigenvalues at index m11 make 66c unique selection impossible, 67c 3*n+2 if multiple eigenvalues at index m22 make 68c unique selection impossible. 69c 70c rv4 and rv5 are temporary storage arrays. 71c 72c note that subroutine tql1, imtql1, or tqlrat is generally faster 73c than tridib, if more than n/4 eigenvalues are to be found. 74c 75c questions and comments should be directed to burton s. garbow, 76c mathematics and computer science div, argonne national laboratory 77c 78c this version dated august 1983. 79c 80c ------------------------------------------------------------------ 81c 82 ierr = 0 83 tag = 0 84 xu = d(1) 85 x0 = d(1) 86 u = 0.0d0 87c .......... look for small sub-diagonal entries and determine an 88c 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 101c 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 108c .......... determine an interval containing exactly 109c 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 135c .......... establish and process next submatrix, refining 136c 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 143c 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 155c 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 159c .......... 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 177c .......... find roots by bisection .......... 178 x0 = ub 179 isturm = 5 180c 181 do 240 i = m1, m2 182 rv5(i) = ub 183 rv4(i) = lb 184 240 continue 185c .......... loop for k-th eigenvalue 186c for k=m2 step -1 until m1 do -- 187c (-do- not used to legalize -computed go to-) .......... 188 k = m2 189 250 xu = lb 190c .......... 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 197c 198 280 if (x0 .gt. rv5(k)) x0 = rv5(k) 199c .......... 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 205c .......... in-line procedure for sturm sequence .......... 206 320 s = p - 1 207 u = 1.0d0 208c 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 218c 219 go to (60,80,200,220,360), isturm 220c .......... 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 231c .......... k-th eigenvalue found .......... 232 420 rv5(k) = x1 233 k = k - 1 234 if (k .ge. m1) go to 250 235c .......... order eigenvalues tagged with their 236c submatrix associations .......... 237 900 s = r 238 r = r + m2 - m1 + 1 239 j = 1 240 k = m1 241c 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 246c 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 252c 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 259c 260 940 if (q .lt. n) go to 100 261 go to 1001 262c .......... set error -- interval cannot be found containing 263c exactly the desired eigenvalues .......... 264 980 ierr = 3 * n + isturm 265 1001 lb = t1 266 ub = t2 267 return 268 end 269