用Jacobi法求实对称矩阵的特征值与特征向量(Fortran源程序)

zzs608

本人编写的用Jacobi法求实对称矩阵的特征值与特征向量，采用Fortran程序，经测试正确无误

1.         PROGRAM TEZHENG_Jacobi
   
2. * 用Jacobi方法求正交矩阵的特征值与特征向量
   
3. *   就是用平面旋转矩阵U不断对矩阵A作正交相似变换把A化为对角矩阵,
   
4. *   从而求出A的特征值与特征向量
   
5. *
   
6. * Made by ZhaoZunsheng,2005.12.20
   
7. *
   
8.         IMPLICIT NONE
   
9.         INTEGER NMAX,N,I,J,P,Q
   
10.         PARAMETER(NMAX=50)
    
11.         REAL A(NMAX,NMAX),AMAX,TEMP,ZEMP,COO,SII,CO,SI,APP,AQQ,APQ,API,AQI
    
12.         REAL R(NMAX,NMAX),RIP,RIQ
    
13.         CHARACTER NAME*12,NAMEO*12,CHR*1
    
14. * 从文件中读入实对称矩阵A
    
15.         WRITE(*,*)'输入实对称矩阵维数n(n<51):'
    
16.         READ(*,*) N
    
17.         WRITE(*,*)'输入矩阵文件:'
    
18.         READ(*,*) NAME
    
19.         OPEN(6,FILE=NAME)
    
20.         DO I=1,N
    
21.          READ(6,*) (A(I,J),J=1,I)
    
22.          DO J=1,I
    
23.          A(J,I)=A(I,J)
    
24.          ENDDO
    
25.         ENDDO
    
26.         CLOSE(6)
    
27. * R矩阵存放正交变换矩阵U,在这先初始化,即单位矩阵
    
28.         DO I=1,N
    
29.         DO J=1,N
    
30.         R(I,J)=0
    
31.         ENDDO
    
32.         R(I,I)=1
    
33.         ENDDO
    
34. * 在矩阵A的非主对角线元素中,找出按模最大的元素Apq
    
35. 100        AMAX=ABS(A(2,1))
    
36.         P=2
    
37.         Q=1
    
38.         DO I=2,N
    
39.          DO J=1,I-1
    
40.          IF(ABS(A(I,J)).GT.AMAX) THEN
    
41.          AMAX=ABS(A(I,J))
    
42.          P=I
    
43.          Q=J
    
44.          ENDIF
    
45.          ENDDO
    
46.         ENDDO
    
47. 
    
48. *        do i=1,n
    
49. *                WRITE(*,*)(A(I,J),J=1,N)
    
50. *        enddo
    
51. 
    
52. * 当非主对角线元素化为0,即小于给定精度时,输出特征值与特征向量
    
53.         IF(AMAX.LE.1.0E-7) THEN
    
54.         WRITE(*,*) 'A的特征值为:'
    
55.         WRITE(*,*) (A(I,I),I=1,N)
    
56.         WRITE(*,*) 'A的特征向量为:'
    
57.         WRITE(*,*) '  X1      X2     X3     ...:'
    
58.         DO I=1,N
    
59.         WRITE(*,*)(R(I,J),J=1,N)
    
60.         ENDDO
    
61. 
    
62.         WRITE(*,*) '是否将结果存入文件(Y/N)?'
    
63.         READ(*,*) CHR
    
64.         IF(CHR.EQ.'Y'.OR.CHR.EQ.'y') THEN
    
65.         WRITE(*,*) '输入文件名(小于12字符):'
    
66.         READ(*,*) NAMEO
    
67.         OPEN(8,FILE=NAMEO)
    
68.         WRITE(8,*) 'A的特征值为:'
    
69.         WRITE(8,*) (A(I,I),I=1,N)
    
70.         WRITE(8,*) 'A的特征向量为:'
    
71.         WRITE(8,*) '  X1      X2     X3     ...:'
    
72.         DO I=1,N
    
73.         WRITE(8,*)(R(I,J),J=1,N)
    
74.         ENDDO
    
75.         CLOSE(8)
    
76.         ENDIF
    
77. 
    
78.         STOP
    
79.         ENDIF
    
80. * 开始准备计算平面旋转矩阵U
    
81.         TEMP=2*A(P,Q)/(A(P,P)-A(Q,Q)+1.0e-30)
    
82.         ZEMP=(A(P,P)-A(Q,Q))/(2*A(P,Q))
    
83.        
    
84.         IF(ABS(TEMP).LT.1.0) THEN
    
85.         COO=(1+TEMP**2)**(-0.5)
    
86.         SII=TEMP*(1+TEMP**2)**(-0.5)
    
87.         ELSE
    
88.         COO=ABS(ZEMP)*(1+ZEMP**2)**(-0.5)
    
89.         SII=SIGN(1.0,ZEMP)*(1+ZEMP**2)**(-0.5)
    
90.         ENDIF
    
91.        
    
92.         CO=SQRT(0.5*(1+COO))
    
93.         SI=SII/(2*CO)
    
94. * 计算平面旋转矩阵U
    
95.         DO I=1,N
    
96.         RIP=R(I,P)*CO+R(I,Q)*SI
    
97.         RIQ=-R(I,P)*SI+R(I,Q)*CO
    
98.         R(I,P)=RIP
    
99.         R(I,Q)=RIQ
    
100.         ENDDO
     
101. * 对A进行变换       
     
102.         APP=A(P,P)*CO**2+A(Q,Q)*SI**2+2*A(P,Q)*CO*SI
     
103.         AQQ=A(P,P)*SI**2+A(Q,Q)*CO**2-2*A(P,Q)*CO*SI
     
104.         APQ=0.5*(A(Q,Q)-A(P,P))*SII+A(P,Q)*COO
     
105.         A(P,P)=APP
     
106.         A(Q,Q)=AQQ
     
107.         A(P,Q)=APQ
     
108.         A(Q,P)=A(P,Q)
     
109.        
     
110.         DO I=1,N
     
111.         IF(I.EQ.P.OR.I.EQ.Q) THEN
     
112.         ELSE
     
113.         API=A(P,I)*CO+A(Q,I)*SI
     
114.         AQI=-A(P,I)*SI+A(Q,I)*CO
     
115.         A(P,I)=API
     
116.         A(Q,I)=AQI
     
117.         A(I,P)=A(P,I)
     
118.         A(I,Q)=A(Q,I)
     
119.         ENDIF
     
120.         ENDDO
     
121.        
     
122.       GOTO 100  
     
123.         END

复制代码

风花雪月

学习学习

zhengdatou

<P>好好学学！！</P>

linqus

<P>谢谢楼主分享[em01]</P>

xinshao

谢谢分享

jiwell

学习一下。谢谢了

runchi0610

<P>非常感谢楼主啦！！我正好有急用啊~~！！！</P>

johhan

谢谢,好好学习一下!!!!!!!

jcp

谢谢楼主了

xuang

<P>谢谢分享,学习中!!</P>

wenzhuhust

nice
3q for sharing

zhangbo1030

perfect! thanks

魏红

非常感谢~~~:handshake