C
C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF COMPLEX GENERAL
C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
C
C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(CGSUMARY).
C
C     THE DIMENSION OF  AR,AI,ZR,ZI,ASAVER,ASAVEI,RM1,  AND  RM2 SHOULD
C     BE  NM  BY  NM.
C     THE DIMENSION OF  WR,WI,WR1,WI1,SELECT,SLHOLD,INT,SCALE,ORTR,ORTI,
C     RV1  AND  RV2  SHOULD BE  NM.
C     THE DIMENSION OF  ARHOLD  AND  AIHOLD  SHOULD BE  NM  BY  NM.
C     HERE NM = 20.
C
      DOUBLE PRECISION AR( 20, 20),AI( 20, 20),ZR( 20, 20),ZI( 20, 20),
     X        ASAVER( 20, 20),ASAVEI( 20, 20),RM1( 20, 20),RM2( 20, 20),
     X        ORTR( 20),ORTI( 20),WR( 20),WI( 20),
     X        WR1( 20),WI1( 20),SCALE( 20),RV1( 20),RV2( 20),
     X        TCRIT( 8),EPSLON,RESDUL,DFL
      DOUBLE PRECISION ARHOLD( 20, 20),AIHOLD( 20, 20)
      INTEGER  INT( 20),IERR( 8),MATCH( 4),LOW,UPP,ERROR,SAME,DIFF
      LOGICAL  SELECT( 20),SLHOLD( 20)
      DATA IREAD1/1/,IREADC/5/,IWRITE/6/
      DATA  SAME,DIFF,LR,QR/4HSAME,4HDIFF,1HL,1HQ/
