1      SUBROUTINE BANDR(NM,N,MB,A,D,E,E2,MATZ,Z)
2C***BEGIN PROLOGUE  BANDR
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C1B1
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Reduces real symmetric band matrix to symmetric
9C            tridiagonal matrix and, optionally, accumulates
10C            orthogonal similarity transformations.
11C***DESCRIPTION
12C
13C     This subroutine is a translation of the ALGOL procedure BANDRD,
14C     NUM. MATH. 12, 231-241(1968) by Schwarz.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-lINEAR ALGEBRA, 273-283(1971).
16C
17C     This subroutine reduces a REAL SYMMETRIC BAND matrix
18C     to a symmetric tridiagonal matrix using and optionally
19C     accumulating orthogonal similarity transformations.
20C
21C     On INPUT
22C
23C        NM must be set to the row dimension of two-dimensional
24C          array parameters as declared in the calling program
25C          dimension statement.
26C
27C        N is the order of the matrix.
28C
29C        MB is the (half) band width of the matrix, defined as the
30C          number of adjacent diagonals, including the principal
31C          diagonal, required to specify the non-zero portion of the
32C          lower triangle of the matrix.
33C
34C        A contains the lower triangle of the symmetric band input
35C          matrix stored as an N by MB array.  Its lowest subdiagonal
36C          is stored in the last N+1-MB positions of the first column,
37C          its next subdiagonal in the last N+2-MB positions of the
38C          second column, further subdiagonals similarly, and finally
39C          its principal diagonal in the N positions of the last column.
40C          Contents of storages not part of the matrix are arbitrary.
41C
42C        MATZ should be set to .TRUE. if the transformation matrix is
43C          to be accumulated, and to .FALSE. otherwise.
44C
45C     On OUTPUT
46C
47C        A has been destroyed, except for its last two columns which
48C          contain a copy of the tridiagonal matrix.
49C
50C        D contains the diagonal elements of the tridiagonal matrix.
51C
52C        E contains the subdiagonal elements of the tridiagonal
53C          matrix in its last N-1 positions.  E(1) is set to zero.
54C
55C        E2 contains the squares of the corresponding elements of E.
56C          E2 may coincide with E if the squares are not needed.
57C
58C        Z contains the orthogonal transformation matrix produced in
59C          the reduction if MATZ has been set to .TRUE.  Otherwise, Z
60C          is not referenced.
61C
62C     Questions and comments should be directed to B. S. Garbow,
63C     Applied Mathematics Division, ARGONNE NATIONAL LABORATORY
64C     ------------------------------------------------------------------
65C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
66C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
67C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
68C                 1976.
69C***ROUTINES CALLED  (NONE)
70C***END PROLOGUE  BANDR
71C
72      INTEGER J,K,L,N,R,I1,I2,J1,J2,KR,MB,MR,M1,NM,N2,R1,UGL,MAXL,MAXR
73      REAL A(NM,MB),D(N),E(N),E2(N),Z(NM,N)
74      REAL G,U,B1,B2,C2,F1,F2,S2,DMIN,DMINRT
75      LOGICAL MATZ
76C
77C***FIRST EXECUTABLE STATEMENT  BANDR
78      DMIN = 2.0E0**(-64)
79      DMINRT = 2.0E0**(-32)
80C     .......... INITIALIZE DIAGONAL SCALING MATRIX ..........
81      DO 30 J = 1, N
82   30 D(J) = 1.0E0
83C
84      IF (.NOT. MATZ) GO TO 60
85C
86      DO 50 J = 1, N
87C
88         DO 40 K = 1, N
89   40    Z(J,K) = 0.0E0
90C
91         Z(J,J) = 1.0E0
92   50 CONTINUE
93C
94   60 M1 = MB - 1
95      IF (M1 - 1) 900, 800, 70
96   70 N2 = N - 2
97C
98      DO 700 K = 1, N2
99         MAXR = MIN0(M1,N-K)
100C     .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- ..........
101         DO 600 R1 = 2, MAXR
102            R = MAXR + 2 - R1
103            KR = K + R
104            MR = MB - R
105            G = A(KR,MR)
106            A(KR-1,1) = A(KR-1,MR+1)
107            UGL = K
108C
109            DO 500 J = KR, N, M1
110               J1 = J - 1
111               J2 = J1 - 1
112               IF (G .EQ. 0.0E0) GO TO 600
113               B1 = A(J1,1) / G
114               B2 = B1 * D(J1) / D(J)
115               S2 = 1.0E0 / (1.0E0 + B1 * B2)
116               IF (S2 .GE. 0.5E0 ) GO TO 450
117               B1 = G / A(J1,1)
118               B2 = B1 * D(J) / D(J1)
119               C2 = 1.0E0 - S2
120               D(J1) = C2 * D(J1)
121               D(J) = C2 * D(J)
122               F1 = 2.0E0 * A(J,M1)
123               F2 = B1 * A(J1,MB)
124               A(J,M1) = -B2 * (B1 * A(J,M1) - A(J,MB)) - F2 + A(J,M1)
125               A(J1,MB) = B2 * (B2 * A(J,MB) + F1) + A(J1,MB)
126               A(J,MB) = B1 * (F2 - F1) + A(J,MB)
127C
128               DO 200 L = UGL, J2
129                  I2 = MB - J + L
130                  U = A(J1,I2+1) + B2 * A(J,I2)
131                  A(J,I2) = -B1 * A(J1,I2+1) + A(J,I2)
132                  A(J1,I2+1) = U
133  200          CONTINUE
134C
135               UGL = J
136               A(J1,1) = A(J1,1) + B2 * G
137               IF (J .EQ. N) GO TO 350
138               MAXL = MIN0(M1,N-J1)
139C
140               DO 300 L = 2, MAXL
141                  I1 = J1 + L
142                  I2 = MB - L
143                  U = A(I1,I2) + B2 * A(I1,I2+1)
144                  A(I1,I2+1) = -B1 * A(I1,I2) + A(I1,I2+1)
145                  A(I1,I2) = U
146  300          CONTINUE
147C
148               I1 = J + M1
149               IF (I1 .GT. N) GO TO 350
150               G = B2 * A(I1,1)
151  350          IF (.NOT. MATZ) GO TO 500
152C
153               DO 400 L = 1, N
154                  U = Z(L,J1) + B2 * Z(L,J)
155                  Z(L,J) = -B1 * Z(L,J1) + Z(L,J)
156                  Z(L,J1) = U
157  400          CONTINUE
158C
159               GO TO 500
160C
161  450          U = D(J1)
162               D(J1) = S2 * D(J)
163               D(J) = S2 * U
164               F1 = 2.0E0 * A(J,M1)
165               F2 = B1 * A(J,MB)
166               U = B1 * (F2 - F1) + A(J1,MB)
167               A(J,M1) = B2 * (B1 * A(J,M1) - A(J1,MB)) + F2 - A(J,M1)
168               A(J1,MB) = B2 * (B2 * A(J1,MB) + F1) + A(J,MB)
169               A(J,MB) = U
170C
171               DO 460 L = UGL, J2
172                  I2 = MB - J + L
173                  U = B2 * A(J1,I2+1) + A(J,I2)
174                  A(J,I2) = -A(J1,I2+1) + B1 * A(J,I2)
175                  A(J1,I2+1) = U
176  460          CONTINUE
177C
178               UGL = J
179               A(J1,1) = B2 * A(J1,1) + G
180               IF (J .EQ. N) GO TO 480
181               MAXL = MIN0(M1,N-J1)
182C
183               DO 470 L = 2, MAXL
184                  I1 = J1 + L
185                  I2 = MB - L
186                  U = B2 * A(I1,I2) + A(I1,I2+1)
187                  A(I1,I2+1) = -A(I1,I2) + B1 * A(I1,I2+1)
188                  A(I1,I2) = U
189  470          CONTINUE
190C
191               I1 = J + M1
192               IF (I1 .GT. N) GO TO 480
193               G = A(I1,1)
194               A(I1,1) = B1 * A(I1,1)
195  480          IF (.NOT. MATZ) GO TO 500
196C
197               DO 490 L = 1, N
198                  U = B2 * Z(L,J1) + Z(L,J)
199                  Z(L,J) = -Z(L,J1) + B1 * Z(L,J)
200                  Z(L,J1) = U
201  490          CONTINUE
202C
203  500       CONTINUE
204C
205  600    CONTINUE
206C
207         IF (MOD(K,64) .NE. 0) GO TO 700
208C     .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW ..........
209         DO 650 J = K, N
210            IF (D(J) .GE. DMIN) GO TO 650
211            MAXL = MAX0(1,MB+1-J)
212C
213            DO 610 L = MAXL, M1
214  610       A(J,L) = DMINRT * A(J,L)
215C
216            IF (J .EQ. N) GO TO 630
217            MAXL = MIN0(M1,N-J)
218C
219            DO 620 L = 1, MAXL
220               I1 = J + L
221               I2 = MB - L
222               A(I1,I2) = DMINRT * A(I1,I2)
223  620       CONTINUE
224C
225  630       IF (.NOT. MATZ) GO TO 645
226C
227            DO 640 L = 1, N
228  640       Z(L,J) = DMINRT * Z(L,J)
229C
230  645       A(J,MB) = DMIN * A(J,MB)
231            D(J) = D(J) / DMIN
232  650    CONTINUE
233C
234  700 CONTINUE
235C     .......... FORM SQUARE ROOT OF SCALING MATRIX ..........
236  800 DO 810 J = 2, N
237  810 E(J) = SQRT(D(J))
238C
239      IF (.NOT. MATZ) GO TO 840
240C
241      DO 830 J = 1, N
242C
243         DO 820 K = 2, N
244  820    Z(J,K) = E(K) * Z(J,K)
245C
246  830 CONTINUE
247C
248  840 U = 1.0E0
249C
250      DO 850 J = 2, N
251         A(J,M1) = U * E(J) * A(J,M1)
252         U = E(J)
253         E2(J) = A(J,M1) ** 2
254         A(J,MB) = D(J) * A(J,MB)
255         D(J) = A(J,MB)
256         E(J) = A(J,M1)
257  850 CONTINUE
258C
259      D(1) = A(1,MB)
260      E(1) = 0.0E0
261      E2(1) = 0.0E0
262      GO TO 1001
263C
264  900 DO 950 J = 1, N
265         D(J) = A(J,MB)
266         E(J) = 0.0E0
267         E2(J) = 0.0E0
268  950 CONTINUE
269C
270 1001 RETURN
271      END
272