声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3800|回复: 7

[Fortran] Jacobi 广义对称矩阵特征值计算子程序

[复制链接]
发表于 2008-5-24 09:11 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
.
    这里给出的是满阵形式存储的Jacobi法广义对称矩阵特征值计算的子程序代码,来自《 NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS(K. J. Bathe  E. L. Wilson)》:

C......................................................................
C                                                                     .
C  PROGRAM TO SOLVE THE GRNERALIZED EIGENPROBLEM USING THE            .
C                                   GENERALIZED JACOBI ITERATION      .
C   INPUT VARIAELES :                                                 .
C           A(N,N) = STIFFNESS MATRIX  (ASSUMPD POSITIVE DEPINITE)    .
C           B(N,N) = MASS MATRIX       (ASSUMPD POSITIVE DEPINITE)    .
C           X(N,N) = MATRIX STORING EIGENVECTORS ON SOLUTION EXIT     .
C           EIGV(N)= VECTOR STORING EIGENVALUES  ON SOLUTION EXIT     .
C           D(N)   = WORKING VERTOR                                   .
C           M(N)   = WORKING VERTOR                                   .
C           N      = ORDER OF NATRICES A AND B                        .
C           RTOL   = CONVERGENCE TOLERANCE (USUALLY SET TO 10^-12)    .
C                                                                     .
C   OUTPUT VARIAELES :                                                .
C           X(N,N) = EIGENVECTORS STORED COLUMNNISE                   .
C           EIGV(N)= EIGENVALUES                                      .
C                                                                     .
C                                                                     .
C                  NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS       .
C                                                                     .
C                             by  K. J. Bathe  E. L. Wilson           .
C                                                                     .
C                                                                     .
C......................................................................
        SUBROUTINE JACOBI(A,B,X,EIGV,D,M,N,RTOL)
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION A(N,N),B(N,N),X(N,N),EIGV(N),D(N),M(N)
C......................................................................
C                                                                     .
C  THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CDC         .
C    EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM OR UNIVAC       .
C    NACHINES . ACTIVATE , DEACTIVATE OR ADJUST ABOVE CARDS FOR       .
C    SIMGLE OR DOUBLE PRECISION ARITHMETIC                            .
C                                                                     .
C......................................................................
C----------------------------------------------------------------------
C               INITIALIZE EIGEVALUE AND EIGENVECTOR BATRICES
C----------------------------------------------------------------------
            DO 10 I=1,N
           IF(A(I,I).GT.0.AND.B(I,I).GT.0.) GOTO 4
           WRITE(*,2020)
2020   FORMAT(10X,' *****   Error  solution   ***** ')
           STOP
    4     D(I)=A(I,I)/B(I,I)
   10    EIGV(I)=D(I)
           DO 30 I=1,N
           DO 20 J=1,N
   20    X(I,J)=0.
   30    X(I,I)=1.
           IF(N.EQ.1) RETURN
C----------------------------------------------------------------------
C               INITIALIZE SWEEP COUNTER AND ITERATION
C----------------------------------------------------------------------
            NSWEEP=0
            NR=N-1
   40     NSWEEP=NSWEEP+1
            WRITE(*,2000) NSWEEP
2000   FORMAT(10X,'Sweep number is',I4)
C--------------------------------------------------------------------------
C  CHECK IF PRESENT OFF-DIAGCNAL ELEMENT IS LARGE ENOUGH TO REQUIRE ZEROING
C--------------------------------------------------------------------------
        EPS=(.01**NSWEEP)**2
        DO 210 J=1,NR
        JJ=J+1
        DO 210 K=JJ,N
        EPTOLA=(A(J,K)*A(J,K))/(A(J,J)*A(K,K))
        EPTOLB=(B(J,K)*B(J,K))/(B(J,J)*B(K,K))
        IF(EPTOLA.LT.EPS.AND.EPTOLB.LT.EPS) GOTO 210