C
      OPEN(UNIT=IREAD1,FILE='FILE41')
      OPEN(UNIT=IREADC,FILE='FILE42')
      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-)/ 26X,9HBALANCING,15X,12HNO BALANCING,10X,
     X16HVECTORS SELECTED/ 18X,3(24H------------------------,1X)/
     X16H ORDER LOW UPP T,2X,
     X2(25H  COM-R2 COM-R   CINVIT     )/1H ,95(1H-)//
     X95H 'BALANCING' IS THE OPTION THAT EMPLOYS SUBROUTINE  CBAL  TO  B
     XALANCE THE MATRIX BEFORE THE    /
     X95H EIGENVALUES ARE COMPUTED AND  CBABK2  TO BACK TRANSFORM THE EI
     XGENVECTORS.                     //
     X95H 'VECTORS COMPUTED' IDENTIFIES BY POSITION, IN THE SET RETURNED
     X BY  COM-R, THOSE EIGENVALUES   /
     X95H FOR WHICH  CINVIT  COMPUTED THE ASSOCIATED EIGENVECTORS.  T  I
     XNDEXES AN EIGENVALUE FOR WHICH  /
     X66H THE EIGENVECTOR WAS COMPUTED AND  F  INDEXES A PASSED EIGENVAL
     XUE.    // 51H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX.   //
     X95H UNDER 'LOW' AND 'UPP' ARE INTEGERS INDICATING THE BOUNDARY IND
     XICES FOR THE BALANCED MATRIX.   //
     X95H UNDER 'T' IS THE LETTER L OR Q INDICATING THE USE OF ELEMENTAR
     XY OR UNITARY TRANSFORMATIONS.   /)
      WRITE(IWRITE,2)
    2 FORMAT(
     X61H PERFORMANCE REPORTED UNDER 'COM-R2 COM-R' IS THAT OF  COMLR2 ,
     X31H  AND  COMLR  ON LINE HEADED BY   /
     X60H 'L' AND THAT OF  COMQR2  AND  COMQR  ON LINE HEADED BY 'Q'.//
     X95H UNDER 'COM-R2 COM-R' ARE TWO NUMBERS AND A KEYWORD.  THE FIRST
     X NUMBER, AN INTEGER, IS THE     /
     X95H ABSOLUTE SUM OF THE ERROR FLAGS RETURNED SEPARATELY FROM  COM-
     XR2  AND  COM-R.  THE SECOND     /
     X95H NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE RESIDUAL C
     XOMPUTED FOR THE  COM-R2  PATH.  /
     X95H THE KEYWORD REPORTS THE DUPLICATION OF THE COMPUTED EIGENVALUE
     XS FROM  COM-R2  AND  COM-R.     /
     X95H 'SAME' MEANS THAT THE EIGENVALUES ARE EXACT DUPLICATES.  'DIFF
     X' MEANS THAT FOR AT LEAST ONE   /
     X95H PAIR OF CORRESPONDING EIGENVALUES, THE MEMBERS OF THE PAIR ARE
     X NOT IDENTICAL.                 //
     X95H UNDER 'CINVIT' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,
     X IS THE ABSOLUTE SUM OF THE     /
     X62H ERROR FLAGS RETURNED FROM THE PATH.  THE SECOND NUMBER IS THE,
     X29H MEASURE OF PERFORMANCE BASED    )
      WRITE(IWRITE,3)
    3 FORMAT(
     X41H UPON THE RESIDUAL COMPUTED FOR THE PATH.//
     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.//
     X45H0FOR REDUCTIONS BY ELEMENTARY TRANSFORMATIONS/
     X95H THE  COMLR2   PATH WITH    BALANCING USES THE EISPACK CODES  C
     XBAL  -COMHES-COMLR2-CBABK2.     /
     X95H THE   COMLR   PATH WITH    BALANCING USES THE EISPACK CODES  C
     XBAL  -COMHES-COMLR .            /
     X60H THE  CINVIT   PATH WITH    BALANCING USES THE EISPACK CODES ,
     X44H  CBAL  -COMHES-COMLR -CINVIT-COMBAK-CBABK2. /
     X95H THE  COMLR2   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XOMHES-COMLR2.                   /
     X95H THE   COMLR   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XOMHES-COMLR .                   /
     X95H THE  CINVIT   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XOMHES-COMLR -CINVIT-COMBAK.     )
      WRITE(IWRITE,4)
    4 FORMAT(42H0FOR REDUCTIONS BY UNITARY TRANSFORMATIONS/
     X95H THE  COMQR2   PATH WITH    BALANCING USES THE EISPACK CODES  C
     XBAL  -CORTH -COMQR2-CBABK2,           /
     X38H AS CALLED FROM DRIVER SUBROUTINE  CG. /
     X95H THE   COMQR   PATH WITH    BALANCING USES THE EISPACK CODES  C
     XBAL  -CORTH -COMQR ,            /
     X38H AS CALLED FROM DRIVER SUBROUTINE  CG. /
     X60H THE  CINVIT   PATH WITH    BALANCING USES THE EISPACK CODES ,
     X44H  CBAL  -CORTH -COMQR -CINVIT-CORTB -CBABK2. /
     X95H THE  COMQR2   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XORTH -COMQR2.                   /
     X95H THE   COMQR   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XORTH -COMQR .                   /
     X95H THE  CINVIT   PATH WITH NO BALANCING USES THE EISPACK CODES  C
     XORTH -COMQR -CINVIT-CORTB .     )
      WRITE(IWRITE,15)
   15 FORMAT(1X,21HD.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 GENERAL MATRICES.      //
     X    26X,9HBALANCING,15X,12HNO BALANCING,10X,16HVECTORS SELECTED/
     X    18X,3(24H------------------------,1X)/
     X    16H ORDER LOW UPP T,2X,2(25H  COM-R2 COM-R   CINVIT       ))
   10 CALL CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,0)
      READ(IREADC,50) MM,(SLHOLD(I),I=1,N)
   50 FORMAT(I4/(72L1))
C
C     MM  AND  SELECT  ARE READ FROM SYSIN AFTER THE MATRIX IS
C     GENERATED.  MM  SPECIFIES TO  INVIT  THE MAXIMUM -UMBER OF
C     EIGENVECTORS THAT WILL BE "OMPUTED.  SELECT  CONTAINS INFORMATION
C     ABOUT WHICH EIGENVECTORS ARE DESIRED.
C
      DO  300  ICALL = 1,12
         ICALL2 = MOD(ICALL,3)
         IF( ICALL2 .NE. 0 ) GO TO 80
         DO  70  I = 1,N
   70    SELECT(I) = SLHOLD(I)
   80    IF( ICALL .NE. 1 )  CALL  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,1)
         GO TO  (90,100,110,120,130,140,160,165,170,175,180,185), ICALL
