C	WESTERN MICHIGAN UNIVERSITY
C	EVEV.F4 (FILE NAME ON LIBRARY DECTAPE)
C	EVEV, 2.4.1 (CALLING NAME, SUBLST NO.)
C	EIGENVALUES, EIGENVECTORS - COMPLEX INPUT AND OUTPUT
C	 ALLOWED.  THIS PROGRAM IS A COMBINATION OF A PROGRAM GIVEN
C	 BY WAYNE STATE UNIVERSITY WITH SOME REVISIONAL AND
C	 ADDITIONAL PROGRAMMING BY B. GRANET AND B. HOUCHARD.
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	FORWMU PROGS. USED:  TTYPTY, PRINTS, DEVCHG, EXISTS, DEVICE
C	INTERNAL SUBR. USED:  ALLMAT, CHECK
C	APLIB PROGS. USED:  IOB, GETFOR
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C---------------EIGENVALUES ARE STORED IN EIGVAL(35), A(35,35) AT FIRST
C--------------- CONTAINS USER'S INPUT MATRIX AND LATER CONTAINS
C---------------EIGENVECTORS.  B(35,35) CONTAINS A COPY OF USER'S
C--------------- INPUT MATRIX AND IS USED IN SUBR. CHECK
      COMPLEX EIGVAL(35),A(35,35),B(35,35)
      DIMENSION ID(16),IFMT(96)
C---------------COMMON /IOBLK/ IS SHARED BY SUBR. IOB
C****AM,2.4.1,#1,WG,13-DEC-77
	COMMON/IOBLK/IDLG,INT,INP,IRP,IDEV,IDEVA,ICODE,IB,NAMI(2)
C****END,MAIN PROG.,STAT.5007-8
      DATA IFMT(1)/'(20G)'/
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICODE)
      IA=35
      IDLG=-1
      INT=5
      IRP=2
      INP=1
	WRITE(IDLG,5007)
5007	FORMAT(1X,'WMU EIGENVALUE AND EIGENVECTORS',/)
C	CALL USAGE('EVEV')
C---------------1 MEANS OUTPUT? PRINTS, 0 - INPUT? PRINTS.
C--------------- THRU COMMON /IOBLK/ IDLG, INT, IRP, IDEV, ICODE
C--------------- ARE INPUT;  IB, NAMI, ARE RETURNED.
	CALL IOB(1)
5006	CALL IOB(0)
200   WRITE(IDLG,100)
100   FORMAT(' ','ENTER OUTPUT IDENTIFICATION'/)
      READ(INT,101) (ID(I),I=1,16)
101   FORMAT(16A5)
15    WRITE(IDLG,1)
1     FORMAT(1X,'ENTER ORDER OF MATRIX.'/)
      READ(INT,2) N
2     FORMAT(I)
      IF(N.GE.1.AND.N.LE.35) GO TO 195
      CALL DEVICE (INT)
      GO TO 15
C---------------IFMT, ISTD ARE RETURNED.  OTHER ARGS. ARE
C--------------- INPUT.  96=NO. OF OBJ. TIME FORMAT WORDS (6 LINES)
C--------------- 4 MEANS WE ARE NOT LIMITING FORMAT.
195   CALL GETFOR(IDLG,INT,IFMT,ISTD,96,4)
      IF(ISTD.EQ.1)IFMT(1)='(20G)'
5009  IF(IDEV.NE.'TTY') GO TO 5001
104   WRITE(IDLG,3)
3     FORMAT(1X,'ENTER DATA,REAL PART FIRST THEN THE IMAGINARY;'/1X,
     1'(AT MOST 10 SETS PER LINE)'/,
     1	' NOTE: USE COMMA(,) INSTEAD OF SEMICOLON(;) TO SEPARATE PAIRS'/
     1	' OF COMPLEX NUMBERS',/)
5003  DO 4 I=1,N
4     READ(INP,IFMT,END=5004)(A(I,J),J=1,N)
      DO 110 I=1,N
      DO 110 J=1,N
110   B(I,J)=A(I,J)
      WRITE (IRP ,111)
