C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL SYMMETRIC
C     BAND MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RSBSUMAR).
C
C     THE DIMENSION OF  ST  AND  AA  SHOULD BE  NM  BY  MMB.
C     THE DIMENSION OF  Z  SHOULD BE  NM  BY  NM.
C     THE DIMENSION OF  W,W1,W2,D,E,E2,RV3,RV4,RV5,RV6,
C     AND  IND  SHOULD BE  NM.
C     THE DIMENSION OF  STHOLD  SHOULD BE  NM  BY  MMB.
C     THE DIMENSION OF  RV1  SHOULD BE  2*MMB**2+4*MMB-3.
C     THE DIMENSION OF  RV2  SHOULD BE  NM*(2*MMB-1).
C     HERE NM = 44  AND  MMB = 5.
C
C     4-28-92:  MODIFIED CALLS TO RSB AND BANDR TO PASS SEPARATE ARRAYS
C               TO THE DUMMY ARGUMENTS E AND E2.  (ECA)
C
      REAL ST( 44, 5),Z( 44, 44),AA(  44, 5),
     X        STHOLD( 44, 5),W( 44),D( 44),E( 44),
     X        E2( 44),W1( 44),W2( 44),RV1( 67),RV2( 396),
     X        RV4( 44),RV5( 44),RV6( 44),TCRIT( 8),EPSLON,RESDUL,
     X        MAXDIF,MAXEIG,U,LB,UB,EPS1,T,R,TSTART,DFL
      REAL  XUB,XLB
      INTEGER  IND( 44),IERR( 6),ERROR
      DATA IREAD1/1/,IREADC/5/,IWRITE/6/
C
      OPEN(UNIT=IREAD1,FILE='FILE37')
      OPEN(UNIT=IREADC,FILE='FILE38')
      REWIND IREAD1
      REWIND IREADC
C
      NM = 44
      MMB = 5
      NV1 = 67
      NV2 = 396
      LCOUNT = 0
      WRITE(IWRITE,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/  96H ORDER B.W. TQL2   TQLRAT IMTQL2 IMTQL
     X1    LB      UB    M  IMTQLV   BISECT M1 NO  TRIDIB   BQR  /1H ,
     X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
     X51H UNDER 'B.W.' IS THE HALF BAND WIDTH OF THE MATRIX. //
     X95H UNDER 'TQL2   TQLRAT' ARE THREE NUMBERS.  THE FIRST NUMBER, AN
     X INTEGER, IS THE ABSOLUTE SUM OF/
     X61H THE ERROR FLAGS RETURNED SEPARATELY FROM  TQL2  AND  TQLRAT.,
     X34H  THE SECOND NUMBER IS THE MEASURE/
     X62H OF PERFORMANCE BASED UPON THE RESIDUAL COMPUTED FOR THE  TQL2,
     X25H  PATH.  THE THIRD NUMBER        /
     X62H MEASURES THE AGREEMENT OF THE EIGENVALUES FROM THE  TQL2  AND,
     X16H  TQLRAT  PATHS.  //
     X95H UNDER 'IMTQL2 IMTQL1' ARE THREE NUMBERS WITH MEANING LIKE THOS
     XE UNDER  'TQL2   TQLRAT'.       //
     X95H UNDER 'LB' AND 'UB' ARE THE INPUT VARIABLES SPECIFYING THE INT
     XERVAL TO  BISECT.               //
     X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY  BISECT,
     X17H  THAT LIE IN THE                 /18H INTERVAL (LB,UB)./)
      WRITE(IWRITE,2)
    2 FORMAT(
     X62H UNDER EACH OF 'IMTQLV', 'BISECT', 'TRIDIB', AND 'BQR' ARE TWO,
     X28H NUMBERS.  THE FIRST NUMBER,     /
     X95H AN INTEGER, IS THE ABSOLUTE SUM OF THE ERROR FLAGS RETURNED FR
     XOM THE RESPECTIVE PATH.         /
     X95H THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE
     X RESIDUAL COMPUTED FOR THE PATH.//
     X95H UNDER 'M1' AND 'NO' ARE THE VARIABLES SPECIFYING THE LOWER BOU
     XNDARY INDEX AND THE NUMBER       /
     X28H OF EIGENVALUES TO  TRIDIB.   //
     X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
     X27H THE CORRESPONDING PATH HAS        /
     X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS. //
     X62H THE  TQL2    PATH USES THE EISPACK CODES   BANDR-TQL2  ,     /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RSB. /
     X62H THE  TQLRAT  PATH USES THE EISPACK CODES   BANDR-TQLRAT,     /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RSB. /
     X62H THE  IMTQL2  PATH USES THE EISPACK CODES   BANDR-IMTQL2.     )
      WRITE(IWRITE,3)
    3 FORMAT(
     X62H THE  IMTQL1  PATH USES THE EISPACK CODES   BANDR-IMTQL1.     /
     X64H THE  IMTQLV  PATH USES THE EISPACK CODES   BANDR-IMTQLV-BANDV.
     X   /
     X64H THE  BISECT  PATH USES THE EISPACK CODES   BANDR-BISECT-BANDV.
     X  /
     X64H THE  TRIDIB  PATH USES THE EISPACK CODES   BANDR-TRIDIB-BANDV.
     X  /
     X62H THE  BQR     PATH USES THE EISPACK CODES   BQR  -BANDV.      )
      WRITE(IWRITE,15)
   15 FORMAT(1X,21HS.P. VERSION 04/15/83 )
    5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
     X    30HMEASURE OF PERFORMANCE  Y  FOR /5X,
     X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
     X    38H FOR REAL SYMMETRIC BAND MATRICES.        /
     X    58H0ORDER B.W. TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M,
     X    38H  IMTQLV   BISECT M1 NO  TRIDIB   BQR  )
   10 CALL RMATIN(NM,MMB,N,MB,ST,STHOLD,0)
      READ(IREADC,50) MM,LB,UB,M11,NO,MBQR,TSTART
   50 FORMAT(I4,2D24.16,3I4,F8.4)