C
C     IF  COM-R  PATH IS TAKEN THEN  COM-R2  PATH MUST ALSO BE TAKEN
C     IN ORDER THAT COMPARISON OF EIGENVALUES BE POSSIBLE.
C
C
C     CGEWZ
C
   90    ICT = 1
         CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
         CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
         CALL  COMLR2(NM,N,LOW,UPP,INT,AR,AI,WR,WI,ZR,ZI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 220
         DO  91  I = 1,N
            WR1(I) = WR(I)
   91       WI1(I) = WI(I)
         CALL  CBABK2(NM,N,LOW,UPP,SCALE,N,ZR,ZI)
         GO TO  190
C
C     CGEW
C
  100    IF( IERR(1) .NE. 0 ) GO TO 300
         CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
         CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
         CALL  COMLR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
         IERR(1) = ERROR
         MATCH(1) = DIFF
         DO  101  I = 1,N
            IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  101    CONTINUE
         MATCH(1) = SAME
         GO TO  300
C
C     CGEW1Z
C
  110    ICT = 2
         CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
         CALL  COMHES(NM,N,LOW,UPP,AR,AI,INT)
         DO  112  I = 1,N
            DO  112  J = 1,N
              ASAVER(I,J) = AR(I,J)
  112         ASAVEI(I,J) = AI(I,J)
         CALL  COMLR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 )  GO TO  220
         CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
     X                ERROR,RM1,RM2,RV1,RV2)
         IERR(ICT) = IABS(ERROR)
         CALL  COMBAK(NM,LOW,UPP,ASAVER,ASAVEI,INT,M,ZR,ZI)
         CALL  CBABK2(NM,N,LOW,UPP,SCALE,M,ZR,ZI)
         GO TO  190
C
C     CGNEWZ
C
  120    ICT = 3
         CALL  COMHES(NM,N,1,N,AR,AI,INT)
         CALL  COMLR2(NM,N,1,N,INT,AR,AI,WR,WI,ZR,ZI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 220
         DO  121  I = 1,N
            WR1(I) = WR(I)
  121       WI1(I) = WI(I)
         GO TO  190
C
C     CGNEW
C
  130    IF( IERR(3) .NE. 0 ) GO TO 300
         CALL  COMHES(NM,N,1,N,AR,AI,INT)
         CALL  COMLR(NM,N,1,N,AR,AI,WR,WI,ERROR)
         IERR(3) = ERROR
         MATCH(2) = DIFF
         DO  131  I = 1,N
            IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  131    CONTINUE
         MATCH(2) = SAME
         GO TO  300
C
C     CGNEW1Z
C
  140    ICT = 4
         CALL  COMHES(NM,N,1,N,AR,AI,INT)
         DO  142  I = 1,N
            DO  142  J = 1,N
              ASAVER(I,J) = AR(I,J)
  142         ASAVEI(I,J) = AI(I,J)
         CALL  COMLR(NM,N,1,N,AR,AI,WR,WI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 )  GO TO  220
         CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
     X                ERROR,RM1,RM2,RV1,RV2)
         IERR(ICT) = IABS(ERROR)
         CALL  COMBAK(NM,1,N,ASAVER,ASAVEI,INT,M,ZR,ZI)
         GO TO 190
