relibc_openlibm/slatec/bqr.f

*DECK BQR
      SUBROUTINE BQR (NM, N, MB, A, T, R, IERR, NV, RV)
C***BEGIN PROLOGUE  BQR
C***PURPOSE  Compute some of the eigenvalues of a real symmetric
C            matrix using the QR method with shifts of origin.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4A6
C***TYPE      SINGLE PRECISION (BQR-S)
C***KEYWORDS  EIGENVALUES, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine is a translation of the ALGOL procedure BQR,
C     NUM. MATH. 16, 85-92(1970) by Martin, Reinsch, and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971).
C
C     This subroutine finds the eigenvalue of smallest (usually)
C     magnitude of a REAL SYMMETRIC BAND matrix using the
C     QR algorithm with shifts of origin.  Consecutive calls
C     can be made to find further eigenvalues.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameter, A, as declared in the calling program
C          dimension statement.  NM is an INTEGER variable.
C
C        N is the order of the matrix A.  N is an INTEGER variable.
C          N must be less than or equal to NM.
C
C        MB is the (half) band width of the matrix, defined as the
C          number of adjacent diagonals, including the principal
C          diagonal, required to specify the non-zero portion of the
C          lower triangle of the matrix.  MB is an INTEGER variable.
C          MB must be less than or equal to N on first call.
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 the last column.
C          Contents of storages not part of the matrix are arbitrary.
C          On a subsequent call, its output contents from the previous
C          call should be passed.  A is a two-dimensional REAL array,
C          dimensioned A(NM,MB).
C
C        T specifies the shift (of eigenvalues) applied to the diagonal
C          of A in forming the input matrix. What is actually determined
C          is the eigenvalue of A+TI (I is the identity matrix) nearest
C          to T.  On a subsequent call, the output value of T from the
C          previous call should be passed if the next nearest eigenvalue
C          is sought.  T is a REAL variable.
C
C        R should be specified as zero on the first call, and as its
C          output value from the previous call on a subsequent call.
C          It is used to determine when the last row and column of
C          the transformed band matrix can be regarded as negligible.
C          R is a REAL variable.
C
C        NV must be set to the dimension of the array parameter RV
C          as declared in the calling program dimension statement.
C          NV is an INTEGER variable.
C
C     On OUTPUT
C
C        A contains the transformed band matrix.  The matrix A+TI
C          derived from the output parameters is similar to the
C          input A+TI to within rounding errors.  Its last row and
C          column are null (if IERR is zero).
C
C        T contains the computed eigenvalue of A+TI (if IERR is zero),
C          where I is the identity matrix.
C
C        R contains the maximum of its input value and the norm of the
C          last column of the input matrix A.
C
C        IERR is an INTEGER flag set to
C          Zero       for normal return,
C          J          if the J-th eigenvalue has not been
C                     determined after a total of 30 iterations.
C
C        RV is a one-dimensional REAL array of dimension NV which is
C          at least (2*MB**2+4*MB-3), used for temporary storage.  The
C          first (3*MB-2) locations correspond to the ALGOL array B,
C          the next (2*MB-1) locations correspond to the ALGOL array H,
C          and the final (2*MB**2-MB) locations correspond to the MB
C          by (2*MB-1) ALGOL array U.
C
C     NOTE. For a subsequent call, N should be replaced by N-1, but
C     MB should not be altered even when it exceeds the current N.
C
C     Calls PYTHAG(A,B) for SQRT(A**2 + B**2).
C
C     Questions and comments should be directed to B. S. Garbow,
C     Applied Mathematics Division, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  PYTHAG
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  BQR
C
      INTEGER I,J,K,L,M,N,II,IK,JK,JM,KJ,KK,KM,LL,MB,MK,MN,MZ
      INTEGER M1,M2,M3,M4,NI,NM,NV,ITS,KJ1,M21,M31,IERR,IMULT
      REAL A(NM,*),RV(*)
      REAL F,G,Q,R,S,T,SCALE
      REAL PYTHAG
C
C***FIRST EXECUTABLE STATEMENT  BQR
      IERR = 0
      M1 = MIN(MB,N)
      M = M1 - 1
      M2 = M + M
      M21 = M2 + 1
      M3 = M21 + M
      M31 = M3 + 1
      M4 = M31 + M2
      MN = M + N
      MZ = MB - M1
      ITS = 0
C     .......... TEST FOR CONVERGENCE ..........
   40 G = A(N,MB)
      IF (M .EQ. 0) GO TO 360
      F = 0.0E0
C
      DO 50 K = 1, M
         MK = K + MZ
         F = F + ABS(A(N,MK))
   50 CONTINUE
C
      IF (ITS .EQ. 0 .AND. F .GT. R) R = F
      IF (R + F .LE. R) GO TO 360
      IF (ITS .EQ. 30) GO TO 1000
      ITS = ITS + 1