C--------------------------------------------------------------------------
C  IF ZEROING IS REQUIRED, CALCULATE THE ROTATION MATRIX ELEMENTS CA AND CG
C--------------------------------------------------------------------------
        AKK=A(K,K)*B(J,K)-B(K,K)*A(J,K)
        AJJ=A(J,J)*B(J,K)-B(J,J)*A(J,K)
        AB=A(J,J)*B(K,K)-A(K,K)*B(J,J)
        CHECK=(AB*AB+4.*AKK*AJJ)/4.
        IF(CHECK) 50,60,60
   50   WRITE(*,2020)
        STOP
   60   SQCH=DSQRT(CHECK)
        D1=AB/2.+SQCH
        D2=AB/2.-SQCH
        DEN=D1
        IF(DABS(D2).GT.DABS(D1)) DEN=D2
        IF(DEN) 80,70,80
   70   CA=0.
        CG=-A(J,K)/A(K,K)
        GOTO 90
   80   CA=AKK/DEN
        CG=-AJJ/DEN
C---------------------------------------------------------------------------
C  PERPORM THE GENERALIZED ROTATION TO ZERO THE PRESENT OFF-DIAGONAL ELEMENT
C---------------------------------------------------------------------------
   90   IF(N-2) 100,190,100
  100   JP1=J+1
        JM1=J-1
        KP1=K+1
        KM1=K-1
        IF(JM1-1) 130,110,110
  110   DO 120 I=1,JM1
        AJ=A(I,J)
        BJ=B(I,J)
        AK=A(I,K)
        BK=B(I,K)
        A(I,J)=AJ+CG*AK
        B(I,J)=BJ+CG*BK
        A(I,K)=AK+CA*AJ
  120   B(I,K)=BK+CA*BJ
  130   IF(KP1-N) 140,140,160
  140   DO 150 I=KP1,N
        AJ=A(J,I)
        BJ=B(J,I)
        AK=A(K,I)
        BK=B(K,I)
        A(J,I)=AJ+CG*AK
        B(J,I)=BJ+CG*BK
        A(K,I)=AK+CA*AJ
  150   B(K,I)=BK+CA*BJ
  160   IF(JP1-KM1) 170,170,190
  170   DO 180 I=JP1,KM1
        AJ=A(J,I)
        BJ=B(J,I)
        AK=A(I,K)
        BK=B(I,K)
        A(J,I)=AJ+CG*AK
        B(J,I)=BJ+CG*BK
        A(I,K)=AK+CA*AJ
  180   B(I,K)=BK+CA*BJ
  190   AK=A(K,K)
        BK=B(K,K)
        A(K,K)=AK+2.*CA*A(J,K)+CA*CA*A(J,J)
        B(K,K)=BK+2.*CA*B(J,K)+CA*CA*B(J,J)
        A(J,J)=A(J,J)+2.*CG*A(J,K)+CG*CG*AK
        B(J,J)=B(J,J)+2.*CG*B(J,K)+CG*CG*BK
        A(J,K)=0.
        B(J,K)=0.
C----------------------------------------------------------------------
C              UPDATE THE EIGENVECTOR MATRIX AFTER EACH ROTATION
C----------------------------------------------------------------------
        DO 200 I=1,N
        XJ=X(I,J)
        XK=X(I,K)
        X(I,J)=XJ+CG*XK
  200   X(I,K)=XK+CA*XJ
  210   CONTINUE
C----------------------------------------------------------------------
C              UPDATE THE EIGENVALUES AFTER EACH SWEEP
C----------------------------------------------------------------------
        DO 220 I=1,N
        IF(A(I,I).GT.0.AND.B(I,I).GT.0) GOTO 220
        WRITE(*,2020)
        STOP
  220   EIGV(I)=A(I,I)/B(I,I)