C
C     MM,LB,UB,M11,NO,MBQR,  AND  TSTART  ARE READ FROM SYSIN AFTER THE
C     MATRIX IS GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE
C     MAXIMUM NUMBER OF EIGENVALUES AND THE BOUNDS FOR THE INTERVAL
C     WHICH IS TO BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE
C     LOWER BOUNDARY INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
C     MBQR  AND  TSTART  SPECIFY TO  BQR  THE NUMBER OF EIGENVALUES
C     AND THE VALUE TO WHICH THEY ARE DESIRED CLOSEST.
C
      DO  230  ICALL = 1,10
         IF( ICALL .NE. 1 )  CALL  RMATIN(NM,MMB,N,MB,ST,STHOLD,1)
C
C     IF  TQLRAT  PATH (LABEL 80) IS TAKEN THEN  TQL2  PATH (LABEL 70)
C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
C     MEANINGFUL.
C     IF  IMTQL1  PATH (LABEL 85) IS TAKEN THEN  IMTQL2  PATH (LABEL 75)
C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
C     MEANINGFUL.
C     IF  TQL2  (IMTQL2)  PATH FAILS, THEN  TQLRAT  (IMTQL1)  PATH IS
C     OMITTED AND PRINTOUT FLAGGED WITH  -1.0.
C
         GO TO  (70,75,80,85,89,230,230,120,125,130),  ICALL
