- #1
draconidz
- 1
- 0
Homework Statement
I have two subroutine syntax program in fotran77 code, and I am stuck to understand it.
this is the first subroutine code
Fortran:
subroutine LZHES(N,A,NA,B,NB,X,NX,WANTX)
c
implicit real*8 (a-h,o-z)
c
complex*16 Y,W,Z,A(NA,N),B(NB,N),X(NX,N)
real*8 C,D
logical WANTX
c
NM1=N-1
DO 30 I=1,NM1
D=0.0
IP1=I+1
DO 10 K=IP1,N
Y=B(K,I)
C=dabs(dreal(Y))+dabs(dimag(Y))
IF (C.LE.D) GOTO 10
D=C
II=K
10 CONTINUE
IF (D.EQ.0.0) GOTO 30
Y=B(I,I)
IF(D.LE.dabs(dreal(Y))+dabs(dimag(Y))) GOTO 15
DO 11 J=1,N
Y=A(I,J)
A(I,J)=A(II,J)
A(II,J)=Y
11 CONTINUE
DO 12 J=I,N
Y=B(I,J)
B(I,J)=B(II,J)
B(II,J)=Y
12 CONTINUE
15 DO 20 J=IP1,N
Y=B(J,I)/B(I,I)
IF (dreal(Y).EQ.0.0.AND.dimag(Y).EQ.0.0) GOTO 20
DO 18 K=1,N
A(J,K)=A(J,K)-Y*A(I,K)
18 CONTINUE
DO 19 K=IP1,N
B(J,K)=B(J,K)-Y*B(I,K)
19 CONTINUE
20 CONTINUE
B(IP1,I)=(0.0,0.0)
30 CONTINUE
IF (.NOT.WANTX) GOTO 40
DO 38 I=1,N
DO 37 J=1,N
X(I,J)=(0.0,0.0)
37 CONTINUE
X(I,I)=(1.0,0.0)
38 CONTINUE
40 NM2=N-2
IF (NM2.LT.1) GOTO 100
DO 90 J=1,NM2
JM2=NM1-J
JP1=J+1
DO 80 II=1,JM2
I=N+1-II
IM1=I-1
IMJ=I-J
W=A(I,J)
Z=A(IM1,J)
IF (dabs(dreal(W))+dabs(dimag(W)).LE.dabs(dreal(Z))+
& dabs(dimag(Z))) GOTO 50
DO 45 K=J,N
Y=A(I,K)
A(I,K)=A(IM1,K)
A(IM1,K)=Y
45 CONTINUE
DO 46 K=IM1,N
Y=B(I,K)
B(I,K)=B(IM1,K)
B(IM1,K)=Y
46 CONTINUE
50 Z=A(I,J)
IF (dreal(Z).EQ.0.0.AND.dimag(Z).EQ.0.0) GOTO 58
Y=Z/A(IM1,J)
DO 52 K=JP1,N
A(I,K)=A(I,K)-Y*A(IM1,K)
52 CONTINUE
DO 54 K=IM1,N
B(I,K)=B(I,K)-Y*B(IM1,K)
54 CONTINUE
58 W=B(I,IM1)
Z=B(I,I)
IF(dabs(dreal(W))+dabs(dimag(W)).LE.dabs(dreal(Z))+
& dabs(dimag(Z))) GOTO 70
DO 60 K=1,I
Y=B(K,I)
B(K,I)=B(K,IM1)
B(K,IM1)=Y
60 CONTINUE
DO 64 K=1,N
Y=A(K,I)
A(K,I)=A(K,IM1)
A(K,IM1)=Y
64 CONTINUE
IF (.NOT.WANTX) GOTO 70
DO 68 K=IMJ,N
Y=X(K,I)
X(K,I)=X(K,IM1)
X(K,IM1)=Y
68 CONTINUE
70 Z=B(I,IM1)
IF (dreal(Z).EQ.0.0.AND.dimag(Z).EQ.0.0) GOTO 80
Y=Z/B(I,I)
DO 72 K=1,IM1
B(K,IM1)=B(K,IM1)-Y*B(K,I)
72 CONTINUE
B(I,IM1)=(0.0,0.0)
DO 74 K=1,N
A(K,IM1)=A(K,IM1)-Y*A(K,I)
74 CONTINUE
IF (.NOT.WANTX) GOTO 80
DO 76 K=IMJ,N
X(K,IM1)=X(K,IM1)-Y*X(K,I)
76 CONTINUE
80 CONTINUE
A(JP1+1,J)=(0.0,0.0)
90 CONTINUE
100 RETURN
END
Fortran:
subroutine LZIT(N,A,NA,B,NB,X,NX,WANTX,ITER,EIGA,EIGB)
c
implicit real*8 (a-h,o-z)
c
integer ITER(N)
real*8 EPSA,EPSB,SS,R,ANORM,BNORM,ANI,BNI,C
real*8 D0,D1,D2,E0,E1
complex*16 A(NA,N),B(NB,N),EIGA(N),EIGB(N)
complex*16 X(NX,N)
complex*16 S,W,Y,Z
complex*16 ANNM1,ALFM,BETM,D,SL,DEN,NUM,ANM1M1
logical WANTX
NN=N
ANORM=0.
BNORM=0.
DO 5 I=1,N
ANI=0.
IF (I.EQ.1) GOTO 2
Y=A(I,I-1)
ANI=ANI+dabs(dreal(Y))+dabs(dimag(Y))
2 BNI=0.
DO 3 J=I,N
ANI=ANI+dabs(dreal(A(I,J)))+dabs(dimag(A(I,J)))
BNI=BNI+dabs(dreal(B(I,J)))+dabs(dimag(B(I,J)))
3 CONTINUE
IF (ANI.GT.ANORM) ANORM=ANI
IF (BNI.GT.BNORM) BNORM=BNI
5 CONTINUE
IF (ANORM.EQ.0.) ANORM=1.0
IF (BNORM.EQ.0.) BNORM=1.0
EPSB=BNORM
EPSA=ANORM
6 EPSA=EPSA/2.0
EPSB=EPSB/2.0
C=ANORM+EPSA
IF (C.GT.ANORM) GOTO 6
IF (N.LE.1) GOTO 100
10 ITS=0
NM1=NN-1
11 D2=dabs(dreal(A(NN,NN)))+dabs(dimag(A(NN,NN)))
DO 12 LB=2,NN
L=NN+2-LB
SS=D2
Y=A(L-1,L-1)
D2=dabs(dreal(Y))+dabs(dimag(Y))
SS=SS+D2
Y=A(L,L-1)
R=SS+dabs(dreal(Y))+dabs(dimag(Y))
IF (R.EQ.SS) GOTO 13
12 CONTINUE
L=1
13 IF (L.EQ.NN) GOTO 100
IF (ITS.LT.30) GOTO 20
ITER(NN)=-1
IF (dabs(dreal(A(NN,NM1)))+dabs(dimag(A(NN,NM1))).GT.0.8*
1 dabs(dreal(ANNM1))+dabs(dimag(ANNM1))) RETURN
20 IF(ITS.EQ.10.OR.ITS.EQ.20) GOTO 28
ANNM1=A(NN,NM1)
ANM1M1=A(NM1,NM1)
S=A(NN,NN)*B(NM1,NM1)-ANNM1*B(NM1,NN)
W=ANNM1*B(NN,NN)*(A(NM1,NN)*B(NM1,NM1)-B(NM1,NN)*ANM1M1)
Y=(ANM1M1*B(NN,NN)-S)/2.
Z=CDSQRT(Y*Y+W)
IF (dreal(Z).EQ.0.0.AND.dimag(Z).EQ.0.0) GOTO 26
D0=Y/Z
IF (D0.LT.0.0) Z=-Z
26 DEN=(Y+Z)*B(NM1,NM1)*B(NN,NN)
IF (dreal(DEN).EQ.0.0.AND.dimag(DEN).EQ.0.0) DEN=EPSA
NUM=(Y+Z)*S-W
GOTO 30
28 Y=A(NM1,NN-2)
NUM=dCMPLX(dabs(dreal(ANNM1))+dabs(dimag(ANNM1)),dabs(dreal(Y))+
1dabs(dimag(Y)))
DEN=(1.0,0.0)
30 IF (NN.EQ.L+1) GOTO 35
D2=dabs(dreal(A(NM1,NM1)))+dabs(dimag(A(NM1,NM1)))
E1=dabs(dreal(ANNM1))+dabs(dimag(ANNM1))
D1=dabs(dreal(A(NN,NN)))+dabs(dimag(A(NN,NN)))
NL=NN-(L+1)
DO 34 MB=1,NL
M=NN-MB
E0=E1
Y=A(M,M-1)
E1=dabs(dreal(Y))+dabs(dimag(Y))
D0=D1
D1=D2
Y=A(M-1,M-1)
D2=dabs(dreal(Y))+dabs(dimag(Y))
Y=A(M,M)*DEN-B(M,M)*NUM
D0=(D0+D1+D2)*(dabs(dreal(Y))+dabs(dimag(Y)))
E0=E0*E1*(dabs(dreal(DEN))+dabs(dimag(DEN)))+D0
IF(E0.EQ.D0) GOTO 36
34 CONTINUE
35 M=L
36 CONTINUE
ITS=ITS+1
W=A(M,M)*DEN-B(M,M)*NUM
Z=A(M+1,M)*DEN
D1=dabs(dreal(Z))+dabs(dimag(Z))
D2=dabs(dreal(W))+dabs(dimag(W))
NP1=N+1
LOR1=L
NNORN=NN
IF (.NOT.WANTX) GOTO 42
LOR1=1
NNORN=N
42 DO 90 I=M,NM1
J=I+1
IF (I.EQ.M) GOTO 50
W=A(I,I-1)
Z=A(J,I-1)
D1=dabs(dreal(Z))+dabs(dimag(Z))
D2=dabs(dreal(W))+dabs(dimag(W))
IF (D1.EQ.0.0) GOTO 11
50 IF (D2.GT.D1) GOTO 60
DO 55 K=I,NNORN
Y=A(I,K)
A(I,K)=A(J,K)
A(J,K)=Y
Y=B(I,K)
B(I,K)=B(J,K)
B(J,K)=Y
55 CONTINUE
IF (I.GT.M) A(I,I-1)=A(J,I-1)
IF(D2.EQ.0.0) GOTO 65
Y=dCMPLX(dreal(W)/D1,dimag(W)/D1)/dCMPLX(dreal(Z)/D1,
& dimag(Z)/D1)
GOTO 62
60 Y=dCMPLX(dreal(Z)/D2,dimag(Z)/D2)/dCMPLX(dreal(W)/D2,dimag(W)/D2)
62 DO 64 K=I,NNORN
A(J,K)=A(J,K)-Y*A(I,K)
B(J,K)=B(J,K)-Y*B(I,K)
64 CONTINUE
IF (I.GT.M) A(J,I-1)=(0.0,0.0)
65 Z=B(J,I)
W=B(J,J)
D2=dabs(dreal(W))+dabs(dimag(W))
D1=dabs(dreal(Z))+dabs(dimag(Z))
IF (D1.EQ.0.0) GOTO 11
IF (D2.GT.D1) GOTO 80
DO 70 K=LOR1,J
Y=A(K,J)
A(K,J)=A(K,I)
A(K,I)=Y
Y=B(K,J)
B(K,J)=B(K,I)
B(K,I)=Y
70 CONTINUE
IF(I.EQ.NM1) GOTO 75
Y=A(J+1,J)
A(J+1,J)=A(J+1,I)
A(J+1,I)=Y
75 IF(.NOT.WANTX) GOTO 79
DO 78 K=1,N
Y=X(K,J)
X(K,J)=X(K,I)
X(K,I)=Y
78 CONTINUE
79 IF(D2.EQ.0.0) GOTO 90
Z=dCMPLX(dreal(W)/D1,dimag(W)/D1)/dCMPLX(dreal(Z)/D1,dimag(Z)/D1)
GOTO 81
80 Z=dCMPLX(dreal(Z)/D2,dimag(Z)/D2)/dCMPLX(dreal(W)/D2,dimag(W)/D2)
81 DO 82 K=LOR1,J
A(K,I)=A(K,I)-Z*A(K,J)
B(K,I)=B(K,I)-Z*B(K,J)
82 CONTINUE
B(J,I)=(0.0,0.0)
IF (I.LT.NM1) A(I+2,I)=A(I+2,I)-Z*A(I+2,J)
IF(.NOT.WANTX) GOTO 90
DO 85 K=1,N
X(K,I)=X(K,I)-Z*X(K,J)
85 CONTINUE
90 CONTINUE
GOTO 11
100 CONTINUE
EIGA(NN)=A(NN,NN)
EIGB(NN)=B(NN,NN)
IF (NN.EQ.1) GOTO 110
ITER(NN)=ITS
NN=NM1
IF (NN.GT.1) GOTO 10
ITER(1)=0
GOTO 100
110 IF(.NOT.WANTX) RETURN
M=N
115 CONTINUE
ALFM=A(M,M)
BETM=B(M,M)
B(M,M)=(1.0,0.0)
L=M-1
IF (L.EQ.0) GOTO 140
120 CONTINUE
L1=L+1
SL=0.
DO 130 J=L1,M
SL=SL+(BETM*A(L,J)-ALFM*B(L,J))*B(J,M)
130 CONTINUE
Y=BETM*A(L,L)-ALFM*B(L,L)
IF(dreal(Y).EQ.0.0.AND.dimag(Y).EQ.0.0) Y=(EPSA+EPSB)/2.0
B(L,M)=-SL/Y
L=L-1
140 IF (L.GT.0) GOTO 120
M=M-1
IF (M.GT.0) GOTO 115
M=N
200 CONTINUE
DO 220 I=1,N
S=0.
DO 210 J=1,M
S=S+X(I,J)*B(J,M)
210 CONTINUE
X(I,M)=S
220 CONTINUE
M=M-1
IF (M.GT.0) GOTO 200
M=N
230 CONTINUE
SS=0.
DO 235 I=1,N
R=dabs(dreal(X(I,M)))+dabs(dimag(X(I,M)))
IF (R.LT.SS) GOTO 235
SS=R
D=X(I,M)
235 CONTINUE
IF(SS.EQ.0.0) GOTO 245
DO 240 I=1,N
X(I,M)=X(I,M)/D
240 CONTINUE
245 M=M-1
IF (M.GT.0) GOTO 230
RETURN
ENDand the question is
I don't know exactly the subroutine doing what,
and what the subroutine output and input
Homework Equations
3. The Attempt at a Solution [/B]
I thought the first sub routine (LZHES) trying to diagonalization matrix B
and I thought the first subroutine input is
N,A,NA,B,NB,X,NX,WANTX
N is the size of matrix
A is the identity matrix
B is the matrix will be diagonalize
NA, NB is the size of matrix A and matrix B
WANTX is the variable to decided to diagonalize matrix or not.
and the subroutine output, I thought is B, the diagonalize matrix B.
but this subroutine when I tried to created in MATLAB code version, matrix B is not diagonalized.
and I thought this subroutine LZHES is connected to the subroutine LZIT.
and the second subroutine trying to find the eigen values of the matrix B.
I thought the first subroutine input is
N,A,NA,B,NB,X,NX,WANTX,ITER,EIGA,EIGB
N is the size of matrix
A is the identity matrix
B is the matrix will be diagonalize
NA, NB is the size of matrix A and matrix B
WANTX is the variable to decided to diagonalize matrix or not.
ITER is the number of the iteration
EIGA is the initialize eigen value
EIGB is the initilaize eigen value
and the output, i thought is the eigen values which is EIGA and EIGB
Last edited by a moderator:
