C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF COMPLEX HERMITIAN
C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(CHSUMARY).
C
C     THE DIMENSION OF  AR,AI,ZR,  AND  ZI  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  TAU  SHOULD BE  2  BY  NM.
C     THE DIMENSION OF  ARHOLD  AND  AIHOLD  SHOULD BE  NM  BY  NM.
C     HERE NM = 20.
C
C     4-28-92:  MODIFIED CALLS TO CH AND HTRIDI TO PASS SEPARATE ARRAYS
C               TO THE DUMMY ARGUMENTS E AND E2.  (ECA)
C
      REAL AR( 20, 20),AI( 20, 20),ZR( 20, 20),ZI( 20, 20),
     X        TAU( 2, 20),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,
     X        MAXEIG,MAXDIF,U,LB,UB,EPS1,DFL
      REAL ARHOLD( 20, 20),AIHOLD( 20, 20)
      REAL XUB,XLB
      INTEGER  IND( 20),IERR( 6),ERROR
      DATA IREAD1/1/,IREADC/5/,IWRITE/6/
C
      OPEN(UNIT=IREAD1,FILE='FILE43')
      OPEN(UNIT=IREADC,FILE='FILE44')
      REWIND IREAD1
      REWIND IREADC
C
      NM = 20
      LCOUNT = 0
      WRITE(IWRITE,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/  34H ORDER TQL2   TQLRAT IMTQL2 IMTQL1,4X,
     X56HLB      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      /
     X27H 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. //
     X63H THE  TQL2    PATH USES THE EISPACK CODES  HTRIDI-TQL2  -HTRIBK
     X,1H, /
     X38H AS CALLED FROM DRIVER SUBROUTINE  CH. /
     X62H THE  TQLRAT  PATH USES THE EISPACK CODES  HTRIDI-TQLRAT,     /
     X38H AS CALLED FROM DRIVER SUBROUTINE  CH. /
     X63H THE  IMTQL2  PATH USES THE EISPACK CODES  HTRIDI-IMTQL2-HTRIBK
     X,1H.  )
      WRITE(IWRITE,3)
    3 FORMAT(
     X62H THE  IMTQL1  PATH USES THE EISPACK CODES  HTRIDI-IMTQL1.     /
     X63H THE  IMTQLV  PATH USES THE EISPACK CODES  HTRID3-IMTQLV-TINVIT
     X ,8H-HTRIB3./
     X64H THE  TSTURM  PATH USES THE EISPACK CODES  HTRIDI-TSTURM-HTRIBK
     X.  /
     X63H THE  BISECT  PATH USES THE EISPACK CODES  HTRIDI-BISECT-TINVIT
     X ,8H-HTRIBK. /
     X63H THE  TRIDIB  PATH USES THE EISPACK CODES  HTRIDI-TRIDIB-TINVIT
     X ,8H-HTRIBK. /)
      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    33H FOR COMPLEX HERMITIAN MATRICES.     /
     X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
     X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
   10 CALL CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,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  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,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     CHWZ  USING  TQL2
C     INVOKED FROM DRIVER SUBROUTINE  CH.
C
   70    ICT = 1
         CALL  CH(NM,N,AR,AI,W,1,ZR,ZI,E,E2,TAU,ERROR)
         IERR(ICT) = ERROR
         M = ERROR - 1
         IF( ERROR .NE. 0 ) GO TO 74
         DO 71 I = 1,N
            W1(I) = W(I)
   71    CONTINUE
         M = N
   74    GO TO  190
C
C     CHWZ  USING  IMTQL2
C
   75    ICT = 2
         DO 77 I = 1,N
            DO  76  J = 1,N
   76         ZR(I,J) = 0.0E0
   77       ZR(I,I) = 1.0E0
         CALL  HTRIDI(NM,N,AR,AI,W,E,E2,TAU)
         CALL  IMTQL2(NM,N,W,E,ZR,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  HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
         GO TO  190
C
C     CHW  USING  TQLRAT
C     INVOKED FROM DRIVER SUBROUTINE  CH.
C
   80    ICT = 7
         IF( IERR(1) .NE. 0 ) GO TO 200
         CALL  CH(NM,N,AR,AI,W,0,AR,AI,E,E2,TAU,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     CHW  USING  IMTQL1
C
   85    ICT = 8
         IF( IERR(2) .NE. 0 ) GO TO 200
         CALL  HTRIDI(NM,N,AR,AI,W,E,E2,TAU)
         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     CHW1Z  ( USAGE HERE COMPUTES ALL THE EIGENVECTORS )
C
   89    ICT = 3
         DO 892 I = 2,N
            IM1 = I - 1
            DO 891 J = 1,IM1
               AR(J,I) = AI(I,J)
  891       CONTINUE
  892    CONTINUE
         CALL  HTRID3(NM,N,AR,D,E,E2,TAU)
         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,ZR,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IERR(ICT) + IABS(ERROR)
         CALL  HTRIB3(NM,N,AR,TAU,M,ZR,ZI)
         CALL  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,1)
         GO TO 190
C
C     CH1W1Z  USING  TSTURM
C
   90    ICT = 4
         EPS1 = 0.0E0
         CALL  HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
         CALL  TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,ZR,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  HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
         GO TO  190
C
C     CH1W1Z  USING  BISECT  AND  TINVIT
C
   95    ICT = 5
         EPS1 = 0.0E0
         CALL  HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
         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,ZR,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IABS(ERROR)
         CALL  HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
         GO TO  190
C
C     CH1W1Z  USING  TRIDIB  AND  TINVIT
C
  110    ICT = 6
         EPS1 = 0.0E0
         CALL  HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
         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,ZR,ERROR,RV1,RV2,RV3,RV4,RV6)
         IERR(ICT) = IABS(ERROR)
         CALL  HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
