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