111   FORMAT('1')
      WRITE(IRP,101)(ID(I),I=1,16)
      CALL ALLMAT(A,EIGVAL,N,IA,NCAL)
      WRITE (IRP ,9)NCAL
9     FORMAT(1H-,'NUMBER OF EIGENVALUES CALCULATED IS ',I10/)
      WRITE (IRP ,10)
10    FORMAT(1X,'EIGENVALUES  ',7X,'REAL',10X,'IMAGINARY'/)
      DO 11 I=1,NCAL
11    WRITE(IRP ,12)EIGVAL(I)
12    FORMAT(15X,2G)
      WRITE(IRP ,13)
13    FORMAT(1H-,'COMPLEX EIGENVECTORS-NORMALIZED TO LENGTH =1.'/)
      DO 50 L=1,NCAL
      SUM=0.0
      DO 51 I=1,N
51    SUM=SUM+REAL(A(I,L))**2+AIMAG(A(I,L))**2
      SUM=SQRT(SUM)
      DO 52 I=1,N
52    A(I,L)=A(I,L)/SUM
50    CONTINUE
C     COLUMNS OF A ARE EIGENVECTORS OF ORIGINAL A
C     SEE 3RD AND 4TH LINES FROM END OF PROGRAM
C     AND LINE WITH STATEMENT NUMBER 38.
      DO 14 J=1,NCAL
14    WRITE(IRP ,109)J,(A(I,J),I=1,N)
109   FORMAT('-','EIGENVECTOR ',I2/(' ',4G))
      CALL CHECK  (B,A,NCAL,N,EIGVAL,IRP)
      GO TO 5006
5001  WRITE(IDLG,5002)
5002  FORMAT(1X,'DATA BEING PROCESSED'/)
      GO TO 5003
5004  WRITE(IDLG,5005)
5005  FORMAT(1X,'YOU PROCESSED THE LAST MATRIX OR YOU ARE EXPECTING'/1X,
     1'MORE DATA WHEN AN END OF FILE MARK APPEARED.'/)
      GO TO 5006
      END
C---------------A, N ARE INPUT;  LAMBDA, IA, NCAL, A ARE OUTPUT. 
C--------------- ON OUTPUT THE FIRST NCAL COLS. OF A CONTAIN
C--------------- EIGENVECTORS OF THE INPUT A.
      SUBROUTINE ALLMAT (A,LAMBDA,M,IA,NCAL)
C
C     PROG. AUTHOURS JOHN RINZEL, R.L.FUNDERLIC, UNION CARBIDE CORP.
C     NUCLEAR DIVISION, CENTRAL DATA PROCESSING FACILITY,
C     OAK RIDGE TENNESSEE
C
      COMPLEX A(IA,1),H(35,35), HL(35,35),LAMBDA(1), VECT (35),
     1MULT(35), SHIFT(3), TEMP, SIN, COS, TEMP1, TEMP2
      LOGICAL INTH(35), TWICE
      INTEGER INT(35), R, RP1, RP2
      N=M
      NCAL = N
      IF (N .NE. 1) GO TO 1
      LAMBDA (1) = A(1,1)
      A(1,1) = 1.
      GO TO 57
1     ICOUNT = 0
      SHIFT (1) = 0.
      IF (N.NE. 2) GO TO 4
2     TEMP = (A(1,1) + A(2,2) + CSQRT((A(1,1) + A(2,2))**2-
     14. * (A(2,2) * A(1,1) - A(2,1) * A(1,2))))/2.
      IF (REAL(TEMP) .NE. 0..OR.AIMAG(TEMP).NE.0.) GO TO 3
      LAMBDA (M) = SHIFT (1)
      LAMBDA (M-1) = A(1,1) + A(2,2) + SHIFT (1)
      GO TO 37
3     LAMBDA (M) = TEMP + SHIFT (1)
      LAMBDA (M-1) = (A(2,2) *A(1,1) - A(2,1) * A(1,2))/ (LAMBDA (M)
     1-SHIFT (1)) + SHIFT (1)
      GO TO 37
