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
