C
C     THIS DRIVER TESTS  EISPACK  FOR THE REAL GENERAL
C     GENERALIZED EIGENPROBLEM  A*X = (LAMBDA)*B*X, SUMMARIZING THE
C     FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RGGSUMAR).
C
C     THE DIMENSION OF  A,  B, AND  Z  SHOULD BE  NM  BY  NM.
C     THE DIMENSION OF  ALFR,ALFI,BETA, AND  NORM  SHOULD BE  NM.
C     THE DIMENSION OF  ALFR1,ALFI1, AND  BETA1  SHOULD BE  NM.
C     THE DIMENSION OF  AHOLD  AND  BHOLD  SHOULD BE  NM  BY  NM.
C     HERE NM = 20.
C
      REAL A( 20, 20),B( 20, 20),Z( 20, 20),
     X        ALFR( 20),ALFI( 20),BETA( 20),NORM( 20),TCRIT( 1),
     X        EPSLON,RESDUL,DFL
      REAL AHOLD( 20, 20),BHOLD( 20, 20)
      REAL ALFR1( 20),ALFI1( 20),BETA1( 20)
      INTEGER IERR( 1),ERROR
      INTEGER MATCH( 1)
      CHARACTER*4 SAME,DIFF
      DATA IREAD1/1/,IREAD2/2/,IWRITE/6/
      DATA SAME,DIFF/'SAME','DIFF'/
C
      OPEN(IREAD1,FILE='FILE45')
      OPEN(IREAD2,FILE='FILE46')
      REWIND IREAD1
      REWIND IREAD2
C
      NM = 20
      LCOUNT = 0
      WRITE(*,1)
    1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
     XTATISTICS//1H ,95(1H-)/24H ORDER  QZIT(T)  QZIT(F) /
     X1H ,95(1H-)//
     X55H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX SYSTEM. //
     X95H UNDER 'QZIT(T)  QZIT(F)' ARE TWO NUMBERS AND A KEYWORD.  THE F
     XIRST NUMBER, AN INTEGER, IS THE /
     X95H ERROR FLAG RETURNED FROM  QZIT.  THE SECOND NUMBER IS THE MEAS
     XURE OF PERFORMANCE BASED UPON   /
     X95H THE RESIDUAL COMPUTED FOR THE  QZIT(T) (BOTH EIGENVALUES AND E
     XIGENVECTORS) PATH.  THE KEYWORD /
     X95H REPORTS THE DUPLICATION OF COMPUTED EIGENVALUES FROM THE  QZIT
     X(T)  AND  QZIT(F) (EIGENVALUES  /
     X95H ONLY) PATHS.  'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPL
     XICATES.  'DIFF' MEANS THAT FOR  /
     X95H AT LEAST ONE PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF
     X THE PAIR ARE NOT IDENTICAL.    )
      WRITE(*,2)
    2 FORMAT(81H0-1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ER
     XROR IN THE PATH HAS        /
     X58H PREVENTED THE COMPUTATION OF THE EIGENVALUES AND VECTORS. //
     X60H THE  QZIT(T)  PATH USES THE EISPACK CODES  QZHES-QZIT-QZVAL
     X,8H-QZVEC, /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RGG. /
     X61H THE  QZIT(F)  PATH USES THE EISPACK CODES  QZHES-QZIT-QZVAL, /
     X39H AS CALLED FROM DRIVER SUBROUTINE  RGG. )
      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    21H FOR THE REAL GENERAL/4X,30H GENERALIZED EIGENPROBLEM  A*X,
     X    15H= (LAMBDA)*B*X./
     X    24H0ORDER  QZIT(T)  QZIT(F) )
   10 CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,0)
      DO  510  ICALL = 1,2
         IF( ICALL .NE. 1 )  CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,1)
         GO TO  (70,100), ICALL
C
C     RGGWZ
C     INVOKED FROM DRIVER SUBROUTINE  RGG.
C
   70    ICT = 1
         CALL  RGG(NM,N,A,B,ALFR,ALFI,BETA,1,Z,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 500
         GO TO 190
C
C     RGGW
C     INVOKED FROM DRIVER SUBROUTINE  RGG.
C
  100    IF( IERR(1) .NE. 0 ) GO TO 510
         CALL  RGG(NM,N,A,B,ALFR1,ALFI1,BETA1,0,A,ERROR)
         IERR(1) = ERROR
         MATCH(1) = DIFF
         DO  101  I = 1,N
            IF( ALFR(I) .NE. ALFR1(I) .OR. ALFI(I) .NE. ALFI1(I) .OR.
     X          BETA(I) .NE. BETA1(I) ) GO TO 510
  101    CONTINUE
         MATCH(1) = SAME
         GO TO 510
C
  190    CALL  RMATIN(NM,N,A,B,AHOLD,BHOLD,1)
         CALL  RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL)
         DFL = 10 * N
         TCRIT(ICT) = RESDUL/EPSLON(DFL)
         GO TO 510
  500    TCRIT(ICT) = -1.0E0
  510 CONTINUE