C     REDUCE MATRIX A TO HESSENBERG FORM
C
4     NM2 = N-2
      DO 15 R = 1,NM2
      RP1 = R+ 1
      RP2 = R + 2
      ABIG= 0.
      INT (R) = RP1
      DO 5 I = RP1, N
      ABSSQ = REAL (A(I,R))**2 + AIMAG(A(I,R))**2
      IF (ABSSQ .LE. ABIG) GO TO 5
      INT (R) = I
      ABIG = ABSSQ
5     CONTINUE
      INTER = INT (R)
      IF (ABIG .EQ. 0.) GO TO 15
      IF (INTER .EQ. RP1) GO TO 8
      DO 6 I = R,N
      TEMP = A(RP1,I)
      A(RP1,I) = A(INTER,I)
6     A(INTER,I) = TEMP
      DO 7 I = 1,N
      TEMP = A(I,RP1)
      A(I,RP1) = A(I,INTER)
7     A(I,INTER) = TEMP
8     DO 9 I = RP2,N
      MULT(I) = A(I,R)/A(RP1,R)
9     A(I,R) = MULT (I)
      DO 11 I = 1,RP1
      TEMP = 0.
      DO 10 J = RP2,N
10    TEMP=TEMP+A(I,J)*MULT(J)
11    A(I,RP1)=A(I,RP1)+TEMP
      DO 13 I = RP2,N
      TEMP = 0.
      DO 12 J = RP2,N
12    TEMP = TEMP + A(I,J) * MULT(J)
13    A(I,RP1) = A(I,RP1) + TEMP - MULT(I) * A(RP1,RP1)
      DO 14 I = RP2,N
      DO 14 J = RP2,N
14    A(I,J) = A(I,J) - MULT(I) * A(RP1,J)
15    CONTINUE
C
C     CALCULATE EPSILON
C
      EPS = 0.
      DO 16 I = 1,N
16    EPS = EPS + CABS (A(1,I))
      DO 18 I = 2,N
      SUM = 0.
      IM1 = I-1
      DO 17 J = IM1,N
17    SUM=SUM + CABS(A(I,J))
18    IF(SUM .GT. EPS) EPS = SUM
      EPS = SQRT(FLOAT(N)) * EPS * 1.E-12
      IF (EPS .EQ. 0.) EPS = 1.E-12
      DO 19 I = 1,N
      DO 19 J = 1,N
19    H(I,J) = A(I,J)
20    IF(N .NE. 1) GO TO 21
      LAMBDA (M) = A(1,1) + SHIFT (1)
      GO TO 37
21    IF(N .EQ. 2)  GO TO 2
22    MN1 = M - N+1
      IF(REAL(A(N,N)).NE.0..OR.AIMAG(A(N,N)).NE.0.)
     1IF(ABS(REAL(A(N,N-1)/A(N,N)))+ABS(AIMAG(A(N,N-1)/A(N,N)))-1.E-9)
     224,24,23
23    IF(ABS(REAL(A(N,N-1)))+ABS(AIMAG(A(N,N-1))).GE.EPS) GO TO 25
24    LAMBDA (MN1) = A(N,N) + SHIFT (1)
      ICOUNT = 0
      N = N - 1
      GO TO 21
C
C     DETERMINE SHIFT
C
25    SHIFT(2) = (A(N-1,N-1)+A(N,N)+CSQRT((A(N-1,N-1)+A(N,N))**2
     1-4.*(A(N,N)*A(N-1,N-1)-A(N,N-1)*A(N-1,N))))/2
      IF(REAL(SHIFT(2)).NE.0.OR.AIMAG(SHIFT(2)).NE.0.)GO TO 26
      SHIFT (3) = A(N-1,N-1) + A(N,N)
      GO TO 27
26    SHIFT(3) = (A(N,N)*A(N-1,N-1)-A(N,N-1)*A(N-1,N))/SHIFT(2)
27    IF(CABS(SHIFT(2)-A(N,N)).LT.CABS(SHIFT(3)-A(N,N)))GO TO 28
      INDEX = 3
      GO TO 29