C     .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR ..........
      IF (F .GT. 0.25E0 * R .AND. ITS .LT. 5) GO TO 90
      F = A(N,MB-1)
      IF (F .EQ. 0.0E0) GO TO 70
      Q = (A(N-1,MB) - G) / (2.0E0 * F)
      S = PYTHAG(Q,1.0E0)
      G = G - F / (Q + SIGN(S,Q))
   70 T = T + G
C
      DO 80 I = 1, N
   80 A(I,MB) = A(I,MB) - G
C
   90 DO 100 K = M31, M4
  100 RV(K) = 0.0E0
C
      DO 350 II = 1, MN
         I = II - M
         NI = N - II
         IF (NI .LT. 0) GO TO 230
C     .......... FORM COLUMN OF SHIFTED MATRIX A-G*I ..........
         L = MAX(1,2-I)
C
         DO 110 K = 1, M3
  110    RV(K) = 0.0E0
C
         DO 120 K = L, M1
            KM = K + M
            MK = K + MZ
            RV(KM) = A(II,MK)
  120    CONTINUE
C
         LL = MIN(M,NI)
         IF (LL .EQ. 0) GO TO 135
C
         DO 130 K = 1, LL
            KM = K + M21
            IK = II + K
            MK = MB - K
            RV(KM) = A(IK,MK)
  130    CONTINUE
C     .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
  135    LL = M2
         IMULT = 0
C     .......... MULTIPLICATION PROCEDURE ..........
  140    KJ = M4 - M1
C
         DO 170 J = 1, LL
            KJ = KJ + M1
            JM = J + M3
            IF (RV(JM) .EQ. 0.0E0) GO TO 170
            F = 0.0E0
C
            DO 150 K = 1, M1
               KJ = KJ + 1
               JK = J + K - 1
               F = F + RV(KJ) * RV(JK)
  150       CONTINUE
C
            F = F / RV(JM)
            KJ = KJ - M1
C
            DO 160 K = 1, M1
               KJ = KJ + 1
               JK = J + K - 1
               RV(JK) = RV(JK) - RV(KJ) * F
  160       CONTINUE
C
            KJ = KJ - M1
  170    CONTINUE
C
         IF (IMULT .NE. 0) GO TO 280
C     .......... HOUSEHOLDER REFLECTION ..........
         F = RV(M21)
         S = 0.0E0
         RV(M4) = 0.0E0
         SCALE = 0.0E0
C
         DO 180 K = M21, M3
  180    SCALE = SCALE + ABS(RV(K))
C
         IF (SCALE .EQ. 0.0E0) GO TO 210
C
         DO 190 K = M21, M3
  190    S = S + (RV(K)/SCALE)**2
C
         S = SCALE * SCALE * S
         G = -SIGN(SQRT(S),F)
         RV(M21) = G
         RV(M4) = S - F * G
         KJ = M4 + M2 * M1 + 1
         RV(KJ) = F - G
C
         DO 200 K = 2, M1
            KJ = KJ + 1
            KM = K + M2
            RV(KJ) = RV(KM)
  200    CONTINUE
C     .......... SAVE COLUMN OF TRIANGULAR FACTOR R ..........
  210    DO 220 K = L, M1
            KM = K + M
            MK = K + MZ
            A(II,MK) = RV(KM)
  220    CONTINUE
C
  230    L = MAX(1,M1+1-I)
         IF (I .LE. 0) GO TO 300
C     .......... PERFORM ADDITIONAL STEPS ..........
         DO 240 K = 1, M21
  240    RV(K) = 0.0E0
C
         LL = MIN(M1,NI+M1)
C     .......... GET ROW OF TRIANGULAR FACTOR R ..........
         DO 250 KK = 1, LL
            K = KK - 1
            KM = K + M1
            IK = I + K
            MK = MB - K
            RV(KM) = A(IK,MK)
  250    CONTINUE
C     .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
         LL = M1
         IMULT = 1
         GO TO 140
C     .......... STORE COLUMN OF NEW A MATRIX ..........
  280    DO 290 K = L, M1
            MK = K + MZ
            A(I,MK) = RV(K)
  290    CONTINUE
C     .......... UPDATE HOUSEHOLDER REFLECTIONS ..........
  300    IF (L .GT. 1) L = L - 1
         KJ1 = M4 + L * M1
C
         DO 320 J = L, M2
            JM = J + M3
            RV(JM) = RV(JM+1)
C
            DO 320 K = 1, M1
               KJ1 = KJ1 + 1
               KJ = KJ1 - M1
               RV(KJ) = RV(KJ1)
  320    CONTINUE
C
  350 CONTINUE
C
      GO TO 40
C     .......... CONVERGENCE ..........
  360 T = T + G
C
      DO 380 I = 1, N
  380 A(I,MB) = A(I,MB) - G
C
      DO 400 K = 1, M1
         MK = K + MZ
         A(N,MK) = 0.0E0
  400 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- NO CONVERGENCE TO
C                EIGENVALUE AFTER 30 ITERATIONS ..........
 1000 IERR = N
 1001 RETURN
      END