C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL TRIDIAGONAL
C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RTSUMARY).
C
C     THE DIMENSION OF  A  SHOULD BE  NM  BY  3.
C     THE DIMENSION OF  Z  SHOULD BE  NM  BY  NM.
C     THE DIMENSION OF  W,D,E,E2,IND,RV1,RV2,RV3,RV4,RV5,RV6,
C     W1,  AND  W2  SHOULD BE  NM.
C     HERE NM = 20.
C
      REAL Z( 20, 20),A( 20, 3),W( 20),D( 20),E( 20),
     X        E2( 20),RV1( 20),RV2( 20),RV3( 20),RV4( 20),RV5( 20),
     X        RV6( 20),W1( 20),W2( 20),TCRIT( 8),EPSLON,RESDUL,MAXEIG,
     X        MAXDIF,U,LB,UB,EPS1,DFL
      REAL  XUB,XLB
      INTEGER  IND( 20),IERR( 6),ERROR
      DATA IREAD1/1/,IREADC/5/,IWRITE/6/
C
      OPEN(IREAD1,FILE='FILE39')
      OPEN(IREADC,FILE='FILE40')
      REWIND IREAD1
      REWIND IREADC
C
      NM = 20
      LCOUNT = 0
      WRITE(*,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/96H ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB
     X      UB    M  IMTQLV   TSTURM   BISECT  M1 NO  TRIDIB  /1H ,
     X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST 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  AND  TSTURM.  //
     X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY  BISECT,
     X30H  AND  TSTURM  THAT LIE IN THE    /18H INTERVAL (LB,UB).//
     X95H UNDER EACH OF 'IMTQLV', 'TSTURM', 'BISECT', AND 'TRIDIB' ARE T
     XWO NUMBERS.  THE FIRST NUMBER,       )
      WRITE(*,2)
    2 FORMAT(
     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   FIGI2-TQL2  .     /
     X62H THE  TQLRAT  PATH USES THE EISPACK CODES   FIGI -TQLRAT.     /
     X62H THE  IMTQL2  PATH USES THE EISPACK CODES   FIGI2-IMTQL2,     /
     X38H AS CALLED FROM DRIVER SUBROUTINE  RT. /
     X62H THE  IMTQL1  PATH USES THE EISPACK CODES   FIGI -IMTQL1,     /
     X38H AS CALLED FROM DRIVER SUBROUTINE  RT. /
     X63H THE  IMTQLV  PATH USES THE EISPACK CODES   FIGI -IMTQLV-TINVIT
     X ,8H-BAKVEC.)
      WRITE(*,3)
    3 FORMAT(
     X64H THE  TSTURM  PATH USES THE EISPACK CODES   FIGI -TSTURM-BAKVEC
     X.  /
     X63H THE  BISECT  PATH USES THE EISPACK CODES   FIGI -BISECT-TINVIT
     X ,8H-BAKVEC. /
     X63H THE  TRIDIB  PATH USES THE EISPACK CODES   FIGI -TRIDIB-TINVIT
     X ,8H-BAKVEC. /)
      WRITE(*,15)
   15 FORMAT(1X,21HS.P. VERSION 04/15/83 )
    5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
     X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
     X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
     X    33H FOR REAL TRIDIAGONAL MATRICES.      /
     X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
     X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
   10 CALL RMATIN(NM,N,A)
      READ(IREADC,50) MM,LB,UB,M11,NO
   50 FORMAT(I4,2D24.16,2(4X,I4))
