%PROCESS MARGINS(1,100);
/* SOLUTION OF SYMMETRIC & UNSYMMETRIC BAND EQUATIONS. */
/* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation",
Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 76-86. */
/* Computes a[i,j], the n x r array giving the r eigenvectors of the
symmetric band matrix, A, corresponding to r approximate eigenvalues
given in the array lambda[i], i = 1(1)r. Each eigenvector is computed
by inverse iteration using BANDED2 and BANDSOL2. In the first iteration
for each eigenvector, a back-substitution only is performed using
BANDSOL2 with e = 1 and a right-hand side consisting of unity elements.
The integer array l[i], i = 1(1)r, gives the number of iterations taken
to determine each vector. A limit, la, is set to the maximum number
of iterations permitted for any one eigenvector. The integer array
c[i], i = 1(1)r, is such that c[i] is the number of eigenvalues of A
greater than lambda[i]. From the r eigenvectors, r improved eigenvalues
are determined using the Rayleigh quotients and are overwritten
on the original eigenvalues in the array lambda[i]. The procedure
INNERPROD is used in the calculation of the Rayleigh quotients. */
/* p. 77 */
/* The band matrix, A, of order n with m1 sub-diagonal
elements and m2 super-diagonal elements in a typical row is
stored as an n x (m1+m2-1) array, a[1:n, -m1:m2].
The matrix (A - lambda x I) is factorised into the product
of a lower-triangular matrix and an upper-triangular matrix
using binary row interchanges. The lower-triangle is
stored as the n x m1 array m[i,j] and the upper triangle is
overwritten on A. Details of the binary row interchanges
are stored in the Boolean n x m1 array, int[i,j]. When A
is symmetric, c gives the number of eigenvalues greater
than lambda. The determinant, d1 x 2^d2, is also computed.
The parameter e influences the path to be taken when a
zero pivot arises. The procedure fails when e = 0 and
(A - lambda x I), modified by the rounding errors, is singular. */
/* p. 79 */
/* When e = 0, this procedure solves (A - lambda x I)x = b,
where A is a band matrix of order n and (A - lambda x I)
has previously been factorized using BANDET2. Here, b is
an n x r matrix consisting of r right-hand sides. Each right-
hand side requires one forward substitution and one backward
substitution. When e = 1, a backward substitution only is used.
This variant should be used as the first step in inverse
iteration for the eigenvector corresponding to lambda. */
/* p. 81 */
/* Computes zr[i,j] and zl[i,j], two n x r arrays giving the r right-hand
eigenvectors and the r left-hand eienvectors of the band matrix, A,
corresponding to the r approximate eigenvalues given in the array
lambda[i], i = 1(1)r. The right-hand and left-hand eigenvectors
corresponding to lambda[i] are both computed by inverse iteration using
the decomposition of (A - lambda x I) given by BANDET1. In the first
iteration for each right-hand eigenvector, a backward substitution only is
performed using BANSOL1 with e = 1 and a right-hand side consisting
of unity elements. Inverse iteration for the left-hand eigenvectors
forms part of the main body of this procedure. The integer arrays l[i]
and c[i] give the number of iterations taken to determine each right-
hand and left-hand eigenvector respectively. A limit, la, is set to the
maximum number of iterations permitted for any one eigenvector.
From the r eigenvectors, r improved eigenvalues are determined using
the generalized Rayleigh quotients and are overwritten on the original
eigenvalues in the array lambda[i]. The procedure INNERPROD is used
in the calculation of the Rayleigh quotients. */
/* p. 82 */
/* The band matrix, A, of order n with m1 sub-diagonal elements
and m2 super-diagonal elements in a typical row is
stored as an n x (m1+m2+1) array, a[a:n, -m1:m2].
The matrix (A - lambda x I) is factorized into the product
of a lower-triangular matrix and an upper-triangular matrix
using partial pivoting. The lower triangle is stored as
the n x m1 array m[i,j] and the upper-triangle is overwritten
on A. Details of the interchanges are stored in the
array int[i], i = 1(1)n. The determinant, d1 x 2^d2, is also
computed. The parameter e influences the path to be taken
when a zero pivot arises. The procedure fails when e = 0
and (A - lambda x I), modified by the rounding errors, is
singular. */
/* p. 83 */
/* When e = 0 this procedure solves (A - lambda x I)x = b,
where A is a band matrix of order n amd (A - lambda x I)
has previously been factorized using BANDET1. Here, b is an
n x r matrix consisting of r right-hand sides. Each right-
hand side requires one forward substitution and one backward
substitution. When e = 1, a backward substitution only is used.
This variant should be used as the first step in inverse
iteration for the eigenvector corresponding to lambda. */