2
SUBROUTINE BQR (NM, N, MB, A, T, R, IERR, NV, RV)
4
C***PURPOSE Compute some of the eigenvalues of a real symmetric
5
C matrix using the QR method with shifts of origin.
6
C***LIBRARY SLATEC (EISPACK)
8
C***TYPE SINGLE PRECISION (BQR-S)
9
C***KEYWORDS EIGENVALUES, EISPACK
10
C***AUTHOR Smith, B. T., et al.
13
C This subroutine is a translation of the ALGOL procedure BQR,
14
C NUM. MATH. 16, 85-92(1970) by Martin, Reinsch, and Wilkinson.
15
C HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971).
17
C This subroutine finds the eigenvalue of smallest (usually)
18
C magnitude of a REAL SYMMETRIC BAND matrix using the
19
C QR algorithm with shifts of origin. Consecutive calls
20
C can be made to find further eigenvalues.
24
C NM must be set to the row dimension of the two-dimensional
25
C array parameter, A, as declared in the calling program
26
C dimension statement. NM is an INTEGER variable.
28
C N is the order of the matrix A. N is an INTEGER variable.
29
C N must be less than or equal to NM.
31
C MB is the (half) band width of the matrix, defined as the
32
C number of adjacent diagonals, including the principal
33
C diagonal, required to specify the non-zero portion of the
34
C lower triangle of the matrix. MB is an INTEGER variable.
35
C MB must be less than or equal to N on first call.
37
C A contains the lower triangle of the symmetric band input
38
C matrix stored as an N by MB array. Its lowest subdiagonal
39
C is stored in the last N+1-MB positions of the first column,
40
C its next subdiagonal in the last N+2-MB positions of the
41
C second column, further subdiagonals similarly, and finally
42
C its principal diagonal in the N positions of the last column.
43
C Contents of storages not part of the matrix are arbitrary.
44
C On a subsequent call, its output contents from the previous
45
C call should be passed. A is a two-dimensional REAL array,
46
C dimensioned A(NM,MB).
48
C T specifies the shift (of eigenvalues) applied to the diagonal
49
C of A in forming the input matrix. What is actually determined
50
C is the eigenvalue of A+TI (I is the identity matrix) nearest
51
C to T. On a subsequent call, the output value of T from the
52
C previous call should be passed if the next nearest eigenvalue
53
C is sought. T is a REAL variable.
55
C R should be specified as zero on the first call, and as its
56
C output value from the previous call on a subsequent call.
57
C It is used to determine when the last row and column of
58
C the transformed band matrix can be regarded as negligible.
59
C R is a REAL variable.
61
C NV must be set to the dimension of the array parameter RV
62
C as declared in the calling program dimension statement.
63
C NV is an INTEGER variable.
67
C A contains the transformed band matrix. The matrix A+TI
68
C derived from the output parameters is similar to the
69
C input A+TI to within rounding errors. Its last row and
70
C column are null (if IERR is zero).
72
C T contains the computed eigenvalue of A+TI (if IERR is zero),
73
C where I is the identity matrix.
75
C R contains the maximum of its input value and the norm of the
76
C last column of the input matrix A.
78
C IERR is an INTEGER flag set to
79
C Zero for normal return,
80
C J if the J-th eigenvalue has not been
81
C determined after a total of 30 iterations.
83
C RV is a one-dimensional REAL array of dimension NV which is
84
C at least (2*MB**2+4*MB-3), used for temporary storage. The
85
C first (3*MB-2) locations correspond to the ALGOL array B,
86
C the next (2*MB-1) locations correspond to the ALGOL array H,
87
C and the final (2*MB**2-MB) locations correspond to the MB
88
C by (2*MB-1) ALGOL array U.
90
C NOTE. For a subsequent call, N should be replaced by N-1, but
91
C MB should not be altered even when it exceeds the current N.
93
C Calls PYTHAG(A,B) for SQRT(A**2 + B**2).
95
C Questions and comments should be directed to B. S. Garbow,
96
C Applied Mathematics Division, ARGONNE NATIONAL LABORATORY
97
C ------------------------------------------------------------------
99
C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
100
C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
101
C system Routines - EISPACK Guide, Springer-Verlag,
103
C***ROUTINES CALLED PYTHAG
104
C***REVISION HISTORY (YYMMDD)
105
C 760101 DATE WRITTEN
106
C 890531 Changed all specific intrinsics to generic. (WRB)
107
C 890831 Modified array declarations. (WRB)
108
C 890831 REVISION DATE from Version 3.2
109
C 891214 Prologue converted to Version 4.0 format. (BAB)
110
C 920501 Reformatted the REFERENCES section. (WRB)
113
INTEGER I,J,K,L,M,N,II,IK,JK,JM,KJ,KK,KM,LL,MB,MK,MN,MZ
114
INTEGER M1,M2,M3,M4,NI,NM,NV,ITS,KJ1,M21,M31,IERR,IMULT
116
REAL F,G,Q,R,S,T,SCALE
119
C***FIRST EXECUTABLE STATEMENT BQR
131
C .......... TEST FOR CONVERGENCE ..........
133
IF (M .EQ. 0) GO TO 360
141
IF (ITS .EQ. 0 .AND. F .GT. R) R = F
142
IF (R + F .LE. R) GO TO 360
143
IF (ITS .EQ. 30) GO TO 1000
145
C .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR ..........
146
IF (F .GT. 0.25E0 * R .AND. ITS .LT. 5) GO TO 90
148
IF (F .EQ. 0.0E0) GO TO 70
149
Q = (A(N-1,MB) - G) / (2.0E0 * F)
151
G = G - F / (Q + SIGN(S,Q))
155
80 A(I,MB) = A(I,MB) - G
157
90 DO 100 K = M31, M4
163
IF (NI .LT. 0) GO TO 230
164
C .......... FORM COLUMN OF SHIFTED MATRIX A-G*I ..........
177
IF (LL .EQ. 0) GO TO 135
185
C .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
188
C .......... MULTIPLICATION PROCEDURE ..........
194
IF (RV(JM) .EQ. 0.0E0) GO TO 170
200
F = F + RV(KJ) * RV(JK)
209
RV(JK) = RV(JK) - RV(KJ) * F
215
IF (IMULT .NE. 0) GO TO 280
216
C .......... HOUSEHOLDER REFLECTION ..........
223
180 SCALE = SCALE + ABS(RV(K))
225
IF (SCALE .EQ. 0.0E0) GO TO 210
228
190 S = S + (RV(K)/SCALE)**2
230
S = SCALE * SCALE * S
234
KJ = M4 + M2 * M1 + 1
242
C .......... SAVE COLUMN OF TRIANGULAR FACTOR R ..........
249
230 L = MAX(1,M1+1-I)
250
IF (I .LE. 0) GO TO 300
251
C .......... PERFORM ADDITIONAL STEPS ..........
256
C .......... GET ROW OF TRIANGULAR FACTOR R ..........
264
C .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
268
C .......... STORE COLUMN OF NEW A MATRIX ..........
273
C .......... UPDATE HOUSEHOLDER REFLECTIONS ..........
274
300 IF (L .GT. 1) L = L - 1
290
C .......... CONVERGENCE ..........
294
380 A(I,MB) = A(I,MB) - G
302
C .......... SET ERROR -- NO CONVERGENCE TO
303
C EIGENVALUE AFTER 30 ITERATIONS ..........
added
removed
deps/openlibm/slatec/bqr.f
