Jacobi 广义对称矩阵特征值计算子程序
.这里给出的是满阵形式存储的Jacobi法广义对称矩阵特征值计算的子程序代码，来自《 NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS（K. J. BatheE. L. Wilson）》：
C......................................................................
C .
CPROGRAM 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 EIGENVALUESON 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 byK. J. BatheE. 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 .
CTHIS 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,' ***** Errorsolution ***** ')
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--------------------------------------------------------------------------
CCHECK 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--------------------------------------------------------------------------
CIF 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---------------------------------------------------------------------------
CPERPORM 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----------------------------------------------------------------------
CCHECK 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 UPDATEDMATRIX 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 编辑 ] 麻烦问一下楼主：这本书哪里能弄到？我目前正在从事这方面的工作。 .
学校图书馆. ..
感谢
非常感谢！以后一定把这种共享精神发扬下去 good:lol :lol :lol laneliu 发表于 2008-5-24 11:43
麻烦问一下楼主：这本书哪里能弄到？我目前正在从事这方面的工作。
.
买了一本原版的，也有PDF文档，非常经典的专著. .. .. 欧阳中华 发表于 2016-8-23 14:20
.
买了一本原版的，也有PDF文档，非常经典的专著. .. ..
差不多十年前从师兄那里抢了一本大工某院士的中译本，老实说读着别扭，还不如后来的Finite Element Procedure原版
国内最近有人翻译了Finite Element Procedure的第二版，价钱贵没买。。。
btw，这个洛阳铲够深的 {:{39}:}
页:
[1]