relibc_openlibm/slatec/bandv.f

*DECK BANDV
      SUBROUTINE BANDV (NM, N, MBW, A, E21, M, W, Z, IERR, NV, RV, RV6)
C***BEGIN PROLOGUE  BANDV
C***PURPOSE  Form the eigenvectors of a real symmetric band matrix
C            associated with a set of ordered approximate eigenvalues
C            by inverse iteration.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4C3
C***TYPE      SINGLE PRECISION (BANDV-S)
C***KEYWORDS  EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine finds those eigenvectors of a REAL SYMMETRIC
C     BAND matrix corresponding to specified eigenvalues, using inverse
C     iteration.  The subroutine may also be used to solve systems
C     of linear equations with a symmetric or non-symmetric band
C     coefficient matrix.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameters, A and Z, as declared in the calling
C          program 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        MBW is the number of columns of the array A used to store the
C          band matrix.  If the matrix is symmetric, MBW is its (half)
C          band width, denoted MB and defined as the number of adjacent
C          diagonals, including the principal diagonal, required to
C          specify the non-zero portion of the lower triangle of the
C          matrix.  If the subroutine is being used to solve systems
C          of linear equations and the coefficient matrix is not
C          symmetric, it must however have the same number of adjacent
C          diagonals above the main diagonal as below, and in this
C          case, MBW=2*MB-1.  MBW is an INTEGER variable.  MB must not
C          be greater than N.
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 column MB.
C          If the subroutine is being used to solve systems of linear
C          equations and the coefficient matrix is not symmetric, A is
C          N by 2*MB-1 instead with lower triangle as above and with
C          its first superdiagonal stored in the first N-1 positions of
C          column MB+1, its second superdiagonal in the first N-2
C          positions of column MB+2, further superdiagonals similarly,
C          and finally its highest superdiagonal in the first N+1-MB
C          positions of the last column.  Contents of storage locations
C          not part of the matrix are arbitrary.  A is a two-dimensional
C          REAL array, dimensioned A(NM,MBW).
C
C        E21 specifies the ordering of the eigenvalues and contains
C            0.0E0 if the eigenvalues are in ascending order, or
C            2.0E0 if the eigenvalues are in descending order.
C          If the subroutine is being used to solve systems of linear
C          equations, E21 should be set to 1.0E0 if the coefficient
C          matrix is symmetric and to -1.0E0 if not.  E21 is a REAL
C          variable.
C
C        M is the number of specified eigenvalues or the number of
C          systems of linear equations.  M is an INTEGER variable.
C
C        W contains the M eigenvalues in ascending or descending order.
C          If the subroutine is being used to solve systems of linear
C          equations (A-W(J)*I)*X(J)=B(J), where I is the identity
C          matrix, W(J) should be set accordingly, for J=1,2,...,M.
C          W is a one-dimensional REAL array, dimensioned W(M).
C
C        Z contains the constant matrix columns (B(J),J=1,2,...,M), if
C          the subroutine is used to solve systems of linear equations.
C          Z is a two-dimensional REAL array, dimensioned Z(NM,M).
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 and W are unaltered.
C
C        Z contains the associated set of orthogonal eigenvectors.
C          Any vector which fails to converge is set to zero.  If the
C          subroutine is used to solve systems of linear equations,
C          Z contains the solution matrix columns (X(J),J=1,2,...,M).
C
C        IERR is an INTEGER flag set to
C          Zero       for normal return,
C          -J         if the eigenvector corresponding to the J-th
C                     eigenvalue fails to converge, or if the J-th
C                     system of linear equations is nearly singular.
C
C        RV and RV6 are temporary storage arrays.  If the subroutine
C          is being used to solve systems of linear equations, the
C          determinant (up to sign) of A-W(M)*I is available, upon
C          return, as the product of the first N elements of RV.
C          RV and RV6 are one-dimensional REAL arrays.  Note that RV
C          is dimensioned RV(NV), where NV must be at least N*(2*MB-1).
C          RV6 is dimensioned RV6(N).
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  (NONE)
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  BANDV
C
      INTEGER I,J,K,M,N,R,II,IJ,JJ,KJ,MB,M1,NM,NV,IJ1,ITS,KJ1,MBW,M21
      INTEGER IERR,MAXJ,MAXK,GROUP
      REAL A(NM,*),W(*),Z(NM,*),RV(*),RV6(*)
      REAL U,V,UK,XU,X0,X1,E21,EPS2,EPS3,EPS4,NORM,ORDER,S