C
C     CGUWZ
C     INVOKED FROM DRIVER SUBROUTINE  CG.
C
  160    ICT = 5
         CALL  CG(NM,N,AR,AI,WR,WI,1,ZR,ZI,SCALE,ORTR,ORTI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 220
         GO TO  190
C
C     CGUW
C     INVOKED FROM DRIVER SUBROUTINE  CG.
C
  165    IF( IERR(5) .NE. 0 ) GO TO 300
         CALL  CG(NM,N,AR,AI,WR1,WI1,0,AR,AI,SCALE,ORTR,ORTI,ERROR)
         IERR(5) = ERROR
         MATCH(3) = DIFF
         DO  166  I = 1,N
            IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  166    CONTINUE
         MATCH(3) = SAME
         GO TO  300
C
C     CGUW1Z
C
  170    ICT = 6
         CALL  CBAL(NM,N,AR,AI,LOW,UPP,SCALE)
         CALL  CORTH(NM,N,LOW,UPP,AR,AI,ORTR,ORTI)
         DO  172  I = 1,N
            DO  172  J = 1,N
              ASAVER(I,J) = AR(I,J)
  172         ASAVEI(I,J) = AI(I,J)
         CALL  COMQR(NM,N,LOW,UPP,AR,AI,WR,WI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 )  GO TO  220
         CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
     X                        ERROR,RM1,RM2,RV1,RV2)
         IERR(ICT) = IABS(ERROR)
         CALL  CORTB(NM,LOW,UPP,ASAVER,ASAVEI,ORTR,ORTI,M,ZR,ZI)
         CALL  CBABK2(NM,N,LOW,UPP,SCALE,M,ZR,ZI)
         GO TO  190
C
C     CGNUWZ
C
  175    ICT = 7
         CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
         CALL  COMQR2(NM,N,1,N,ORTR,ORTI,AR,AI,WR,WI,ZR,ZI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 ) GO TO 220
         DO  176  I = 1,N
            WR1(I) = WR(I)
  176       WI1(I) = WI(I)
         GO TO  190
C
C     CGNUW
C
  180    IF( IERR(7) .NE. 0 ) GO TO 300
         CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
         CALL  COMQR(NM,N,1,N,AR,AI,WR,WI,ERROR)
         IERR(7) = ERROR
         MATCH(4) = DIFF
         DO  181  I = 1,N
            IF( WR(I) .NE. WR1(I) .OR. WI(I) .NE. WI1(I) )  GO TO  300
  181    CONTINUE
         MATCH(4) = SAME
         GO TO  300
C
C     CGNUW1Z
C
  185    ICT = 8
         CALL  CORTH(NM,N,1,N,AR,AI,ORTR,ORTI)
         DO  187  I = 1,N
            DO  187  J = 1,N
              ASAVER(I,J) = AR(I,J)
  187         ASAVEI(I,J) = AI(I,J)
         CALL  COMQR(NM,N,1,N,AR,AI,WR,WI,ERROR)
         IERR(ICT) = ERROR
         IF( ERROR .NE. 0 )  GO TO  220
         CALL  CINVIT(NM,N,ASAVER,ASAVEI,WR,WI,SELECT,MM,M,ZR,ZI,
     X                ERROR,RM1,RM2,RV1,RV2)
         IERR(ICT) = IABS(ERROR)
         CALL  CORTB(NM,1,N,ASAVER,ASAVEI,ORTR,ORTI,M,ZR,ZI)
C
  190    CALL  CMATIN(NM,N,AR,AI,ARHOLD,AIHOLD,1)
         IF( ICALL2 .EQ. 0 ) CALL CGW1ZR(NM,N,AR,AI,WR,WI,SELECT,
     X                                   M,ZR,ZI,RV1,RESDUL)
         IF( ICALL2 .EQ. 1 ) CALL CGWZR(NM,N,AR,AI,WR,WI,ZR,ZI,
     X                                  RV1,RESDUL)
         DFL = 10 * N
         TCRIT(ICT) = RESDUL/EPSLON(DFL)
         GO TO 300
  220    TCRIT(ICT) = -1.0D0
  300 CONTINUE