C
      IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(*,5)
      LCOUNT = LCOUNT + 1
      WRITE(*,520) N,IERR(1),TCRIT(1),MATCH(1)
  520 FORMAT(I4,I5,F6.3,4X,A4)
      GO TO 10
      END
      SUBROUTINE RGGWZR(NM,N,A,B,ALFR,ALFI,BETA,Z,NORM,RESDUL)
C
      REAL NORM(N),ALFR(N),ALFI(N),BETA(N),A(NM,N),
     X       B(NM,N), Z(NM,N), TNORM, XR, XI, S, SUM,  SUMZ, SUMR, SUMI,
     X       SUMQ,RESDUL,NORMAB,SUM1,SUMA,SUMB,NORMA,NORMB,
     X       SUMR2,SUMI2,SUMR3,SUMI3,PYTHAG
      LOGICAL CPLX
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-B*Z*DIAG(W) WHERE A AND B ARE REAL GENERAL MATRICES, Z IS
C     A MATRIX WHICH CONTAINS THE EIGENVECTORS OF THE EIGENPROBLEM
C     A*Z - B*Z*DIAG(W), AND  W  STANDS FOR A VECTOR OF CORRESPONDING
C     EIGENVALUES OF THE EIGENPROBLEM OBTAINED FROM THE VECTORS  ALFR,
C     ALFI,  AND  BETA  BY THE CORRESPONDENCES
C                    W(J) = (ALFR(J) + I*ALFI(J)) / BETA(J) .
C     ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RGGWZR).
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 MATRICES  A  AND  B;
C
C        A(NM,N)  AND  B(NM,N)  ARE REAL GENERAL MATRICES;
C
C        Z(NM,N) IS A MATRIX WHICH CONTAINS THE REAL AND IMAGINARY PARTS
C           OF THE EIGENVECTORS OF THE SYSTEM.  THE I-TH COLUMN OF Z IS
C           A REAL EIGENVECTOR IF THE I-TH EIGENVALUE IS REAL.  IF THE
C           I-TH EIGENVALUE IS COMPLEX AND THE FIRST OF A CONJUGATE
C           PAIR, THE I-TH AND (I+1)-TH COLUMNS OF  Z  CONTAIN THE REAL
C           AND IMAGINARY PARTS OF ITS EIGENVECTORS;
C
C        ALFR(N), ALFI(N)  AND BETA(N)  ARE ARRAYS CONTAINING THE
C           COMPONENTS OF THE EIGENVALUES OF THE SYSTEM.
C
C     OUTPUT.
C
C        Z(NM,N) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
C           APPROXIMATE EIGENVECTORS OF THE SYSTEM.  THE EIGENVECTORS
C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST
C           ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
C           EIGENVECTOR DIVIDED BY N IS REAL AND POSITIVE;
C
C        NORM(N) IS AN ARRAY SUCH THAT FOR EACH K,
C
C                          !!BETA(K)*A*Z(K)-ALFA*B*Z(K)!!
C         NORM(K) = -----------------------------------------------
C                     !!Z(K)!!*(BETA(K)*!!A!!+!ALFA(K)!*!!B!!)
C
C         WHERE Z(K) IS THE K-TH EIGENVECTOR AND ALFA = (ALFR + I*ALFI);
C
C        RESDUL IS THE REAL NUMBER
C           !!BETA*A*Z-ALFA*B*Z!!/(!!Z!!*(BETA*!!A!!+!ALFA!*!!B!!))
C
C     ----------------------------------------------------------------
C
      NORMB = 0.0E0
      NORMA = 0.0E0
      RESDUL = 0.0E0
      CPLX = .FALSE.
C
      DO 20 I = 1,N
        SUMA = 0.0E0
        SUMB = 0.0E0
        DO 10 L = 1,N
          SUMA = SUMA + ABS(A(L,I))
   10     SUMB = SUMB + ABS(B(L,I))
        NORMA = AMAX1(NORMA,SUMA)
   20   NORMB = AMAX1(NORMB,SUMB)
C
      DO 160 I=1,N
         IF(CPLX) GO TO 80
         S = 0.0E0
         SUMZ = 0.0E0
         IF(ALFI(I) .NE. 0.0E0) GO TO 90
C        ..........THIS IS FOR REAL EIGENVALUES..........
         DO 40 L=1,N
            SUMQ = 0.0E0
            SUM1 = 0.0E0
            SUMZ = SUMZ + ABS(Z(L,I))
C
            DO 35 K=1,N
               SUMQ = SUMQ + B(L,K)*Z(K,I)
   35          SUM1 = SUM1 + A(L,K)*Z(K,I)
C
            SUM = -ALFR(I)*SUMQ
   40       S = S + ABS(SUM + SUM1*BETA(I))
C
         NORMAB = NORMA*BETA(I) + NORMB*ABS(ALFR(I))
         IF( NORMAB .EQ. 0.0E0 ) NORMAB = 1.0E0
         NORM(I) = SUMZ
         IF( SUMZ .EQ. 0.0E0 ) GO TO 150
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 50 L=1,N
            IF(ABS(Z(L,I)) .GE. NORM(I)/N) GO TO 60
   50       CONTINUE
