# Understanding fortran77 syntax

1. Jun 27, 2015

### draconidz

1. The problem statement, all variables and given/known data
I have two subroutine syntax program in fotran77 code, and I am stuck to understand it.
this is the first subroutine code
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

and this is the second subroutine code
Code (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
END

and the question is
I don't know exactly the subroutine doing what,
and what the subroutine output and input

2. Relevant equations

3. The attempt at a solution

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: Jun 27, 2015
2. Jun 27, 2015

### SteamKing

Staff Emeritus
An article was published in 1975 in the ACM journal Transactions on Mathematical Software (TOMS), as cited here:

http://dl.acm.org/citation.cfm?doid=355644.355652

The original Fortran source code (with comments) for these two subroutines LZHES and LZIT can be recovered from the ACM Collected Algorithms page located here:

http://netlib.org/toms/