subroutine imtqlv(n,d,e,e2,w,ind,ierr,rv1) c integer i,j,k,l,m,n,ii,mml,tag,ierr double precision d(n),e(n),e2(n),w(n),rv1(n) double precision b,c,f,g,p,r,s,tst1,tst2,pythag integer ind(n) c c this subroutine is a variant of imtql1 which is a translation of c algol procedure imtql1, num. math. 12, 377-383(1968) by martin and c wilkinson, as modified in num. math. 15, 450(1970) by dubrulle. c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). c c this subroutine finds the eigenvalues of a symmetric tridiagonal c matrix by the implicit ql method and associates with them c their corresponding submatrix indices. c c on input c c n is the order of the matrix. 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 on output 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 w contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c ind contains the submatrix indices associated with the c corresponding eigenvalues in w -- 1 for eigenvalues c belonging to the first submatrix from the top, c 2 for those belonging to the second submatrix, etc.. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c rv1 is a temporary storage array. c c calls pythag for dsqrt(a*a + b*b) . 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 k = 0 tag = 0 c do 100 i = 1, n w(i) = d(i) if (i .ne. 1) rv1(i-1) = e(i) 100 continue c e2(1) = 0.0d0 rv1(n) = 0.0d0 c do 290 l = 1, n j = 0 c .......... look for small sub-diagonal element .......... 105 do 110 m = l, n if (m .eq. n) go to 120 tst1 = dabs(w(m)) + dabs(w(m+1)) tst2 = tst1 + dabs(rv1(m)) if (tst2 .eq. tst1) go to 120 c .......... guard against underflowed element of e2 .......... if (e2(m+1) .eq. 0.0d0) go to 125 110 continue c 120 if (m .le. k) go to 130 if (m .ne. n) e2(m+1) = 0.0d0 125 k = m tag = tag + 1 130 p = w(l) if (m .eq. l) go to 215 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... g = (w(l+1) - p) / (2.0d0 * rv1(l)) r = pythag(g,1.0d0) g = w(m) - p + rv1(l) / (g + dsign(r,g)) s = 1.0d0 c = 1.0d0 p = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii f = s * rv1(i) b = c * rv1(i) r = pythag(f,g) rv1(i+1) = r if (r .eq. 0.0d0) go to 210 s = f / r c = g / r g = w(i+1) - p r = (w(i) - g) * s + 2.0d0 * c * b p = s * r w(i+1) = g + p g = c * r - b 200 continue c w(l) = w(l) - p rv1(l) = g rv1(m) = 0.0d0 go to 105 c .......... recover from underflow .......... 210 w(i+1) = w(i+1) - p rv1(m) = 0.0d0 go to 105 c .......... order eigenvalues .......... 215 if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. w(i-1)) go to 270 w(i) = w(i-1) ind(i) = ind(i-1) 230 continue c 250 i = 1 270 w(i) = p ind(i) = tag 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end