C
C     RSBWZ  USING  TQL2
C     INVOKED FROM DRIVER SUBROUTINE  RSB.
C
   70    ICT = 1
         CALL  RSB(NM,N,MB,ST,W,1,Z,E,E2,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 200
         DO 71 I = 1,N
            W1(I) = W(I)
   71    CONTINUE
         M = N
         GO TO  195
C
C     RSBWZ  USING  IMTQL2
C
   75    ICT = 2
         CALL  BANDR(NM,N,MB,ST,W,E,E2,.TRUE.,Z)
         CALL  IMTQL2(NM,N,W,E,Z,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 200
         DO  78  I=1,N
   78      W2(I) = W(I)
         M = N
         GO TO  195
C
C     RSBW  USING  TQLRAT
C     INVOKED FROM DRIVER SUBROUTINE  RSB.
C
   80    ICT = 7
         IF( IERR(1) .NE. 0 ) GO TO 200
         CALL  RSB(NM,N,MB,ST,W,0,Z,E,E2,ERROR)
         IF( ERROR .NE. 0 ) GO TO 200
         IERR(1) = ERROR
         MAXEIG = 0.0E0
         MAXDIF = 0.0E0
         DO 81 I = 1,N
            IF( ABS(W(I)) .GT. MAXEIG ) MAXEIG = ABS(W(I))
            U = ABS(W1(I) - W(I))
            IF( U .GT. MAXDIF ) MAXDIF = U
   81    CONTINUE
         IF( MAXEIG .EQ. 0.0E0 ) MAXEIG = 1.0E0
         DFL = 10 * N
         TCRIT(ICT) = MAXDIF/EPSLON(MAXEIG*DFL)
         GO TO  230
C
C     RSBW  USING  IMTQL1
C
   85    ICT = 8
         IF( IERR(2) .NE. 0 ) GO TO 200
         CALL  BANDR(NM,N,MB,ST,W,E,E2,.FALSE.,Z)
         CALL  IMTQL1(N,W,E,ERROR)
         IERR(2) = ERROR
         MAXEIG = 0.0E0
         MAXDIF = 0.0E0
         DO 86 I = 1,N
            IF( ABS(W(I)) .GT. MAXEIG ) MAXEIG = ABS(W(I))
            U = ABS(W2(I) - W(I))
            IF( U .GT. MAXDIF ) MAXDIF = U
   86    CONTINUE
         IF( MAXEIG .EQ. 0.0E0 ) MAXEIG = 1.0E0
         DFL = 10 * N
         TCRIT(ICT) = MAXDIF/EPSLON(MAXEIG*DFL)
         GO TO  230
C
C     RSBW1Z  ( USAGE HERE COMPUTES ALL THE EIGENVECTORS )
C
   89    ICT = 3
         DO 894 I = 1,MB
            DO 894 J = 1,N
  894          AA(J,I) = ST(J,I)
         CALL  BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z)
         CALL  IMTQLV(N,D,E,E2,W,IND,ERROR,RV1)
         IERR(ICT) = ERROR
         M = N
         IF( ERROR .NE. 0 ) M = ERROR - 1
         CALL  BANDV(NM,N,MB,AA,0.0E0,M,W,Z,ERROR,NV2,RV2,RV6)
         IERR(ICT) = IERR(ICT) + IABS(ERROR)
         GO TO 195
C
C     RSB1W1Z  USING  BISECT  AND  BANDV
C
  120    ICT = 4
         EPS1 = 0.0E0
         DO 122 I=1,MB
           DO 122 J=1,N
  122        AA(J,I) = ST(J,I)
         CALL  BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z)
         CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV4,RV5)
         MBISCT = M
         XLB = LB
         XUB = UB
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 200
         CALL  BANDV(NM,N,MB,AA,0.0E0,M,W,Z,ERROR,NV2,RV2,RV6)
         IERR(ICT) = IABS(ERROR)
         GO TO 195
C
C     RSB1W1Z  USING  TRIDIB  AND  BANDV
C
  125    ICT = 5
         EPS1 = 0.0E0
         DO 126 I=1,MB
           DO 126 J=1,N
  126        AA(J,I) = ST(J,I)
         CALL  BANDR(NM,N,MB,ST,D,E,E2,.FALSE.,Z)
         CALL  TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,NO,W,IND,ERROR,RV4,RV5)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 200
         M = NO
         CALL  BANDV(NM,N,MB,AA,0.0E0,M,W,Z,ERROR,NV2,RV2,RV6)
         IERR(ICT) = IABS(ERROR)
         GO TO 195
