C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL SYMMETRIC
C     PACKED MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RSPSUMAR).
C
C     THE DIMENSION OF  A  SHOULD BE  NNN  AND THE DIMENSION
C     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     THE DIMENSION OF  AHOLD  SHOULD BE  NNN.
C     HERE NM = 20, AND  NNN = NM*(NM+1)/2.
C
C     4-28-92:  MODIFIED CALLS TO RSP AND TRED3 TO PASS SEPARATE ARRAYS
C               TO THE DUMMY ARGUMENTS E AND E2.  (ECA)
C
      REAL A(210),Z(20,20),AHOLD(210),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(UNIT=IREAD1,FILE='FILE35')
      OPEN(UNIT=IREADC,FILE='FILE36')
      REWIND IREAD1
      REWIND IREADC
C
      NM = 20
      NNN = (NM*(NM+1))/2
      LCOUNT = 0
      WRITE(IWRITE,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(IWRITE,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. //
     X64H THE  TQL2    PATH USES THE EISPACK CODES   TRED3-TQL2  -TRBAK3
     X,  /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RSP. /
     X62H THE  TQLRAT  PATH USES THE EISPACK CODES   TRED3-TQLRAT,     /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RSP. /
     X64H THE  IMTQL2  PATH USES THE EISPACK CODES   TRED3-IMTQL2-TRBAK3
     X.  )
      WRITE(IWRITE,3)
    3 FORMAT(
     X62H THE  IMTQL1  PATH USES THE EISPACK CODES   TRED3-IMTQL1.     /
     X63H THE  IMTQLV  PATH USES THE EISPACK CODES   TRED3-IMTQLV-TINVIT
     X ,8H-TRBAK3./
     X64H THE  TSTURM  PATH USES THE EISPACK CODES   TRED3-TSTURM-TRBAK3
     X.  /
     X63H THE  BISECT  PATH USES THE EISPACK CODES   TRED3-BISECT-TINVIT
     X ,8H-TRBAK3. /
     X63H THE  TRIDIB  PATH USES THE EISPACK CODES   TRED3-TRIDIB-TINVIT
     X ,8H-TRBAK3. /)
      WRITE(IWRITE,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    40H FOR REAL SYMMETRIC PACKED MATRICES.        /
     X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
     X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
   10 CALL RMATIN(NNN,N,A,AHOLD,0)
      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 THE BOUNDS FOR THE INTERVAL WHICH IS
C     TO BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE LOWER
C     BOUNDARY INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
C
      DO  230  ICALL = 1,10
         IF( ICALL .NE. 1 )  CALL  RMATIN(NNN,N,A,AHOLD,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,90,95,230,110,230),  ICALL
C
C     RSPWZ  USING  TQL2
C     INVOKED FROM DRIVER SUBROUTINE  RSP.
C
   70    ICT = 1
         CALL  RSP(NM,N,NNN,A,W,1,Z,E,E2,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 74
         M = N
         DO 71 I = 1,N
            W1(I) = W(I)
   71    CONTINUE
   74    GO TO  190
C
C     RSPWZ  USING  IMTQL2
C
   75    ICT = 2
         DO 77 I=1,N
           DO 76 J=1,N
   76        Z(I,J)=0.0E0
   77      Z(I,I)=1.0E0
         CALL  TRED3(N,NNN,A,W,E,E2)
         CALL  IMTQL2(NM,N,W,E,Z,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 79
         DO 78 I = 1,N
   78       W2(I) = W(I)
         M = N
   79    CALL  TRBAK3(NM,N,NNN,A,M,Z)
         GO TO  190
C
C     RSPW  USING  TQLRAT
C     INVOKED FROM DRIVER SUBROUTINE  RSP.
C
   80    ICT = 7
         IF( IERR(1) .NE. 0 ) GO TO 200
         CALL  RSP(NM,N,NNN,A,W,0,Z,E,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     RSPW  USING  IMTQL1
C
   85    ICT = 8
         IF( IERR(2) .NE. 0 ) GO TO 200
         CALL  TRED3(N,NNN,A,W,E,E2)
         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(8) = MAXDIF/EPSLON(MAXEIG*DFL)
         GO TO  230
C
C     RSPW1Z  ( USAGE HERE COMPUTES ALL THE EIGENVECTORS )
C
   89    ICT = 3
         CALL  TRED3(N,NNN,A,D,E,E2)
         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  TRBAK3(NM,N,NNN,A,M,Z)
         GO TO 190
C
C     RSP1W1Z  USING  TSTURM
C
   90    ICT = 4
         EPS1 = 0.0E0
         CALL  TRED3(N,NNN,A,D,E,E2)
         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
         IF( ERROR .EQ. 3*N + 1 ) GO TO 200
         IF( ERROR .GT. 4*N ) M = ERROR - 4*N + 1
         CALL  TRBAK3(NM,N,NNN,A,M,Z)
         GO TO  190
C
C     RSP1W1Z  USING  BISECT  AND  TINVIT
C
   95    ICT = 5
         EPS1 = 0.0E0
         CALL  TRED3(N,NNN,A,D,E,E2)
         CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV1,RV2)
         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  TRBAK3(NM,N,NNN,A,M,Z)
         GO TO  190
C
C     RSP1W1Z  USING  TRIDIB  AND  TINVIT
C
  110    ICT = 6
         EPS1 = 0.0E0
         CALL  TRED3(N,NNN,A,D,E,E2)
         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  TRBAK3(NM,N,NNN,A,M,Z)
C
  190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
         CALL  RMATIN(NNN,N,A,AHOLD,1)
         CALL  RSPWZR(NM,N,NNN,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(IWRITE,5)
      LCOUNT = LCOUNT + 1
      WRITE(IWRITE,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 RSPWZR(NM,N,NNN,M,A,W,Z,NORM,RESDUL)
C
      REAL NORM(M), W(M), A(NNN), Z(NM,M), NORMA, S, SUM,
     X       SUMA, SUMZ, RESDUL, TNORM
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-Z*DIAG(W)   WHERE  A  IS A REAL SYMMETRIC MATRIX STORED IN
C     PACKED FORM,  W  IS A VECTOR WHICH CONTAINS M EIGENVALUES OF
C     A , AND  Z  IS AN ARRAY WHICH CONTAINS THE CORRESPONDING EIGEN-
C     VECTORS OF  A . ALL NORMS APPEARING IN THE COMMENTS BELOW ARE
C     1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSPWZR).
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        NNN  IS THE DIMENSION OF THE ARRAY PARAMETER  A;
C
C        M  IS THE NUMBER OF EIGENVECTORS FOR WHICH RESIDUALS ARE
C           DESIRED;
C
C        A(N*(N+1)/2) IS A VECTOR WHICH CONTAINS THE ELEMENTS OF THE
C           LOWER TRIANGULAR PART OF THE MATRIX  A  (AS MENTIONED ABOVE)
C           IN PACKED FORM.  BY PACKED FORM, WE MEAN THAT THE FIRST ROW
C           OF THE TRIANGLE IS STORED IN THE FIRST POSITION OF  A , THE
C           SECOND ROW OF THE TRIANGLE IS STORED IN THE NEXT TWO
C           POSITIONS, AND SO FORTH UNTIL WE HAVE THE N-TH ROW STORED
C           IN THE LAST N POSITIONS OF  A;
C
C        W(M)  IS AN ARRAY CONTAINING THE EIGENVALUES OF  A;
C
C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE EIGENVECTORS OF  A.
C
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
      INCMT = 0
C
      DO 40 I=1,N
         J=I-1
         SUMA = 0.0E0
         ISP = INCMT
         INCMT = INCMT + I
         LSTOP = N+1-I
C
         IF(I .EQ. 1) GO TO 25
C
         DO 20 L=1,J
           L1 = ISP + L
   20      SUMA = SUMA + ABS(A(L1))
   25    ISP = ISP + 1
         DO 30 L=1,LSTOP
            ISP = ISP + J
            SUMA = SUMA + ABS(A(ISP))
   30       J = J+1
C
   40    NORMA = AMAX1(NORMA,SUMA)
C
      IF(NORMA .EQ. 0.0E0) NORMA = 1.0E0
C
      DO 100 I=1,M
         S = 0.0E0
         SUMZ = 0.0E0
         INCMT = 0
C
         DO 65 L=1,N
            SUMZ = SUMZ + ABS(Z(L,I))
            SUM = -W(I)*Z(L,I)
            J = L-1
            ISP = INCMT
            INCMT = INCMT + L
            KSTOP = N+1-L
C
            IF(L .EQ. 1) GO TO 55
C
            DO 50 K=1,J
              K1 = ISP + K
   50         SUM = SUM + A(K1)*Z(K,I)
   55       ISP = ISP + 1
            DO 60 K=1,KSTOP
              ISP =ISP + J
              K2 = K - 1 + L
              SUM =SUM + A(ISP)*Z(K2,I)
   60       J = J + 1
   65     S = S + ABS(SUM)
C
         NORM(I) = SUMZ
C
         IF(SUMZ .EQ. 0.0E0) GO TO 100
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 70 L=1,N
            IF(ABS(Z(L,I)) .GE. NORM(I)/N) GO TO 80
   70       CONTINUE
C
   80    TNORM = SIGN(NORM(I),Z(L,I))
C
         DO 90 L=1,N
   90       Z(L,I) = Z(L,I)/TNORM
C
         NORM(I) = S/(NORM(I)*NORMA)
  100    RESDUL = AMAX1(NORM(I),RESDUL)
C
      RETURN
      END
      SUBROUTINE RMATIN(NNM,N,A,AHOLD,INITIL)
C
C     THIS INPUT SUBROUTINE READS A REAL SYMMETRIC MATRIX FROM SYSIN OF
C     ORDER N AND STORES IT IN PACKED FORM. THE INPUT MATRIX IS READ
C     AS A FULL MATRIX  (MM IS EQUAL TO ZERO)  OR AS ROWS OF AN UPPER
C     TRIANGULAR MATRIX  (MM IS NOT EQUAL TO ZERO) .
C     TO GENERATE THE MATRIX  A  INITIALLY,  INITIL  IS TO BE 0.
C     TO REGENERATE THE MATRIX  A  FOR THE PURPOSE OF THE RESIDUAL
C     CALCULATION,  INITIL  IS TO BE  1.
C
C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSPREADI).
C
      REAL A(NNM),AHOLD(NNM) 
      INTEGER  IA( 20)
      DATA IREADA/1/,IWRITE/6/
C
      IF( INITIL .EQ. 1 )  GO TO  30
      READ(IREADA,5) N,MM
    5 FORMAT(I6,6X,I6)
      IF( N .EQ. 0 )  GO TO  70
      INCMT=1
      DO 15 L=1,N
         JST = L
         IF( MM .EQ. 0 )  JST = 1
         READ(IREADA,10) (IA(J),J=JST,N) 
   10    FORMAT(6I12)
         J = L - 1
         ISP = INCMT
         INCMT =INCMT + L
         DO 12 K = L,N
           ISP = ISP + J
           A(ISP) = IA(J + 1)
   12      J = J + 1
   15    CONTINUE
      NNN = (N*(N+1))/2
      DO 20 I = 1,NNN
   20   AHOLD(I) = A(I)
      RETURN
   30 NNN = (N*(N+1))/2
      DO 40 I = 1,NNN
   40    A(I) = AHOLD(I)
      RETURN
   70 WRITE(IWRITE,80)
   80 FORMAT(47H0END OF DATA FOR SUBROUTINE  RMATIN(RSPREADI). /1H1)
      STOP
      END