C
   60    TNORM = SIGN(NORM(I),Z(L,I))
C
         DO 70 L=1,N
   70       Z(L,I) = Z(L,I)/TNORM
C
         NORM(I) = S/(NORM(I)*NORMAB)
   80    CPLX = .FALSE.
         GO TO 150
C        ..........THIS IS FOR COMPLEX EIGENVALUES..........
   90    DO 110 L = 1,N
            SUMR = 0.0E0
            SUMI = 0.0E0
            SUMR2 = 0.0E0
            SUMI2 = 0.0E0
            SUMZ = SUMZ + PYTHAG(Z(L,I),Z(L,I+1))
C
            DO 100 K=1,N
               SUMR = SUMR + B(L,K)*Z(K,I)
               SUMI = SUMI + B(L,K)*Z(K,I+1)
               SUMR2 = SUMR2 + A(L,K)*Z(K,I)
  100          SUMI2 = SUMI2 + A(L,K)*Z(K,I+1)
            SUMR3 = -ALFR(I)*SUMR + ALFI(I)*SUMI
            SUMI3 = -ALFI(I)*SUMR - ALFR(I)*SUMI
C
  110       S = S+PYTHAG(SUMR3+SUMR2*BETA(I),SUMI3+SUMI2*BETA(I))
            NORMAB = NORMA*BETA(I)+NORMB*PYTHAG(ALFR(I),ALFI(I))
            IF( NORMAB .EQ. 0.0E0 ) NORMAB = 1.0E0
C
         NORM(I) = SUMZ
         IF( SUMZ .EQ. 0.0E0 ) GO TO 150
C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE WILL
C                  ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I) LARGER
C                  THAN !!Z(I)!!/N..........
         DO 120 L=1,N
            IF(PYTHAG(Z(L,I),Z(L,I+1)) .GE. NORM(I)/N)
     1      GO TO 130
  120       CONTINUE
C
  130    XR = NORM(I)*Z(L,I)/PYTHAG(Z(L,I),Z(L,I+1))
         XI = NORM(I)*Z(L,I+1)/PYTHAG(Z(L,I),Z(L,I+1))
C
         DO 140 L= 1,N
            CALL CDIV(Z(L,I),Z(L,I+1),XR,XI,Z(L,I),Z(L,I+1))
  140    CONTINUE
C
         NORM(I) = S/(NORM(I)*NORMAB)
         NORM(I+1) = NORM(I)
         CPLX = .TRUE.
  150    RESDUL = AMAX1(NORM(I),RESDUL)
  160    CONTINUE
C
      RETURN
      END
      SUBROUTINE RMATIN(NM,N,A,B,AHOLD,BHOLD,INITIL)
C
C     THIS INPUT SUBROUTINE READS TWO REAL MATRICES  A  AND  B  FROM
C     SYSIN OF ORDER N.
C     TO GENERATE THE MATRICES INITIALLY,  INITIL  IS TO BE 0.
C     TO REGENERATE THE MATRICES FOR THE PURPOSE OF THE RESIDUAL
C     CALCULATION,  INITIL  IS TO BE  1.
C
C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RSGREADI).
C
      REAL A(NM,NM),B(NM,NM),AHOLD(NM,NM),BHOLD(NM,NM)
      INTEGER  IA( 20), IB( 20)
      DATA IREADA/1/,IREADB/2/,IWRITE/6/
C
      IF( INITIL .EQ. 1 )  GO TO  30
      READ(IREADA,5) N, M
    5 FORMAT(I6,6X,I6)
      IF( N .EQ. 0 )  GO TO  70
      IF (M .NE. 1) GO TO 16
      DO  10  I = 1,N
         READ(IREADA,17) (IA(J), J=I,N)
         DO   9  J = I,N
           A(I,J) = IA(J)
    9      A(J,I) =  A(I,J)
   10 CONTINUE
   11 READ(IREADB,5) N,M
      IF( M .NE. 1 ) GO TO 20
      DO 15 I = 1,N
         READ(IREADB,17) (IB(J), J=I,N)
         DO 14 J = I,N
            B(I,J)=IB(J)
   14       B(J,I)=B(I,J)
   15 CONTINUE
      GO TO 22
   16 DO  18  I = 1,N
         READ(IREADA,17) (IA(J), J=1,N)
   17    FORMAT(6I12)
         DO  18  J = 1,N
   18      A(I,J) = IA(J)
      GO TO 11
   20 DO 25 I = 1,N
         READ(IREADB,17) (IB(J),J=1,N)
         DO 25 J = 1,N
   25    B(I,J) = IB(J)
   22 DO  23  I = 1,N
         DO  23  J = 1,N
           BHOLD(I,J) = B(I,J)
   23      AHOLD(I,J) = A(I,J)
      RETURN
   30 DO  40  I = 1,N
         DO  40  J = 1,N
           B(I,J) = BHOLD(I,J)
   40      A(I,J) = AHOLD(I,J)
      RETURN
   70 WRITE(*,80)
   80 FORMAT(47H0END OF DATA FOR SUBROUTINE  RMATIN(RSGREADI). /1H1)
      STOP
      END