C
C     RSB1W1Z  USING  BQR  AND  BANDV
C
  130    ICT = 6
         DO 132 I = 1,MB
            DO 132 J = 1,N
  132          AA(J,I) = ST(J,I)
         N1 = N
         T = TSTART
         R = 0.0E0
         DO 133 J = 1,N
  133       ST(J,MB) = ST(J,MB) - T
         DO 138 I = 1,MBQR
            CALL  BQR(NM,N1,MB,ST,T,R,ERROR,NV1,RV1)
            IF( ERROR .EQ. 0 ) GO TO 134
            M = I - 1
            GO TO 139
  134       IF( I .EQ. 1 ) GO TO 136
            DO 135 JJ = 2,I
               J = I + 2 - JJ
               IF( T .GT. W(J - 1) ) GO TO 137
               W(J) = W(J - 1)
  135       CONTINUE
  136       J = 1
  137       W(J) = T
            N1 = N1 - 1
  138    CONTINUE
         M = MBQR
  139    IERR(ICT) = ERROR
         CALL  BANDV(NM,N,MB,AA,0.0E0,M,W,Z,ERROR,NV2,RV2,RV6)
         IERR(ICT) = IERR(ICT) + IABS(ERROR)
         GO TO 195
C
  195    CALL  RMATIN(NM,MMB,N,MB,ST,STHOLD,1)
         CALL  RSBWZR(NM,N,M,MB,ST,W,Z,RV1,RESDUL)
         DFL = 10 * N
         TCRIT(ICT) = RESDUL/EPSLON(DFL)
         GO TO 230
  200    TCRIT(ICT) = -1.0E0
  230 CONTINUE
C
      IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
      LCOUNT = LCOUNT + 1
      WRITE(IWRITE,520) N,MB,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2),
     X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,4),
     X             M11,NO,(IERR(I),TCRIT(I),I=5,6)
  520 FORMAT(I4,2X,I3,2(I3,2F6.3),2(1PE8.0),I3,2(I3,0PF6.3),
     X       I2,1X,I2,2(I3,0PF6.3))
      GO TO  10
      END
      SUBROUTINE RSBWZR(NM,N,M,MB,A,W,Z,NORM,RESDUL)
C
      INTEGER  MB,MB1
      REAL NORM(M), W(M), A(NM,MB), Z(NM,M), NORMA, TNORM,
     X       S, SUM, SUMA, SUMZ, RESDUL
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-Z*DIAG(W)  WHERE  A  IS A REAL SYMMETRIC BAND MATRIX,
C     W  IS A VECTOR WHICH CONTAINS  M  EIGENVALUES OF  A, AND  Z
C     IS AN ARRAY WHICH CONTAINS THE  M  CORRESPONDING EIGENVECTORS OF
C     A.  ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSBWZR).
C
C     INPUT.
C
C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX  A;
C
C        M IS THE NUMBERS OF EIGENVECTORS WHOSE RESIDUALS ARE DESIRED;
C
C        MB  IS THE BAND WIDTH OF THE INPUT MATRIX  A .  BAND WIDTH IS
C            DEFINED AS THE NUMBER OF ADJACENT DIAGONALS, INCLUDING THE
C            PRINCIPAL DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO
C            PORTION OF THE LOWER TRIANGLE OF THE MATRIX;
C
C        A(N,MB) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
C           SUBDIAGONAL AND DIAGONAL OF THE SYMMETRIC BAND
C           MATRIX;
C
C        W(M) IS A VECTOR WHOSE FIRST  M  COMPONENTS CONTAIN EIGENVALUES
C           OF  A;
C
C        Z(NM,M) IS AN ARRAY WHOSE FIRST  M  COLUMNS CONTAIN THE
C           EIGENVECTORS OF  A  CORRESPONDING TO THE EIGENVALUES IN  W.
C
C     OUTPUT.
C
C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
C           APPROXIMATE EIGENVECTORS OF  A.  THE EIGENVECTORS
C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY
C           THAT THE FIRST ELEMENT WHOSE MAGNITUDE IS LARGER
C           THAN THE NORM OF THE EIGENVECTOR DIVIDED BY  N  IS
C           POSITIVE;
C
C        NORM(M) IS AN ARRAY SUCH THAT FOR EACH  K
C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
C
C        RESDUL IS THE REAL NUMBER
C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
C
C     ----------------------------------------------------------------
C
      RESDUL = 0.0E0
      IF( M .EQ. 0 ) RETURN
      NORMA = 0.0E0
