home *** CD-ROM | disk | FTP | other *** search
- C ******************************************************************
- C * *
- C * E I G E N *
- C * *
- C ******************************************************************
- PROGRAM EIGEN
- IMPLICIT INTEGER*4(I-N)
- IMPLICIT REAL*8 (A-H,O-Z)
- COMMON MAVAIL,IA(30000)
- MAVAIL=30000
- C ... SAMPLE: EIGEN PROBLEM SOLVER
- PRINT *,' '
- PRINT *,'EIGEN PROBLEM SOLVER'
- PRINT *,' '
- PRINT *,'ENTER EQUATION NUMBER '
- READ *, N
- PRINT *,'NOW ENTER STIFFNESS MATRIX & MASS MATRIX'
- PRINT *,' '
- C ... TEST IN-CORE DATA MANAGEMENT
- CALL DBOPEN(1,'TSTDAT.DAT','NEW')
- PRINT *,'DATABASE OPENED.'
- CALL DEFINE(1,'STIF',0,1,N,N,0,NA)
- CALL DEFINE(1,'MASS',0,1,N,N,0,NB)
- CALL DEFINE(1,'VCTR',0,1,N,N,0,NC)
- CALL DEFINE(1,'EIGN',0,1,N,1,0,ND)
- CALL DEFINE(1,'WORK',0,1,N,1,0,NE)
- print *,'DEFINE MATRICES DONE.'
- CALL MATINP(1,'STIF')
- CALL MATINP(1,'MASS')
- CALL MATOUT(1,'STIF')
- CALL MATOUT(1,'MASS')
- RTOL = 1.0D-12
- CALL JACOBI(IA(NA),IA(NB),IA(NC),IA(ND),IA(NE),N,RTOL,15,0,6)
- CALL MATOUT(1,'EIGN')
- CALL MATOUT(1,'VCTR')
- CALL DIR(0)
- CALL DBCLOS(1,'DELETE')
- STOP 'DONE.'
- END
- C ******************************************************************
- C * *
- C * J A C O B I *
- C * *
- C ******************************************************************
- SUBROUTINE JACOBI (A,B,X,EIGV,D,N,RTOL,NSMAX,IFPR,IOUT)
- IMPLICIT INTEGER*4(I-N)
- IMPLICIT REAL*8 (A-H,O-Z)
- C
- DIMENSION A(N,N),B(N,N),X(N,N),EIGV(N),D(N)
- C
- C TO SOLVE THE GENERALIZED EIGENPROBLEM USING THE GENERALIZED
- C JACOBI ITERATION
- C
- C INPUT VARIABLES :
- C A(N,N) = STIFFNESS MATRIX (ASSUME POSITIVE DEFINITE)
- C B(N,N) = MASS MATRIX (ASSUME POSITIVE DEFINITE)
- C D(N) = WORKING VAECTOR
- C N = ORDER OF MATRIX A AND B
- C RTOL = CONVERGENCE TOLERANCE (SET TO 1.E-12)
- C NSMAX = MAXIMUM NUMBER OF SWEEP ALLOWED (SET TO 15)
- C IFPR = FLAG FOR PRINTING DURING ITERATION
- C ( EQ.0 NO PRINTING, EQ.1 PRINT INTERMEDIATE RESULTS)
- C IOUT = OUTPUT LUN
- C OUTPUT VARIABLES :
- C A(N,N) = DIAGONALIZED STIFFNESS MATRIX
- C B(N,N) = DIAGONALIZED MASS MATRIX
- C X(N,N) = EIGENVECTORS STORED COLUMNWISE
- C EIGV(N) = EIGENVALUES
- C
-
- C INITIALIZE EIGENVALUE AND EIGENVECTOR MATRICES
- DO 10 I=1,N
- IF (A(I,I).LE.0..OR.B(I,I).LT.0.) THEN
- WRITE(IOUT,2020)
- STOP 'JACOBI'
- END IF
- 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.0
- IF(N.EQ.1) RETURN
- C INITIALIZE SWEEP COUNTER AND EIGEN ITERATION
- NSWEEP=0
- NR=N-1
- C
- C WE START ITERATION
- 40 NSWEEP=NSWEEP+1
- IF(IFPR.EQ.1) WRITE(IOUT,1000) NSWEEP
- C CHECK IF PRESENT OFF-DIAGONAL ELEMENT IS LARGE ENOUGH TO
- C REQUIRE ZEROING
- EPS=(0.01**NSWEEP)**2
- DO 210 J=1,NR
- JJ=J+1
- DO 210 K=JJ,N
- TT=A(J,K)*A(J,K)
- TB=A(J,J)*A(K,K)
- EPTOLA=DABS(TT/TB)
- TT=B(J,K)*B(J,K)
- TB=B(J,J)*B(K,K)
- EPTOLB=TT/TB
- IF ((EPTOLA.LT.EPS).AND.(EPTOLB.LT.EPS)) GO TO 210
- C IF ZEROING IS REQUIRED, CALCULATE THE ROTATION MATRIX ELEMENT CA, CG
- 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.0*AKK*AJJ)/4.0
- IF (CHECK) 60,60,70
- 60 WRITE (IOUT,1004) CHECK
- STOP
- 70 SQCH=DSQRT(CHECK)
- D1=AB/2.0+SQCH
- D2=AB/2.0-SQCH
- DEN=D1
- IF (DABS(D2).GT.DABS(D1)) DEN=D2
- IF (DEN) 90,80,90
- 80 CA=0.
- CG=-A(J,K)/A(K,K)
- GO TO 100
- 90 CA=AKK/DEN
- CG=-AJJ/DEN
- C
- C WE PERFORM THE GENERALIZED ROTATION
- 100 IF (N-2) 101,180,101
- 101 JP1=J+1
- JM1=J-1
- KP1=K+1
- KM1=K-1
- C
- IF (JM1-1) 120,110,110
- 110 DO 105 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
- 105 B(I,K)=BK+CA*BJ
- C
- 120 IF (KP1-N) 130,130,140
- 130 DO 125 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
- 125 B(K,I)=BK+CA*BJ
- C
- 140 IF (JP1-KM1) 150,150,180
- 150 DO 160 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
- 160 B(I,K)=BK+CA*BJ
- 180 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.0
- B(J,K)=0.0
- C
- C UPDATE EIGENVECTORS
- DO 190 I=1,N
- XJ=X(I,J)
- XK=X(I,K)
- X(I,J)=XJ+CG*XK
- 190 X(I,K)=XK+CA*XJ
- C
- 210 CONTINUE
- C
- DO 220 I=1,N
- 220 EIGV(I)=A(I,I)/B(I,I)
- IF(IFPR.EQ.0) GO TO 227
- WRITE (IOUT,1005)
- WRITE (IOUT,1002) (EIGV(I),I=1,N)
- 227 CONTINUE
- C
- C CHECK FOR CONVERGENCE
- DO 240 I=1,N
- TOL=RTOL*D(I)
- DIF=DABS(EIGV(I)-D(I))
- IF(DIF.GT.TOL) GO TO 300
- 240 CONTINUE
- C
- C CHECK IF ALL OFF-DIAG ELEMENTS ARE SATISFACTORILY SMALL
- EPS=RTOL**2
- DO 260 J=1,NR
- JJ=J+1
- DO 260 K=JJ,N
- TT=A(J,K)*A(J,K)
- TB=A(J,J)*A(K,K)
- EPSA=DABS(TT/TB)
- TT=B(J,K)*B(J,K)
- TB=B(J,J)*B(K,K)
- EPSB=TT/TB
- IF ((EPSA.LT.EPS).AND.(EPSB.LT.EPS)) GO TO 260
- GO TO 300
- 260 CONTINUE
- C
- DO 310 I=1,N
- DO 310 J=I,N
- B(J,I)=B(I,J)
- 310 A(J,I)=A(I,J)
- RETURN
- C
- 300 DO 320 I=1,N
- 320 D(I)=EIGV(I)
- IF (NSWEEP.LT.NSMAX) GO TO 40
- DO 330 I=1,N
- DO 330 J=I,N
- B(J,I)=B(I,J)
- 330 A(J,I)=A(I,J)
- RETURN
- C
- 1000 FORMAT (27H0SWEEP NUMBER IN *JACOBI* =, I4)
- 1002 FORMAT (1H ,12E11.4)
- 1004 FORMAT (37H0***ERROR SOLUTION STOP IN *JACOBI*, / 12X,
- 1 8HCHECK = , E20.14 / 1X)
- 1005 FORMAT (36H0CURRENT EIGENVALUES IN *JACOBI* ARE, / 1X)
- 2020 FORMAT ('0**ERROR** MATRIX NOT POSITIVE DEFINITE IN JACOBI'/)
- C
- END