C
C***FIRST EXECUTABLE STATEMENT  BANDV
      IERR = 0
      IF (M .EQ. 0) GO TO 1001
      MB = MBW
      IF (E21 .LT. 0.0E0) MB = (MBW + 1) / 2
      M1 = MB - 1
      M21 = M1 + MB
      ORDER = 1.0E0 - ABS(E21)
C     .......... FIND VECTORS BY INVERSE ITERATION ..........
      DO 920 R = 1, M
         ITS = 1
         X1 = W(R)
         IF (R .NE. 1) GO TO 100
C     .......... COMPUTE NORM OF MATRIX ..........
         NORM = 0.0E0
C
         DO 60 J = 1, MB
            JJ = MB + 1 - J
            KJ = JJ + M1
            IJ = 1
            S = 0.0E0
C
            DO 40 I = JJ, N
               S = S + ABS(A(I,J))
               IF (E21 .GE. 0.0E0) GO TO 40
               S = S + ABS(A(IJ,KJ))
               IJ = IJ + 1
   40       CONTINUE
C
            NORM = MAX(NORM,S)
   60    CONTINUE
C
         IF (E21 .LT. 0.0E0) NORM = 0.5E0 * NORM
C     .......... EPS2 IS THE CRITERION FOR GROUPING,
C                EPS3 REPLACES ZERO PIVOTS AND EQUAL
C                ROOTS ARE MODIFIED BY EPS3,
C                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ..........
         IF (NORM .EQ. 0.0E0) NORM = 1.0E0
         EPS2 = 1.0E-3 * NORM * ABS(ORDER)
         EPS3 = NORM
   70    EPS3 = 0.5E0*EPS3
         IF (NORM + EPS3 .GT. NORM) GO TO 70
         UK = SQRT(REAL(N))
         EPS3 = UK * EPS3
         EPS4 = UK * EPS3
   80    GROUP = 0
         GO TO 120
C     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
  100    IF (ABS(X1-X0) .GE. EPS2) GO TO 80
         GROUP = GROUP + 1
         IF (ORDER * (X1 - X0) .LE. 0.0E0) X1 = X0 + ORDER * EPS3
C     .......... EXPAND MATRIX, SUBTRACT EIGENVALUE,
C                AND INITIALIZE VECTOR ..........
  120    DO 200 I = 1, N
            IJ = I + MIN(0,I-M1) * N
            KJ = IJ + MB * N
            IJ1 = KJ + M1 * N
            IF (M1 .EQ. 0) GO TO 180
C
            DO 150 J = 1, M1
               IF (IJ .GT. M1) GO TO 125
               IF (IJ .GT. 0) GO TO 130
               RV(IJ1) = 0.0E0
               IJ1 = IJ1 + N
               GO TO 130
  125          RV(IJ) = A(I,J)
  130          IJ = IJ + N
               II = I + J
               IF (II .GT. N) GO TO 150
               JJ = MB - J
               IF (E21 .GE. 0.0E0) GO TO 140
               II = I
               JJ = MB + J
  140          RV(KJ) = A(II,JJ)
               KJ = KJ + N
  150       CONTINUE
C
  180       RV(IJ) = A(I,MB) - X1
            RV6(I) = EPS4
            IF (ORDER .EQ. 0.0E0) RV6(I) = Z(I,R)
  200    CONTINUE
C
         IF (M1 .EQ. 0) GO TO 600
C     .......... ELIMINATION WITH INTERCHANGES ..........
         DO 580 I = 1, N
            II = I + 1
            MAXK = MIN(I+M1-1,N)
            MAXJ = MIN(N-I,M21-2) * N
C
            DO 360 K = I, MAXK
               KJ1 = K
               J = KJ1 + N
               JJ = J + MAXJ
C
               DO 340 KJ = J, JJ, N
                  RV(KJ1) = RV(KJ)
                  KJ1 = KJ
  340          CONTINUE