C----------------------------------------------------------------------
C             CHECK FOR CONVERGENCE
C----------------------------------------------------------------------
  230   DO 240 I=1,N
        TOL=RTOL*D(I)
        DIF=DABS(EIGV(I)-D(I))
        IF(DIF.GT.TOL) GOTO 400
  240   CONTINUE
C----------------------------------------------------------------------
C  CHECK ALL OFF-DIAGONAL ELEMENTS TO SET IF ANOTHER SWEEP IS REQUIRED
C----------------------------------------------------------------------
        EPS=RTOL**2
        DO 250 J=1,NR
        JJ=J+1
        DO 250 K=JJ,N
        EPSA=(A(J,K)*A(J,K))/(A(J,J)*A(K,K))
        EPSB=(B(J,K)*B(J,K))/(B(J,J)*B(K,K))
        IF(EPSA.LT.EPS.AND.EPSB.LT.EPS) GOTO 250
        GOTO 400
  250   CONTINUE
C
C   FILL OUT BOTTOM TRIANGLE OF RESULTANT MATRICES AND SCALE EIGENVERTORS
C
  255   DO 260 I=1,N
        DO 260 J=1,N
        A(J,I)=A(I,J)
  260   B(J,I)=B(I,J)
        DO 270 J=1,N
        BB=DSQRT(B(J,J))
        DO 270 K=1,N
  270   X(K,J)=X(K,J)/BB
        DO 280 I=1,N
        D(I)=EIGV(I)
        DO 280 J=1,N
  280   A(I,J)=X(I,J)
        DO 290 I=1,N
        M(I)=N
        DO 290 J=1,N
  290   IF(D(I).LT.D(J)) M(I)=M(I)-1
        DO 300 I=1,N
        K=M(I)
        EIGV(K)=D(I)
        DO 300 J=1,N
  300   X(J,K)=A(J,I)
        RETURN
C----------------------------------------------------------------------
C               UPDATE  D  MATRIX AND START NEW SWEEP , IF ALLOWED
C----------------------------------------------------------------------
  400   DO 500 I=1,N
  500   D(I)=EIGV(I)
        IF(NSWEEP.LT.500) GOTO 40
        GOTO 255
        END
C======================================================================

[ 本帖最后由 欧阳中华 于 2008-5-24 09:17 编辑 ]

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2008-5-24 11:43 | 显示全部楼层
麻烦问一下楼主:这本书哪里能弄到?我目前正在从事这方面的工作。
 楼主| 发表于 2008-5-24 14:34 | 显示全部楼层
.
    学校图书馆. ..
发表于 2008-5-26 09:21 | 显示全部楼层

感谢

非常感谢!以后一定把这种共享精神发扬下去
发表于 2010-4-6 20:56 | 显示全部楼层
good
:lol :lol :lol
 楼主| 发表于 2016-8-23 14:20 | 显示全部楼层
laneliu 发表于 2008-5-24 11:43
麻烦问一下楼主:这本书哪里能弄到?我目前正在从事这方面的工作。

.
    买了一本原版的,也有PDF文档,非常经典的专著. .. ..

点评

差不多十年前从师兄那里抢了一本大工某院士的中译本,老实说读着别扭,还不如后来的Finite Element Procedure原版 国内最近有人翻译了Finite Element Procedure的第二版,价钱贵没买。。。 btw,这个洛阳铲够深的  详情 回复 发表于 2016-8-27 19:29
发表于 2016-8-27 19:29 | 显示全部楼层
欧阳中华 发表于 2016-8-23 14:20
.
    买了一本原版的,也有PDF文档,非常经典的专著. .. ..

差不多十年前从师兄那里抢了一本大工某院士的中译本,老实说读着别扭,还不如后来的Finite Element Procedure原版
国内最近有人翻译了Finite Element Procedure的第二版,价钱贵没买。。。
btw,这个洛阳铲够深的
发表于 2016-9-1 16:16 | 显示全部楼层
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-1 03:53 , Processed in 0.072136 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表