C
  190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
         DO 195 I = 1,N
            AI(I,I) = 0.0E0
  195    CONTINUE
         CALL  CHWZR(NM,N,M,AR,AI,W,ZR,ZI,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 CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,INITIL)
C
C     THIS INPUT SUBROUTINE READS A COMPLEX MATRIX  A = (AR,AI)
C     FROM SYSIN OF ORDER N.
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(CGREADI).
C
      REAL AR(NM,NM),AI(NM,NM),ARHOLD(NM,NM),AIHOLD(NM,NM)
      INTEGER  IAR( 20), IAI( 20)
      DATA IREADA/1/,IWRITE/6/
C
      IF( INITIL .EQ. 1 )  GO TO  30
      READ(IREADA,5) N
    5 FORMAT(I6)
      IF( N .EQ. 0 )  GO TO  70
      DO  15  I = 1,N
         READ(IREADA,10) (IAR(J),IAI(J),J=1,N)
   10    FORMAT(2I18)
         DO  15  J = 1,N
           AR(I,J) = IAR(J)
   15      AI(I,J) = IAI(J)
      DO  20  I = 1,N
         DO  20  J = 1,N
           ARHOLD(I,J) = AR(I,J)
   20      AIHOLD(I,J) = AI(I,J)
      RETURN
   30 DO  40  I = 1,N
         DO  40  J = 1,N
           AR(I,J) = ARHOLD(I,J)
   40      AI(I,J) = AIHOLD(I,J)
      RETURN
   70 WRITE(IWRITE,80)
   80 FORMAT(44H0END OF DATA FOR SUBROUTINE CMATIN(CGREADI).  /1H1)
      STOP
      END
      SUBROUTINE CHWZR(NM,N,M,AR,AI,W,ZR,ZI,NORM,RESDUL)
C
      REAL NORM(M),W(M),AR(NM,N),AI(NM,N),
     X       ZR(NM,M),ZI(NM,M),NORMA,XR,XI,S,SUMA,SUMZ,SUMR,SUMI,RESDUL
      REAL PYTHAG
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-Z*DIAG(W)  WHERE  A  IS A HERMITIAN MATRIX,  W IS
C     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(CHWZR).
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 NUMBER OF EIGENVECTORS FOR WHICH RESIDUALS ARE
C           DESIRED;
C
C        AI(NM,N), AR(NM,N) ARE ARRAYS CONTAINING THE REAL AND
C           IMAGINARY PARTS OF  A.  ONLY THE FULL UPPER TRIANGLE
C           NEED BE SUPPLIED;
C
C        W(M) IS A VECTOR WHOSE FIRST  M  COMPONENTS CONTAIN  M
C           EIGENVALUES OF  A;
C
C        ZR(NM,M), ZI(NM,M) ARE ARRAYS WHOSE FIRST  M  COLUMNS CONTAIN
C           THE REAL AND IMAGINARY PARTS OF THE ELEMENTS OF  Z.
C
C     OUTPUT.
C
C        ZR(NM,M), ZI(NM,M) ARE ARRAYS WHOSE COLUMNS CONTAIN THE
C           REAL AND IMAGINARY PARTS OF THE NORMALIZED APPROXIMATE
C           EIGENVECTORS OF  A.  THE EIGENVECTORS ARE NORMALIZED BY
C           THE 1-NORM IN SUCH A WAY THAT THE FIRST ELEMENT WHOSE
C           MAGNITUDE IS LARGER THAN THE NORM OF THE EIGENVECTOR
C           DIVIDED BY  N  IS REAL AND POSITIVE;
C
C        NORM(N) 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
      NORMA = 0.0E0
      RESDUL = 0.0E0
      IF( M .EQ. 0 ) RETURN
C
      DO 40 I=1,N
         SUMA = 0.0E0
         IF(I .EQ. 1) GO TO 20
C
         DO 10 L=2,I
            AR(I,L-1) = AR(L-1,I)
            AI(I,L-1) = -AI(L-1,I)
   10    CONTINUE
C
   20    DO 30 L=1,N
   30       SUMA = SUMA + PYTHAG(AR(I,L),AI(I,L))
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
C
         DO 60 L=1,N
            SUMZ = SUMZ + PYTHAG(ZR(L,I),ZI(L,I))
            SUMR = -W(I)*ZR(L,I)
            SUMI = -W(I)*ZI(L,I)
C
            DO 50 K=1,N
               SUMR = SUMR + AR(L,K)*ZR(K,I) - AI(L,K)*ZI(K,I)
   50          SUMI = SUMI + AR(L,K)*ZI(K,I) + AI(L,K)*ZR(K,I)
C
   60       S = S + PYTHAG(SUMR,SUMI)
C
         NORM(I) = SUMZ
         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(PYTHAG(ZR(L,I),ZI(L,I)) .GE. NORM(I)/N)
     1      GO TO 80
   70       CONTINUE
C
   80    XR = NORM(I)*ZR(L,I)/PYTHAG(ZR(L,I),ZI(L,I))
         XI = NORM(I)*ZI(L,I)/PYTHAG(ZR(L,I),ZI(L,I))
C
         DO 90 L=1,N
            CALL CDIV(ZR(L,I),ZI(L,I),XR,XI,ZR(L,I),ZI(L,I))
   90    CONTINUE
C
         NORM(I) = S/(NORM(I)*NORMA)
  100    RESDUL = AMAX1(NORM(I),RESDUL)
C
      RETURN
      END
