Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/evev/evev.for
There are 2 other files named evev.for in the archive. Click here to see a list.
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