/* WR-276 */
/* TRIDIAGONALIZATION OF A SYMMETRIC BAND MATRIX */
/* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation",
Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 276-278. */
/* Input:
n = order of the matrix A (and of the matrix V).
m = band width of the matrix A (= number of diagonals above the main
diagonal).
matv = boolean variable to decide whether the transformation matrix V is
to be computed (matv = TRUE) or not (matv = FALSE).
a = array a[1:n, 0:m] are the elements of the given band matrix A arranged
in a rectangular array [see text].
The elements a[i,k] with i+k > n are irrelevant.
Output:
a = the elements of the resulting tridiagonal matrix J. The diagonal
elements of J are given by a[i,0] (i = 1, 2, ..., n), the elements of the
superdiagonal are given by a[i,1] (i = 1, 2, ..., n-1). The other elements
a[i,2], a[i,3], ..., a[i,m] are meaningless.
v = if matv = TRUE, the array v[1:n, 1:n] contains the elements of the
orthogonal transformation matrix V. This is a full square matrix.
If matv = FALSE, no values are assigned to v[i,k]. */
/* Transforms a real and symmetric band matrix A of order n
and band width m (i.e. with a[i,k] = 0 for all i and k with
abs(i-k) > m) to tridiagonal form by an approprite finite
sequence of Jacobi rotations. The elements of A in and above
the main diagonal are supposed as a rectangular array (array
a[1:n, 0:m]), where a[i,j] denotes the element in the i-th
row and (i+j)-th column of the usual representations of a
matrix. During the transformation the property of the band
matrix is maintained. The method yields the tridiagonal
matrix J = (Vt)AV (where Vt is the transpose of V), the
elements of which are a(i,0) (i = 1(1)n) and a[i,1] (i = 1(1)n-1).
If the parameter matv is given the value TRUE, the orthogonal
matrix V is delivered as array v[1:n, 1:n], otherwise not. */
BANDRD:PROC(N,M,MATV,A,V) OPTIONS (REORDER);
DCL (N,M,J,K,LBA1,LBA2,UBA1,UBA2,MAXR,R,UGL,L,MAXL) BIN FIXED(31,0);
DCL (A(*,*),V(*,*),B,S,C,C2,S2,CS,U,U1,G) BIN FLOAT(21);
DCL MATV BIT(1);
DCL R1 BIN FIXED(31,0);
LBA1=LBOUND(A,1); UBA1=LBA1+N-1;
LBA2=LBOUND(A,2); UBA2=LBA2+M;
IF MATV THEN DO; V=0E0;
DO K=LBA1 TO UBA1;
V(K,K)=1E0;
END; END;
DO K=LBA1 TO UBA1-2;
IF UBA1-K