Question

用fortran语言编写： 用雅克比方法求矩阵特征值和特征向量

用雅克比方法求矩阵特征值和特征向量

Answer 1

SUBROUTINE JCB(N,A,S,EPS)
C THIS SUBROUTINE COMPUTS EIGENVALUES AND standard EIGENVECTORS OF A
DIMENSION A(N,N),S(N,N)
DO 30 I=1,N
DO 30 J=1,I
IF(I-J) 20,10,20
10 S(I,J)=1.
GO TO 30
20 S(I,J)=0.
S(J,I)=0.
30 CONTINUE
G=0.
DO 40 I=2,N
I1=I-1
DO 40 J=1,I1
40 G=G+2.*A(I,J)*A(I,J)
S1=SQRT(G)
S2=EPS/FLOAT(N)*S1
S3=S1
L=0
50 S3=S3/FLOAT(N)
60 DO 130 IQ=2,N
IQ1=IQ-1
DO 130 IP=1,IQ1
IF(ABS(A(IP,IQ)).LT.S3) GOTO 130
L=1
V1=A(IP,IP)
V2=A(IP,IQ)
V3=A(IQ,IQ)
U=0.5*(V1-V3)
IF(U.EQ.0.0) G=1.
IF(ABS(U).GE.1E-10) G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U)
ST=G/SQRT(2.*(1.+SQRT(1.-G*G)))
CT=SQRT(1.-ST*ST)
DO 110 I=1,N
G=A(I,IP)*CT-A(I,IQ)*ST
A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
A(I,IP)=G
G=S(I,IP)*CT-S(I,IQ)*ST
S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT
110 S(I,IP)=G
DO 120 I=1,N
A(IP,I)=A(I,IP)
120 A(IQ,I)=A(I,IQ)
G=2.*V2*ST*CT
A(IP,IP)=V1*CT*CT+V3*ST*ST-G
A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
A(IQ,IP)=A(IP,IQ)
130 CONTINUE
IF(L-1) 150,140,150
140 L=0
GO TO 60
150 IF(S3.GT.S2) GOTO 50
RETURN
END