subroutine bandv(nm,n,mbw,a,e21,m,w,z,ierr,nv,rv,rv6) c integer i,j,k,m,n,r,ii,ij,jj,kj,mb,m1,nm,nv,ij1,its,kj1,mbw,m21, x ierr,maxj,maxk,group double precision a(nm,mbw),w(m),z(nm,m),rv(nv),rv6(n) double precision u,v,uk,xu,x0,x1,e21,eps2,eps3,eps4,norm,order, x epslon,pythag c c this subroutine finds those eigenvectors of a real symmetric c band matrix corresponding to specified eigenvalues, using inverse c iteration. the subroutine may also be used to solve systems c of linear equations with a symmetric or non-symmetric band c coefficient matrix. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c mbw is the number of columns of the array a used to store the c band matrix. if the matrix is symmetric, mbw is its (half) c band width, denoted mb and defined as the number of adjacent c diagonals, including the principal diagonal, required to c specify the non-zero portion of the lower triangle of the c matrix. if the subroutine is being used to solve systems c of linear equations and the coefficient matrix is not c symmetric, it must however have the same number of adjacent c diagonals above the main diagonal as below, and in this c case, mbw=2*mb-1. c c a contains the lower triangle of the symmetric band input c matrix stored as an n by mb array. its lowest subdiagonal c is stored in the last n+1-mb positions of the first column, c its next subdiagonal in the last n+2-mb positions of the c second column, further subdiagonals similarly, and finally c its principal diagonal in the n positions of column mb. c if the subroutine is being used to solve systems of linear c equations and the coefficient matrix is not symmetric, a is c n by 2*mb-1 instead with lower triangle as above and with c its first superdiagonal stored in the first n-1 positions of c column mb+1, its second superdiagonal in the first n-2 c positions of column mb+2, further superdiagonals similarly, c and finally its highest superdiagonal in the first n+1-mb c positions of the last column. c contents of storages not part of the matrix are arbitrary. c c e21 specifies the ordering of the eigenvalues and contains c 0.0d0 if the eigenvalues are in ascending order, or c 2.0d0 if the eigenvalues are in descending order. c if the subroutine is being used to solve systems of linear c equations, e21 should be set to 1.0d0 if the coefficient c matrix is symmetric and to -1.0d0 if not. c c m is the number of specified eigenvalues or the number of c systems of linear equations. c c w contains the m eigenvalues in ascending or descending order. c if the subroutine is being used to solve systems of linear c equations (a-w(r)*i)*x(r)=b(r), where i is the identity c matrix, w(r) should be set accordingly, for r=1,2,...,m. c c z contains the constant matrix columns (b(r),r=1,2,...,m), if c the subroutine is used to solve systems of linear equations. c c nv must be set to the dimension of the array parameter rv c as declared in the calling program dimension statement. c c on output c c a and w are unaltered. c c z contains the associated set of orthogonal eigenvectors. c any vector which fails to converge is set to zero. if the c subroutine is used to solve systems of linear equations, c z contains the solution matrix columns (x(r),r=1,2,...,m). c c ierr is set to c zero for normal return, c -r if the eigenvector corresponding to the r-th c eigenvalue fails to converge, or if the r-th c system of linear equations is nearly singular. c c rv and rv6 are temporary storage arrays. note that rv is c of dimension at least n*(2*mb-1). if the subroutine c is being used to solve systems of linear equations, the c determinant (up to sign) of a-w(m)*i is available, upon c return, as the product of the first n elements of rv. 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 if (m .eq. 0) go to 1001 mb = mbw if (e21 .lt. 0.0d0) mb = (mbw + 1) / 2 m1 = mb - 1 m21 = m1 + mb order = 1.0d0 - dabs(e21) c .......... find vectors by inverse iteration .......... do 920 r = 1, m its = 1 x1 = w(r) if (r .ne. 1) go to 100 c .......... compute norm of matrix .......... norm = 0.0d0 c do 60 j = 1, mb jj = mb + 1 - j kj = jj + m1 ij = 1 v = 0.0d0 c do 40 i = jj, n v = v + dabs(a(i,j)) if (e21 .ge. 0.0d0) go to 40 v = v + dabs(a(ij,kj)) ij = ij + 1 40 continue c norm = dmax1(norm,v) 60 continue c if (e21 .lt. 0.0d0) norm = 0.5d0 * norm c .......... eps2 is the criterion for grouping, c eps3 replaces zero pivots and equal c roots are modified by eps3, c eps4 is taken very small to avoid overflow .......... if (norm .eq. 0.0d0) norm = 1.0d0 eps2 = 1.0d-3 * norm * dabs(order) eps3 = epslon(norm) uk = n uk = dsqrt(uk) eps4 = uk * eps3 80 group = 0 go to 120 c .......... look for close or coincident roots .......... 100 if (dabs(x1-x0) .ge. eps2) go to 80 group = group + 1 if (order * (x1 - x0) .le. 0.0d0) x1 = x0 + order * eps3 c .......... expand matrix, subtract eigenvalue, c and initialize vector .......... 120 do 200 i = 1, n ij = i + min0(0,i-m1) * n kj = ij + mb * n ij1 = kj + m1 * n if (m1 .eq. 0) go to 180 c do 150 j = 1, m1 if (ij .gt. m1) go to 125 if (ij .gt. 0) go to 130 rv(ij1) = 0.0d0 ij1 = ij1 + n go to 130 125 rv(ij) = a(i,j) 130 ij = ij + n ii = i + j if (ii .gt. n) go to 150 jj = mb - j if (e21 .ge. 0.0d0) go to 140 ii = i jj = mb + j 140 rv(kj) = a(ii,jj) kj = kj + n 150 continue c 180 rv(ij) = a(i,mb) - x1 rv6(i) = eps4 if (order .eq. 0.0d0) rv6(i) = z(i,r) 200 continue c if (m1 .eq. 0) go to 600 c .......... elimination with interchanges .......... do 580 i = 1, n ii = i + 1 maxk = min0(i+m1-1,n) maxj = min0(n-i,m21-2) * n c do 360 k = i, maxk kj1 = k j = kj1 + n jj = j + maxj c do 340 kj = j, jj, n rv(kj1) = rv(kj) kj1 = kj 340 continue c rv(kj1) = 0.0d0 360 continue c if (i .eq. n) go to 580 u = 0.0d0 maxk = min0(i+m1,n) maxj = min0(n-ii,m21-2) * n c do 450 j = i, maxk if (dabs(rv(j)) .lt. dabs(u)) go to 450 u = rv(j) k = j 450 continue c j = i + n jj = j + maxj if (k .eq. i) go to 520 kj = k c do 500 ij = i, jj, n v = rv(ij) rv(ij) = rv(kj) rv(kj) = v kj = kj + n 500 continue c if (order .ne. 0.0d0) go to 520 v = rv6(i) rv6(i) = rv6(k) rv6(k) = v 520 if (u .eq. 0.0d0) go to 580 c do 560 k = ii, maxk v = rv(k) / u kj = k c do 540 ij = j, jj, n kj = kj + n rv(kj) = rv(kj) - v * rv(ij) 540 continue c if (order .eq. 0.0d0) rv6(k) = rv6(k) - v * rv6(i) 560 continue c 580 continue c .......... back substitution c for i=n step -1 until 1 do -- .......... 600 do 630 ii = 1, n i = n + 1 - ii maxj = min0(ii,m21) if (maxj .eq. 1) go to 620 ij1 = i j = ij1 + n jj = j + (maxj - 2) * n c do 610 ij = j, jj, n ij1 = ij1 + 1 rv6(i) = rv6(i) - rv(ij) * rv6(ij1) 610 continue c 620 v = rv(i) if (dabs(v) .ge. eps3) go to 625 c .......... set error -- nearly singular linear system .......... if (order .eq. 0.0d0) ierr = -r v = dsign(eps3,v) 625 rv6(i) = rv6(i) / v 630 continue c xu = 1.0d0 if (order .eq. 0.0d0) go to 870 c .......... orthogonalize with respect to previous c members of group .......... if (group .eq. 0) go to 700 c do 680 jj = 1, group j = r - group - 1 + jj xu = 0.0d0 c do 640 i = 1, n 640 xu = xu + rv6(i) * z(i,j) c do 660 i = 1, n 660 rv6(i) = rv6(i) - xu * z(i,j) c 680 continue c 700 norm = 0.0d0 c do 720 i = 1, n 720 norm = norm + dabs(rv6(i)) c if (norm .ge. 0.1d0) go to 840 c .......... in-line procedure for choosing c a new starting vector .......... if (its .ge. n) go to 830 its = its + 1 xu = eps4 / (uk + 1.0d0) rv6(1) = eps4 c do 760 i = 2, n 760 rv6(i) = xu c rv6(its) = rv6(its) - eps4 * uk go to 600 c .......... set error -- non-converged eigenvector .......... 830 ierr = -r xu = 0.0d0 go to 870 c .......... normalize so that sum of squares is c 1 and expand to full order .......... 840 u = 0.0d0 c do 860 i = 1, n 860 u = pythag(u,rv6(i)) c xu = 1.0d0 / u c 870 do 900 i = 1, n 900 z(i,r) = rv6(i) * xu c x0 = x1 920 continue c 1001 return end