zzs608 发表于 2005-12-19 20:13

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

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

        PROGRAM TEZHENG_Jacobi
* 用Jacobi方法求正交矩阵的特征值与特征向量
*   就是用平面旋转矩阵U不断对矩阵A作正交相似变换把A化为对角矩阵,
*   从而求出A的特征值与特征向量
*
* Made by ZhaoZunsheng,2005.12.20
*
        IMPLICIT NONE
        INTEGER NMAX,N,I,J,P,Q
        PARAMETER(NMAX=50)
        REAL A(NMAX,NMAX),AMAX,TEMP,ZEMP,COO,SII,CO,SI,APP,AQQ,APQ,API,AQI
        REAL R(NMAX,NMAX),RIP,RIQ
        CHARACTER NAME*12,NAMEO*12,CHR*1
* 从文件中读入实对称矩阵A
        WRITE(*,*)'输入实对称矩阵维数n(n<51):'
        READ(*,*) N
        WRITE(*,*)'输入矩阵文件:'
        READ(*,*) NAME
        OPEN(6,FILE=NAME)
        DO I=1,N
       READ(6,*) (A(I,J),J=1,I)
       DO J=1,I
       A(J,I)=A(I,J)
       ENDDO
        ENDDO
        CLOSE(6)
* R矩阵存放正交变换矩阵U,在这先初始化,即单位矩阵
        DO I=1,N
        DO J=1,N
        R(I,J)=0
        ENDDO
        R(I,I)=1
        ENDDO
* 在矩阵A的非主对角线元素中,找出按模最大的元素Apq
100        AMAX=ABS(A(2,1))
        P=2
        Q=1
        DO I=2,N
       DO J=1,I-1
       IF(ABS(A(I,J)).GT.AMAX) THEN
       AMAX=ABS(A(I,J))
       P=I
       Q=J
       ENDIF
       ENDDO
        ENDDO

*        do i=1,n
*                WRITE(*,*)(A(I,J),J=1,N)
*        enddo

* 当非主对角线元素化为0,即小于给定精度时,输出特征值与特征向量
        IF(AMAX.LE.1.0E-7) THEN
        WRITE(*,*) 'A的特征值为:'
        WRITE(*,*) (A(I,I),I=1,N)
        WRITE(*,*) 'A的特征向量为:'
        WRITE(*,*) 'X1      X2   X3   ...:'
        DO I=1,N
        WRITE(*,*)(R(I,J),J=1,N)
        ENDDO

        WRITE(*,*) '是否将结果存入文件(Y/N)?'
        READ(*,*) CHR
        IF(CHR.EQ.'Y'.OR.CHR.EQ.'y') THEN
        WRITE(*,*) '输入文件名(小于12字符):'
        READ(*,*) NAMEO
        OPEN(8,FILE=NAMEO)
        WRITE(8,*) 'A的特征值为:'
        WRITE(8,*) (A(I,I),I=1,N)
        WRITE(8,*) 'A的特征向量为:'
        WRITE(8,*) 'X1      X2   X3   ...:'
        DO I=1,N
        WRITE(8,*)(R(I,J),J=1,N)
        ENDDO
        CLOSE(8)
        ENDIF

        STOP
        ENDIF
* 开始准备计算平面旋转矩阵U
        TEMP=2*A(P,Q)/(A(P,P)-A(Q,Q)+1.0e-30)
        ZEMP=(A(P,P)-A(Q,Q))/(2*A(P,Q))
       
        IF(ABS(TEMP).LT.1.0) THEN
        COO=(1+TEMP**2)**(-0.5)
        SII=TEMP*(1+TEMP**2)**(-0.5)
        ELSE
        COO=ABS(ZEMP)*(1+ZEMP**2)**(-0.5)
        SII=SIGN(1.0,ZEMP)*(1+ZEMP**2)**(-0.5)
        ENDIF
       
        CO=SQRT(0.5*(1+COO))
        SI=SII/(2*CO)
* 计算平面旋转矩阵U
        DO I=1,N
        RIP=R(I,P)*CO+R(I,Q)*SI
        RIQ=-R(I,P)*SI+R(I,Q)*CO
        R(I,P)=RIP
        R(I,Q)=RIQ
        ENDDO
* 对A进行变换       
        APP=A(P,P)*CO**2+A(Q,Q)*SI**2+2*A(P,Q)*CO*SI
        AQQ=A(P,P)*SI**2+A(Q,Q)*CO**2-2*A(P,Q)*CO*SI
        APQ=0.5*(A(Q,Q)-A(P,P))*SII+A(P,Q)*COO
        A(P,P)=APP
        A(Q,Q)=AQQ
        A(P,Q)=APQ
        A(Q,P)=A(P,Q)
       
        DO I=1,N
        IF(I.EQ.P.OR.I.EQ.Q) THEN
        ELSE
        API=A(P,I)*CO+A(Q,I)*SI
        AQI=-A(P,I)*SI+A(Q,I)*CO
        A(P,I)=API
        A(Q,I)=AQI
        A(I,P)=A(P,I)
        A(I,Q)=A(Q,I)
        ENDIF
        ENDDO
       
      GOTO 100
        END

风花雪月 发表于 2005-12-22 08:14

回复：（zzs608）用Jacobi法求实对称矩阵的特征值与特...

学习学习

zhengdatou 发表于 2005-12-23 15:42

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

linqus 发表于 2005-12-24 19:25

<P>谢谢楼主分享</P>

xinshao 发表于 2005-12-25 11:25

回复：（zzs608）用Jacobi法求实对称矩阵的特征值与特...

谢谢分享

jiwell 发表于 2005-12-30 14:37

学习一下。谢谢了

runchi0610 发表于 2006-1-2 14:48

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

johhan 发表于 2006-1-3 21:29

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

jcp 发表于 2006-2-20 09:14

谢谢楼主了

xuang 发表于 2006-3-16 09:40

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

wenzhuhust 发表于 2010-4-6 21:15

nice
3q for sharing

zhangbo1030 发表于 2010-7-6 08:33

perfect! thanks

魏红 发表于 2010-7-13 12:04

回复 楼主 zzs608 的帖子

非常感谢~~~:handshake
页: [1]
查看完整版本: 用Jacobi法求实对称矩阵的特征值与特征向量(Fortran源程序)