C
C     MM,LB,UB,M11,  AND  NO  ARE READ FROM SYSIN AFTER THE MATRIX IS
C     GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE MAXIMUM
C     NUMBER OF EIGENVALUES AND BOUNDS FOR THE INTERVAL WHICH IF TO
C     BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE LOWER BOUNDARY
C     INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
C
      DO  230  ICALL = 1,10
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,90,95,230,110,230),  ICALL
C
C     RTWZ  USING  TQL2
C
   70    ICT = 1
         CALL  FIGI2(NM,N,A,W,E,Z,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 200
         CALL  TQL2(NM,N,W,E,Z,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 190
         DO 71 I = 1,N
            W1(I) = W(I)
   71    CONTINUE
         M = N
         GO TO  190
C
C     RTWZ  USING  IMTQL2
C     INVOKED FROM DRIVER SUBROUTINE  RT.
C
   75    ICT = 2
         CALL  RT(NM,N,A,W,1,Z,E,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .GT. N )  GO TO  200
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 190
         DO 78 I = 1,N
            W2(I) = W(I)
   78    CONTINUE
         M = N
         GO TO  190
C
C     RTW  USING  TQLRAT
C
   80    ICT = 7
         IF( IERR(1) .NE. 0 ) GO TO 200
         CALL  FIGI(NM,N,A,W,E,E2,ERROR)
         CALL TQLRAT(N,W,E2,ERROR)
         IERR(1) = ERROR
         IF( ERROR .NE. 0 ) GO TO 200
         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(7) = MAXDIF/EPSLON(MAXEIG*DFL)
         GO TO  230
C
C     RTW  USING  IMTQL1
C     INVOKED FROM DRIVER SUBROUTINE  RT.
C
   85    ICT = 8
         IF( IERR(2) .NE. 0 ) GO TO 200
         CALL  RT(NM,N,A,W,0,Z,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(8) = MAXDIF/EPSLON(MAXEIG*DFL)
         GO TO  230
C
C     RTW1Z  USING  ( USAGE HERE COMPUTES ALL THE EIGENVALUES )
C
   89    ICT = 3
         CALL  FIGI(NM,N,A,D,E,E2,ERROR)
         IERR(ICT) = IABS(ERROR)
         IF( ERROR .NE. 0 ) GO TO 200
         CALL IMTQLV(N,D,E,E2,W,IND,ERROR,RV1)
         IERR(ICT) = ERROR
         M = N
         IF( ERROR .NE. 0 ) M = ERROR - 1
         CALL TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IERR(ICT) + IABS(ERROR)
         CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
         IERR(ICT) = IERR(ICT) + ERROR
         GO TO 190
C
C     RT1W1Z  USING  TSTURM
C
   90    ICT = 4
         EPS1 = 0.0E0
         CALL  FIGI(NM,N,A,D,E,E2,ERROR)
         IERR(ICT) = IABS(ERROR)
         IF( ERROR .NE. 0 )  GO TO  200
         CALL  TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,ERROR,
     X                RV1,RV2,RV3,RV4,RV5,RV6)
         IERR(ICT) = ERROR
         XLB = LB
         XUB = UB
         IERR(ICT) = ERROR
         IF( ERROR .EQ. 3*N + 1 ) GO TO 200
         IF( ERROR .GT. 4*N ) M = ERROR - 4*N - 1
         CALL  BAKVEC(NM,N,A,E,M,Z,ERROR)
         IERR(ICT) = IERR(ICT) + ERROR
         GO TO  190
C
C     RT1W1Z  USING  BISECT  AND  TINVIT
C
   95    ICT = 5
         EPS1 = 0.0E0
         CALL  FIGI(NM,N,A,D,E,E2,ERROR)
         IERR(ICT) = IABS(ERROR)
         IF( ERROR .NE. 0 )  GO TO  200
         CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV4,RV5)
         IERR(ICT) = ERROR
         MBISCT = M
         XLB = LB
         XUB = UB
         IF( ERROR .NE. 0 ) GO TO 200
         CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IABS(ERROR)
         CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
         IERR(ICT) = IERR(ICT) + ERROR
         GO TO  190
C
C     RT1W1Z  USING  TRIDIB  AND  TINVIT
C
  110    ICT = 6
         EPS1 = 0.0E0
         CALL  FIGI(NM,N,A,D,E,E2,ERROR)
         IERR(ICT) = IABS(ERROR)
         IF( ERROR .NE. 0 )  GO TO  200
         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  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IABS(ERROR)
         CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
         IERR(ICT) = IERR(ICT) + ERROR
C
  190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
         CALL  RTWZR(NM,N,M,A,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(*,5)
      LCOUNT = LCOUNT + 1
      WRITE(*,240) N,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2),
     X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,5),
     X             M11,NO,IERR(6),TCRIT(6)
  240 FORMAT(I4,2(I3,2F6.3),2(1PE8.0),I3,3(I3,0PF6.3),3I3,F6.3)
      GO TO  10
      END
      SUBROUTINE RTWZR(NM,N,M,A,W,Z,NORM,RESDUL)
C
      REAL NORM(M), W(M), A(NM,3), 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 NON-SYMMETRIC TRIDIAGONAL
C     MATRIX,  W    IS A VECTOR WHICH CONTAINS  M  EIGENVALUES OF  A,
C     AND  Z  IS AN ARRAY WHICH CONTAINS THE  M  CORRESPONDING
C     EIGENVECTORS OF  A.  ALL NORMS APPEARING IN THE COMMENTS BELOW
C     ARE 1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RTWZR).
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        A(N,3) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
C           SUBDIAGONAL,DIAGONAL AND SUPERDIAGONAL OF THE SYMMETRIC
C           TRIDIAGONAL MATRIX.  THE SUBDIAGONAL BEGINS AT  A(2,1), THE
C           SUPERDIAGONAL BEGINS AT  A(1,3)  AND  A(1,1)  AND  A(N,3)
C           ARE ARBITRARY;
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 = MIN0(N,I+1)
         SUMA = 0.0E0
         LSTART = MAX0(1,I+2-N)
         LSTOP = MIN0(3,I+1)
C
         DO 10 L=LSTART,LSTOP
            SUMA = SUMA + ABS(A(J,L))
   10       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-2)
            KSTOP = MIN0(3,N+2-L)
            KSTART = MAX0(1,3-L)
C
            DO 50 K=KSTART,KSTOP
               J = J+1
   50          SUM = SUM + A(L,K)*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,N,A)
C
C     THIS INPUT SUBROUTINE READS A REAL TRIDIAGONAL MATRIX
C     FROM SYSIN OF ORDER N.
C
C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RTREADI).
C
      REAL A(NM,3)
      INTEGER  IA(3)
      DATA IREADA/1/,IWRITE/6/
C
      READ(IREADA,5) N
    5 FORMAT(I6)
      IF( N .EQ. 0 )  GO TO  70
      DO  15  I = 1,N
         READ(IREADA,10) (IA(J),J=1,3)
   10    FORMAT(3I12)
         DO  15  J = 1,3
   15      A(I,J) = IA(J)
      RETURN
   70 WRITE(*,80)
   80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RTREADI). /1H1)
      STOP
      END