28    INDEX = 2
29    IF(CABS(A(N-1,N-2)).GE.EPS) GO TO 30
      LAMBDA (MN1) = SHIFT (2) + SHIFT (1)
      LAMBDA(MN1+1)=SHIFT(3)+SHIFT(1)
      ICOUNT = 0
      N = N - 2
      GO TO 20
30    SHIFT (1) = SHIFT (1) + SHIFT (INDEX)
      DO 31 I = 1,N
31    A(I,I)=A(I,I)-SHIFT(INDEX)
C     PERFORM GIVEN ROTATIONS, OR ITERATES
C
      IF (ICOUNT.LE.10) GO TO 32
      NCAL=M-N
      GO TO 37
32    NM1 = N - 1
      TEMP1 = A(1,1)
      TEMP2 = A(2,1)
      DO 36 R = 1,NM1
      RP1 = R + 1
      RHO = SQRT(REAL(TEMP1)**2+AIMAG(TEMP1)**2+
     1REAL(TEMP2)**2+AIMAG(TEMP2)**2)
      IF(RHO.NE.0)GO TO 60
      TEMP1=A(RP1,RP1)
      TEMP2=A(R+2,RP1)
      GO TO 36
60    COS=TEMP1/RHO
      SIN = TEMP2/RHO
      INDEX = MAX0 (R-1,1)
      DO 33 I = INDEX,N
      TEMP = CONJG(COS)*A(R,I)+CONJG(SIN)*A(RP1,I)
      A(RP1,I) = -SIN*A(R,I) + COS*A(RP1,I)
33    A(R,I) = TEMP
      TEMP1 = A(RP1,RP1)
      TEMP2 = A(R+2, R+1)
      DO 34 I = 1,R
      TEMP = COS* A(I,R) + SIN*A(I,RP1)
      A(I,RP1) = -CONJG(SIN)*A(I,R) + CONJG(COS)*A(I,RP1)
34    A(I,R)= TEMP
      INDEX = MIN0 (R+2,N)
      DO 35 I = RP1,INDEX
      A(I,R) = SIN * A(I,RP1)
35    A(I,RP1) = CONJG(COS) * A(I,RP1)
36    CONTINUE
      ICOUNT = ICOUNT + 1
      GO TO 22
C
C     CALCULATE VECTORS
C
37    IF(NCAL.EQ.0) GO TO 57
      N = M
      NM1 = N- 1
      IF (N.NE.2) GO TO 38
      EPS = AMAX1(CABS(LAMBDA(1)),CABS(LAMBDA(2)))*1.E-8
      IF(EPS.EQ.0) EPS = 1.E-12
      H(1,1) = A(1,1)
      H(1,2) = A(1,2)
      H(2,1) = A(2,1)
      H(2,2) = A(2,2)
38    DO 56 L = 1,NCAL
      DO 40 I = 1,N
      DO 39 J = 1,N
39    HL(I,J) = H(I,J)
40    HL(I,I) = HL(I,I) - LAMBDA(L)
      DO 44 I = 1,NM1
      MULT (I) = 0.
      INTH(I) = .FALSE.
      IP1 = I + 1
      IF(CABS(HL(I+1,I)).LE.CABS(HL(I,I)))GO TO 42
      INTH (I) = .TRUE.
      DO 41 J=1,N
      TEMP=HL(I+1,J)
      HL(I+1,J) = HL(I,J)
41    HL(I,J) = TEMP
42    IF(REAL(HL(I,I)).EQ.0..AND.AIMAG(HL(I,I)).EQ.0.)GO TO 44
      MULT(I) = -HL(I+1,I)/HL(I,I)
      DO 43 J = IP1,N
43    HL(I+1,J) = HL(I+1,J) + MULT(I) * HL(I,J)
44    CONTINUE
      DO 45 I = 1,N
45    VECT(I) = 1.
      TWICE = .FALSE.
46    IF(REAL(HL(N,N)).EQ.0..AND.AIMAG(HL(M,N)).EQ.0.)HL(N,N)=EPS
      VECT(N) = VECT(N)/HL(N,N)
      DO 48 I = 1,NM1
      K = N - I
      DO 47 J = K,NM1