C
               RV(KJ1) = 0.0E0
  360       CONTINUE
C
            IF (I .EQ. N) GO TO 580
            U = 0.0E0
            MAXK = MIN(I+M1,N)
            MAXJ = MIN(N-II,M21-2) * N
C
            DO 450 J = I, MAXK
               IF (ABS(RV(J)) .LT. ABS(U)) GO TO 450
               U = RV(J)
               K = J
  450       CONTINUE
C
            J = I + N
            JJ = J + MAXJ
            IF (K .EQ. I) GO TO 520
            KJ = K
C
            DO 500 IJ = I, JJ, N
               V = RV(IJ)
               RV(IJ) = RV(KJ)
               RV(KJ) = V
               KJ = KJ + N
  500       CONTINUE
C
            IF (ORDER .NE. 0.0E0) GO TO 520
            V = RV6(I)
            RV6(I) = RV6(K)
            RV6(K) = V
  520       IF (U .EQ. 0.0E0) GO TO 580
C
            DO 560 K = II, MAXK
               V = RV(K) / U
               KJ = K
C
               DO 540 IJ = J, JJ, N
                  KJ = KJ + N
                  RV(KJ) = RV(KJ) - V * RV(IJ)
  540          CONTINUE
C
               IF (ORDER .EQ. 0.0E0) RV6(K) = RV6(K) - V * RV6(I)
  560       CONTINUE
C
  580    CONTINUE
C     .......... BACK SUBSTITUTION
C                FOR I=N STEP -1 UNTIL 1 DO -- ..........
  600    DO 630 II = 1, N
            I = N + 1 - II
            MAXJ = MIN(II,M21)
            IF (MAXJ .EQ. 1) GO TO 620
            IJ1 = I
            J = IJ1 + N
            JJ = J + (MAXJ - 2) * N
C
            DO 610 IJ = J, JJ, N
               IJ1 = IJ1 + 1
               RV6(I) = RV6(I) - RV(IJ) * RV6(IJ1)
  610       CONTINUE
C
  620       V = RV(I)
            IF (ABS(V) .GE. EPS3) GO TO 625
C     .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM ..........
            IF (ORDER .EQ. 0.0E0) IERR = -R
            V = SIGN(EPS3,V)
  625       RV6(I) = RV6(I) / V
  630    CONTINUE
C
         XU = 1.0E0
         IF (ORDER .EQ. 0.0E0) GO TO 870
C     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
C                MEMBERS OF GROUP ..........
         IF (GROUP .EQ. 0) GO TO 700
C
         DO 680 JJ = 1, GROUP
            J = R - GROUP - 1 + JJ
            XU = 0.0E0
C
            DO 640 I = 1, N
  640       XU = XU + RV6(I) * Z(I,J)
C
            DO 660 I = 1, N
  660       RV6(I) = RV6(I) - XU * Z(I,J)
C
  680    CONTINUE
C
  700    NORM = 0.0E0
C
         DO 720 I = 1, N
  720    NORM = NORM + ABS(RV6(I))
C
         IF (NORM .GE. 0.1E0) GO TO 840
C     .......... IN-LINE PROCEDURE FOR CHOOSING
C                A NEW STARTING VECTOR ..........
         IF (ITS .GE. N) GO TO 830
         ITS = ITS + 1
         XU = EPS4 / (UK + 1.0E0)
         RV6(1) = EPS4
C
         DO 760 I = 2, N
  760    RV6(I) = XU
C
         RV6(ITS) = RV6(ITS) - EPS4 * UK
         GO TO 600
C     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
  830    IERR = -R
         XU = 0.0E0
         GO TO 870
C     .......... NORMALIZE SO THAT SUM OF SQUARES IS
C                1 AND EXPAND TO FULL ORDER ..........
  840    U = 0.0E0
C
         DO 860 I = 1, N
  860    U = U + RV6(I)**2
C
         XU = 1.0E0 / SQRT(U)
C
  870    DO 900 I = 1, N
  900    Z(I,R) = RV6(I) * XU
C
         X0 = X1
  920 CONTINUE
C
 1001 RETURN
      END