C
      IF( MOD(LCOUNT,16) .EQ. 0 ) WRITE(IWRITE,5)
      LCOUNT = LCOUNT + 1
      WRITE(IWRITE,520) N,LOW,UPP,LR,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
     X              IERR(2*I),TCRIT(2*I),I=1,2),(SELECT(I),I=1,N)
  520 FORMAT(I4,1X,2I4,2X,A1,2(I4,F6.3,1X,A4,I4,F6.3),
     X       4X,21L1/(72X,20L1))
      WRITE(IWRITE,530) QR,(IERR(2*I-1),TCRIT(2*I-1),MATCH(I),
     X              IERR(2*I),TCRIT(2*I),I=3,4)
  530 FORMAT(15X,A1,2(I4,F6.3,1X,A4,I4,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
      DOUBLE PRECISION 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 CGWZR(NM,N,AR,AI,WR,WI,ZR,ZI,NORM,RESDUL)
C
      DOUBLE PRECISION NORM(N),WR(N),WI(N),AR(NM,N),AI(NM,N),
     X       ZR(NM,N), ZI(NM,N), NORMA, XR, XI, S, SUMA, SUMZ,
     X       SUMR, SUMI, RESDUL, PYTHAG
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-Z*DIAG(W)  WHERE  A  IS A COMPLEX GENERAL MATRIX,  Z  IS A
C     MATRIX WHOSE COLUMNS ARE APPROXIMATE EIGENVECTORS OF  A,
C     AND  W  IS A VECTOR WHOSE ELEMENTS ARE APPROXIMATE EIGENVALUES
C     CORRESPONDING TO THE VECTORS IN  Z.  ALL NORMS USED IN THE
C     COMMENTS BELOW ARE 1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(CGWZR)
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        AR(NM,N), AI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
C           IMAGINARY PARTS OF THE ELEMENTS OF  A;
C
C        ZR(NM,N), ZI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
C           IMAGINARY PARTS OF THE ELEMENTS OF  Z;
C
C        WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
C           PARTS OF  W.
C
C     OUTPUT.
C
C        ZR(NM,N), ZI(NM,N) 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
      RESDUL = 0.0D0
      NORMA = 0.0D0
C
      DO 20 I=1,N
         SUMA = 0.0D0
C
         DO 10 L=1,N
   10       SUMA = SUMA + PYTHAG(AR(L,I),AI(L,I))
C
   20    NORMA = DMAX1(NORMA,SUMA)
C
      IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
C
      DO 80 I=1,N
         SUMZ = 0.0D0
         S = 0.0D0
C
         DO 40 L=1,N
            SUMZ = SUMZ + PYTHAG(ZR(L,I),ZI(L,I))
            SUMR = -WR(I)*ZR(L,I) + WI(I)*ZI(L,I)
            SUMI = -WR(I)*ZI(L,I) - WI(I)*ZR(L,I)
C
            DO 30 K=1,N
               SUMR = SUMR + AR(L,K)*ZR(K,I) - AI(L,K)*ZI(K,I)
   30          SUMI = SUMI + AR(L,K)*ZI(K,I) + AI(L,K)*ZR(K,I)
C
   40       S = S + PYTHAG(SUMR,SUMI)
C
         NORM(I) = SUMZ
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(PYTHAG(ZR(L,I),ZI(L,I)) .GE. NORM(I)/N) GO
     1      TO 60
   50        CONTINUE
C
   60    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 70 L=1,N
            CALL CDIV(ZR(L,I),ZI(L,I),XR,XI,ZR(L,I),ZI(L,I))
   70    CONTINUE
C
         NORM(I) = S/(NORM(I)*NORMA)
         RESDUL = DMAX1(NORM(I),RESDUL)
   80    CONTINUE
C
      RETURN
      END
      SUBROUTINE CGW1ZR(NM,N,AR,AI,WR,WI,SELECT,M,ZR,ZI,NORM,RESDUL)
C
      DOUBLE PRECISION NORM(N),WR(N),WI(N),AR(NM,N),AI(NM,N),
     X       ZR(NM,M), ZI(NM,M), NORMA, XR, XI, S, SUMA, SUMZ,
     X       SUMR, SUMI, RESDUL, PYTHAG
      LOGICAL SELECT(N)
C
C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
C     A*Z-Z*DIAG(W)   WHERE  A  IS A COMPLEX GENERAL MATRIX,  Z  IS A
C     MATRIX WHOSE COLUMNS ARE APPROXIMATE EIGENVECTORS OF  A,
C     AND  W  IS A VECTOR WHOSE ELEMENTS ARE APPROXIMATE EIGENVALUES
C     CORRESPONDING TO THE VECTORS IN  Z.  ALL NORMS USED IN THE
C     COMMENTS BELOW ARE 1-NORMS.
C
C     THIS SUBROUTINE IS CATALOGUED AS  EISPDRV4(CGW1ZR).
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        AR(NM,N), AI(NM,N) ARE ARRAYS CONTAINING THE REAL AND
C           IMAGINARY PARTS OF THE ELEMENTS OF  A;
C
C        ZR(NM,M), ZI(NM,M) ARE ARRAYS CONTAINING THE REAL AND
C           IMAGINARY PARTS OF THE ELEMENTS OF  Z;
C
C        WR(N), WI(N) ARE ARRAYS CONTAINING THE REAL AND IMAGINARY
C           PARTS OF  W;
C
C        SELECT(N)  IS A BOOLEAN ARRAY WHICH SPECIFIES THE EIGENVALUES
C           IN THE VECTOR  W  THAT ARE ASSOCIATED WITH EIGENVECTORS
C           IN THE MATRIX  Z.  THE I-TH VALUE IN  W  IS ASSOCIATED WITH
C           A VECTOR IN  Z  IF SELECT(I) = .TRUE..  FURTHERMORE IF
C           SELECT(I) IS THE J-TH  .TRUE.  ELEMENT, THEN THE I-TH
C           ELEMENT OF  W  IS ASSOCIATED WITH THE J-TH VECTOR OF  Z;
C
C        M  IS THE NUMBER OF ELEMENTS IN SELECT WHOSE VALUE IS  .TRUE..
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 WHENEVER SELECT(K) =  .TRUE.
C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
C           WHERE  Z(K)  IS THE EIGENVECTOR CORRESPONDING TO THE
C           K-TH EIGENVALUE (WR(K),WI(K));
C
C        RESDUL  IS A REAL NUMBER DEFINED BY THE FOLLOWING.
C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
C
C     ----------------------------------------------------------------
C
      RESDUL = 0.0D0
      IF( M .EQ. 0 ) RETURN
      NORMA = 0.0D0
      II=0
C
      DO 20 I=1,N
         SUMA = 0.0D0
C
         DO 10 L=1,N
   10       SUMA = SUMA + PYTHAG(AR(L,I),AI(L,I))
C
         NORMA = DMAX1(NORMA,SUMA)
         IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
   20    CONTINUE
C
      IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
C
      DO 80 I=1,N
         SUMZ = 0.0D0
         S = 0.0D0
         IF( .NOT. SELECT(I) ) GO TO 80
         II = II+1
C
         DO 40 L=1,N
            SUMZ = SUMZ + PYTHAG(ZR(L,II),ZI(L,II))
            SUMR = -WR(I)*ZR(L,II) + WI(I)*ZI(L,II)
            SUMI = -WR(I)*ZI(L,II) - WI(I)*ZR(L,II)
C
            DO 30 K=1,N
               SUMR = SUMR + AR(L,K)*ZR(K,II) - AI(L,K)*ZI(K,II)
   30          SUMI = SUMI + AR(L,K)*ZI(K,II) + AI(L,K)*ZR(K,II)
C
   40       S = S + PYTHAG(SUMR,SUMI)
C
         IF(SUMZ .NE. 0.0D0) GO TO 45
         NORM(I) = -1
         GO TO 80
   45    NORM(I) = SUMZ
C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR  Z(II)
C                  LARGER THAN !!Z(II)!!/N..........
         DO 50 L=1,N
            IF(PYTHAG(ZR(L,II),ZI(L,II)) .GE. NORM(I)/N)
     1      GO TO 60
   50        CONTINUE
C
   60    XR = NORM(I)*ZR(L,II)/PYTHAG(ZR(L,II),ZI(L,II))
         XI = NORM(I)*ZI(L,II)/PYTHAG(ZR(L,II),ZI(L,II))
C
         DO 70 L=1,N
            CALL CDIV(ZR(L,II),ZI(L,II),XR,XI,ZR(L,II),ZI(L,II))
   70    CONTINUE
C
         NORM(I) = S/(NORM(I)*NORMA)
         RESDUL = DMAX1(NORM(I),RESDUL)
   80    CONTINUE
C
      RETURN
      END