Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50144/eig1.f4
There are no other files named eig1.f4 in the archive.
C THE PROGRAM EIG1D.F4 IS A SAMPLE PROGRAM, WHEN USED IN
C CONJUNCTION WITH EIG1(LOADED AND SAVED AS EIG1D.SAV)
C WILL PRODUCE A RUNNING EXAMPLE OF THE USE OF EIG1.
SUBROUTINE EIG1(A,B,N,EPS,AM,IR,NA,NB)
DIMENSION A(500),B(500),AM(500),IR(500)
IJSF(I,J)=N*(I-1)+J+(I*(1-I))/2
IJDF(I,J)=NA*(J-1)+I
1 EN=N
3 IF(NA-1)7,5,7
5 MF=1
GO TO 10
7 MF=0
10 IF(EPS)601,11,12
11 SE=1.E-6
GO TO 15
12 SE=SQRT(EPS)
15 DO 60 I=1,N
DO 50 J=1,N
IJ=NB*(J-1)+I
50 B(IJ)=0.
AM(I)=0.
II=NB*(I-1)+I
60 B(II)=1.
80 DO 90 J=2,N
K=J-1
DO 90 I=1,K
IF(MF)601,82,83
82 L=IJDF(I,J)
GO TO 84
83 L=IJSF(I,J)
84 T=ABS (A(L))
IF(T-AM(J))90,85,85
85 AM(J)=T
IR(J)=I
90 CONTINUE
100 BA=0.
DO 125 J=2,N
IF(AM(J)-BA)125,120,120
120 BA=AM(J)
L=IR(J)
M=J
125 CONTINUE
IF(MF)601,126,127
126 LM=IJDF(L,M)
MM=LM-L+M
LL=IJDF(L,L)
GO TO 138
127 LM=IJSF(L,M)
LL=LM-M+L
MM=IJSF(M,M)
138 IF(EN*ABS(A(LM))-SE)590,150,150
150 LM1=L-1
LP1=L+1
MM1=M-1
MP1=M+1
200 R=SQRT((A(LL)-A(MM))**2+4.*A(LM)*A(LM))
T=A(LL)-A(MM)
400 C=SQRT(.5*(1.+ABS(T)/R))
S=-A(LM)/(R*C)
IF(T)401,402,402
401 S=-S
402 SS=S*S
CC=C*C
SC=S*C
P=2.*A(LM)*SC
A(LM)=0
AM(L)=0
AM(M)=0
IF(LM1)411,450,411
411 DO 440 I=1,LM1
IF(MF)601,412,413
412 IL=IJDF(I,L)
IM=IJDF(I,M)
GO TO 414
413 IL=IJSF(I,L)
IM=IL-L+M
414 TE=A(IL)
A(IL)=A(IL)*C-S*A(IM)
A(IM)=A(IM)*C+TE*S
TEM=ABS(A(IL))
IF(TEM-AM(L))430,415,415
415 AM(L)=TEM
IR(L)=I
430 TEM=ABS(A(IM))
IF(TEM-AM(M))440,435,435
435 AM(M)=TEM
IR(M)=I
440 CONTINUE
450 IF(LP1-M)451,500,451
451 DO 490 K=LP1,MM1
IF(MF)601,462,463
462 LK=IJDF(L,K)
LT=LK-L
KM=IJDF(K,M)
GO TO 464
463 LK=IJSF(L,K)
KM=IJSF(K,M)
464 TE=A(LK)
A(LK)=A(LK)*C-A(KM)*S
A(KM)=A(KM)*C+TE*S
IF(IR(K)-L)465,470,465
465 TEM=ABS(A(LK))
IF(TEM-AM(K))480,467,467
467 AM(K)=TEM
IR(K)=L
GO TO 480
470 AM(K)=0.
KM1=K-1
DO 475 IT=1,KM1
IF(MF)601,471,472
471 ITK=LT+IT
GO TO 473
472 ITK=IJSF(IT,K)
473 TEM=ABS(A(ITK))
IF(TEM-AM(K))475,474,474
474 AM(K)=TEM
IR(K)=IT
475 CONTINUE
480 TEM=ABS(A(KM))
IF(TEM-AM(M))490,485,485
485 AM(M)=TEM
IR(M)=K
490 CONTINUE
500 IF(M-N)511,580,511
511 DO 540 J=MP1,N
IF(MF)601,512,513
512 LJ=IJDF(L,J)
LT=LJ-L
MJ=LT+M
GO TO 514
513 LJ=IJSF(L,J)
MJ=IJSF(M,J)
514 TE=A(LJ)
A(LJ)=A(LJ)*C-A(MJ)*S
A(MJ)=A(MJ)*C+TE*S
IF(IR(J)-L)515,525,515
515 IF(IR(J)-M)517,525,517
517 TEM=ABS(A(LJ))
IF(TEM-AM(J))520,518,518
518 AM(J)=TEM
IR(J)=L
520 TEM=ABS(A(MJ))
IF(TEM-AM(J))540,521,521
521 AM(J)=TEM
IR(J)=M
GO TO 540
525 AM(J)=0.
JM1=J-1
DO 530 IT=1,JM1
IF(MF)601,526,527
526 ITJ=LT+IT
GO TO 528
527 ITJ=IJSF(IT,J)
528 TEM=ABS(A(ITJ))
IF(TEM-AM(J))530,529,529
529 AM(J)=TEM
IR(J)=IT
530 CONTINUE
540 CONTINUE
580 TE=A(LL)
A(LL)=A(LL)*CC-P+A(MM)*SS
A(MM)=A(MM)*CC+P+TE*SS
581 DO 585 I=1,N
IL=NB*(L-1)+I
IM=NB*(M-1)+I
TE=B(IL)
B(IL)=B(IL)*C-B(IM)*S
B(IM)=B(IM)*C+TE*S
585 CONTINUE
GO TO 100
590 IF(MF) 601,601,602
602 DO 600 I=2,N
II=IJSF(I,I)
600 A(I)=A(II)
601 END