47    VECT(K) = VECT(K)-HL(K,J+1)*VECT(J+1)
      IF(REAL(HL(K,K)).EQ.0..AND.AIMAG(HL(K,K)).EQ.0.)HL(K,K)=EPS
48    VECT(K) = VECT(K)/HL(K,K)
      BIG = 0.
      DO 49 I= 1,N
      SUM = ABS(REAL(VECT(I)))+ABS(AIMAG(VECT(I)))
49    IF(SUM.GT.BIG)BIG = SUM
      DO 50 I = 1,N
50    VECT(I) = VECT(I)/BIG
      IF(TWICE) GO TO 52
      DO 51 I = 1,NM1
      IF(.NOT.INTH(I)) GO TO 51
      TEMP = VECT(I)
      VECT(I) = VECT(I+1)
      VECT(I+1) = TEMP
51    VECT(I+1) = VECT(I+1)+MULT(I)*VECT(I)
      TWICE = .TRUE.
      GO TO 46
52    IF (N.EQ.2) GO TO 55
      MN2 = N - 2
      DO 54 I = 1,NM2
      N1I = N - 1 - I
      NI1=N-I+1
      DO 53 J = NI1,N
53    VECT(J) = H(J,N1I) * VECT(N1I+1) + VECT(J)
      INDEX = INT(N1I)
      TEMP = VECT(N1I+1)
      VECT(N1I+1) = VECT(INDEX)
54    VECT(INDEX) = TEMP
55    DO 56 I = 1,N
56    A(I,L) = VECT(I)
57    RETURN
      END
C---------------ALL PARAMETERS ARE INPUT
      SUBROUTINE CHECK(B,A,NCAL,N,EIGVAL,IRP)
      COMPLEX B(35,35),A(35,35),EIGVAL(35),AX(35),LX(35),DELTA(35),SUM
      DIMENSION RERR(35)
      REAL MAXRE
      DUMMY='1'
      DO 15 L=1,NCAL
      DO 1 I=1,N
      SUM=0.0
      DO 2 J=1,N
2     SUM=SUM+B(I,J)*A(J,L)
1     AX(I)=SUM
      DO 3 J=1,N
3     LX(J)=EIGVAL(L)*A(J,L)
      DO 4 J=1,N
4     DELTA(J)=AX(J)-LX(J)
      DO 8 J=1,N
8     RERR(J)=CABS(DELTA(J))/CABS(AX(J))
      K=1
      MAXRE=0.0
      DO 5 I=1,N
      SUM1=RERR(I)
      IF(SUM1-MAXRE)5,5,6
6     K=I
      MAXRE=SUM1
5     CONTINUE
      SUM1=0.0
      SUM2=0.0
      DO 9 I=1,N
      SUM1=SUM1+REAL(DELTA(I))**2+AIMAG(DELTA(I))**2
9     SUM2=SUM2+REAL(   AX(I))**2+AIMAG(   AX(I))**2
      IF(SUM2.LE.(.0000001)) GO TO 50
      GO TO 52
50    WRITE(IRP, 51)
51    FORMAT(' ','LENGTH OF AX IS SMALL;WE WILL NOT DIVIDE.'/)
      GO TO 15
52    TOTRE=SQRT(SUM1/SUM2)
      IF(L.NE.1)DUMMY='0'
      WRITE(IRP,7)DUMMY,K,AX(K),LX(K),DELTA(K),MAXRE,TOTRE
15    CONTINUE
7     FORMAT(A1 ,7X,'ROW WITH THE MAX. RELATIVE ERROR = 'I4/1X,4X,
     1'ROW WITH THE MAX. REL. ERROR HAS AX = ',G15.8,1X,G15.8/1X,4X,
     2'ROW WITH THE MAX. REL. ERROR HAS LX= ',G15.8,1X,G15.8/1X,5X,
     2'ROW WITH THE MAX. REL. ERROR HAS AX-LX= ',G15.8,1X,G15.8/1X,13X,
     3'THE MAXIMUM RELATIVE ERROR = ',G15.8/1X,
     4'EUCL. NORM OF (AX-LX)/EUCL. NORM OF AX = ',G15.8/)
      RETURN
      END