C
      DO 40 I=1,N
         J=I
         SUMA = 0.0E0
         IF(I .EQ. 1) GO TO 20
         MB1 = MB-1
         LSTART = MAX0(1,MB+1-I)
C
         DO 10 L=LSTART,MB1
   10       SUMA = SUMA + ABS(A(I,L))
C
   20    LSTOP = MIN0(MB,N+1-J)
C
         DO 30 L=1,LSTOP
            L1 = MB + 1 - L
            SUMA = SUMA + ABS(A(J,L1))
   30       J = J+1
C
   40    NORMA = AMAX1(SUMA,NORMA)
C
      IF(NORMA .EQ. 0.0E0) NORMA = 1.0E0
C
      DO 120 I=1,M
         S = 0.0E0
         SUMZ = 0.0E0
C
         DO 80 L=1,N
            SUM = -W(I)*Z(L,I)
            SUMZ = SUMZ + ABS(Z(L,I))
            J = MAX0(0,L-MB)
            IF(L .EQ. 1) GO TO 60
            MB1 = MB-1
            KSTART = MAX0(1,MB+1-L)
C
            DO 50 K=KSTART,MB1
               J = J+1
   50          SUM = SUM + A(L,K)*Z(J,I)
C
   60       KSTOP = MIN0(MB,N+1-L)
C
            DO 70 K=1,KSTOP
               J = J+1
               K1 = MB + 1 - K
   70          SUM = SUM + A(J,K1)*Z(J,I)
C
   80       S = S + ABS(SUM)
C
         NORM(I) = SUMZ
         IF (SUMZ .EQ. 0.0E0) GO TO 120
C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
C                  LARGER THAN !!Z(I)!!/N..........
         DO 90 L=1,N
            IF(ABS(Z(L,I)) .GE. NORM(I)/N) GO TO 100
   90       CONTINUE
C
  100    TNORM = SIGN(NORM(I),Z(L,I))
C
         DO 110 L=1,N
  110       Z(L,I) = Z(L,I)/TNORM
C
         NORM(I) = S/(NORM(I)*NORMA)
  120    RESDUL = AMAX1(NORM(I),RESDUL)
C
      RETURN
      END
      SUBROUTINE RMATIN(NM,MMB,N,MB,ST,STHOLD,INITIL)
C
C     THIS INPUT SUBROUTINE READS A REAL SYMMETRIC BAND MATRIX
C     FROM SYSIN OF ORDER  N,  AND BAND WIDTH  MB .
C     TO GENERATE THE MATRIX  ST  INITIALLY,  INITIL  IS TO BE 0.
C     TO REGENERATE THE MATRIX  ST  FOR THE PURPOSE OF THE RESIDUAL
C     CALCULATION,  INITIL  IS TO BE  1.
C
C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSBREADI).
C
      REAL ST(NM,MMB),STHOLD(NM,MMB)
      INTEGER  IA( 5)
      DATA IREADA/1/,IWRITE/6/
C
      IF( INITIL .EQ. 1 )  GO TO  30
      READ(IREADA,5) N,MB
    5 FORMAT(2I6)
      IF( N .EQ. 0 )  GO TO  70
      DO 8 I = 1,N
         DO 7 J = 1,MB
            ST(I,J) = 0.0E0
    7    CONTINUE
    8 CONTINUE
      DO 15 I=1,N
         MBB = MIN0(MB,N-I+1)
         READ(IREADA,10) (IA(J),J=1,MBB)
   10    FORMAT(6I12)
         DO 11 J=1,MBB
           M = MB+1-J
           K = I+J-1
   11      ST(K,M) = IA(J)
   15 CONTINUE
      DO  20  I = 1,N
         DO  20  J = 1,MB
   20      STHOLD(I,J) = ST(I,J)
      RETURN
   30 DO  40  I = 1,N
         DO  40  J = 1,MB
   40      ST(I,J) = STHOLD(I,J)
      RETURN
   70 WRITE(IWRITE,80)
   80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RSBREADI).  /1H1)
      STOP
      END
