PARAMETER (NP=30)
    IMPLICIT REAL*8 (a-h, o-z)
    DOUBLE PRECISION A(NP,NP),GK(NP,NP),GM(NP,NP)

    OPEN(5,FILE=\' \')
    OPEN(10,FILE=\' \')
    READ(5,*) OPT    
    IF(OPT .EQ. 1 .OR. OPT .EQ. 2) THEN
        READ(5,*) N, TOL
      READ(5,*) ((A(I,J), J=1,N), I=1,N)
      WRITE (10,11)N,N
      WRITE (10,12)
      WRITE (10,13) ((I,J,A(I,J),J=1,N),I=1,N)
    ELSEIF(OPT .EQ. 3 .OR. OPT .EQ. 4) THEN
      READ(5,*) N, TOL
      READ(5,*) ((GK(I,J), J=1,N), I=1,N)
      READ(5,*) ((GM(I,J), J=1,N), I=1,N)
      WRITE(10,14) N,N,N,N
      WRITE(10,15)
      WRITE(10,16) ((I,J,GK(I,J),J=1,N),I=1,N)
      WRITE(10,17)
      WRITE(10,18) ((I,J,GM(I,J),J=1,N),I=1,N)
    ENDIF

   11    FORMAT (50(\'=\'),/,5X,\'Eigen Problem\',/,
     1  10X,\'- Power Method & Inverse Method\',/,50(\'=\'),///,
     2  5X,\'[A][X]=λ[X]\',//,
     3  5X,\'where [A] is the Coefficient Matrix : size=\',I3,\' X\',I3,/,
    4  11X,\'[X] is the Eigen Vector\',/,         
     5  11X,\' λ is the Eigen Value\',//)
   12    FORMAT (3X,9(\'*\'),\'  MATRIX [A]  \',9(\'*\'),/)
   13    FORMAT (5X,\'A(\',I3,\',\',I3,\') = \',F15.5)

   14    FORMAT (50(\'=\'),/,5X,\'Eigen Problem\',/,
     1  10X,\'- Power Method & Inverse Method\',/,50(\'=\'),///,
     2  5X,\'[K]Φ=λ[M]Φ\',//,
     3  5X,\'where [K] is the Stiffness Matrix : size=\',I3,\' X\',I3,/,
    4  5X,\'      [M] is the Mass Matrix : size=\',I3,\' X\',I3,/,
     4  11X,\' Φ is the Eigen Vector\',/,         
     5  11X,\' λ is the Eigen Value\',//)
   15    FORMAT (3X,9(\'*\'),\'  MATRIX [K]  \',9(\'*\'),/)
   16    FORMAT (5X,\'K(\',I3,\',\',I3,\') = \',F15.5)
   17    FORMAT (/,3X,9(\'*\'),\'  MATRIX [M]  \',9(\'*\'),/)
   18    FORMAT (5X,\'M(\',I3,\',\',I3,\') = \',F15.5)

    IF(OPT .EQ. 1) THEN
       CALL POWERS (A,N,NP,TOL)
    ELSEIF(OPT .EQ. 2) THEN
     CALL INVERS (A,N,NP,TOL)
    ELSEIF(OPT .EQ. 3) THEN
     CALL POWERG (GK,GM,N,NP,TOL)
    ELSEIF(OPT .EQ. 4) THEN
     CALL INVERG (GK,GM,N,NP,TOL)  
    ENDIF
        
    PRINT*,\'EINGEN PROBLEM Program is terminated\'
    PRINT*,\'Check Your OUT File !!!\'
     STOP
    END

C STANDARD POWER METHOD SUBROUTINE

    SUBROUTINE POWERS (A,N,NP,TOL)
    IMPLICIT REAL*8 (a-h, o-z)
     DOUBLE PRECISION A(NP,NP),X(NP),Y(NP),Y1(NP),DFF(NP)
    DO J=1,N
      X(J)=1.0
    ENDDO

  100    DFF=0.0
    DFF1=0.0
    SQ=0.0
    DO K=1,N
      Q=0.0
      DO L=1,N
        Q=Q+A(K,L)*X(L)
      ENDDO
      Y(K)=Q
      SQ=SQ+Q**2
    ENDDO
    EVL=SQRT(SQ)
    DO K=1,N
      Y1(K)=Y(K)/EVL
    DFF(k)=(Y1(K)-X(K))
    X(K)=Y1(K)
    ENDDO

    CALL NORM(DFF, DFF1, NP, N)
    IF(DFF1 .LE. TOL) THEN
      GOTO 199
    ELSE
      GOTO 100
    ENDIF

  199    WRITE (10,111) EVL
  111    FORMAT (/,3x,\'STANDARD POWER METHOD\',//,
     1      3X,\'***** EIGEN VALUE *****\',/,5X,F10.4)
    WRITE (10,112)
  112    FORMAT (/,3X,\'***** EIGEN VECTOR *****\')
      DO J=1,N
      WRITE (10,113) X(J)
  113    FORMAT (3X,F10.4)
    ENDDO
    RETURN
    END

C STANDARD INVERSE METHOD SOURUOTINE

    SUBROUTINE INVERS (A,N,NP,TOL)
    IMPLICIT REAL*8 (a-h, o-z)
     DOUBLE PRECISION A(NP,NP),X(NP),Y(NP),Y1(NP),AIN(NP,NP),DFF(NP)
    DO J=1,N
      X(J)=1.0
    ENDDO
C    CALL SHIFT (A,N,NP,SH)
c      DO I=1,SH           ! SHIFT ITERATION FOR SHIFT NUMBER
    
      DO K=1,N
      DO L=1,N
        AIN(K,L)=A(K,L)
      ENDDO
    ENDDO
      CALL LUMINV (AIN,N,NP)
  200    DFF=0.0
    DFF1=0.0
    SQ=0.0
    DO K=1,N
      Q=0.0
      DO L=1,N
        Q=Q+AIN(K,L)*X(L)
      ENDDO
      Y(K)=Q
      SQ=SQ+Q**2
    ENDDO
    EVLINV=SQRT(SQ)
    EVL=1.0/EVLINV
    DO K=1,N
      Y1(K)=Y(K)/EVLINV
    DFF(K)=(Y1(K)-X(K))
    X(K)=Y1(K)
    ENDDO

    CALL NORM(DFF, DFF1, NP, N)
    IF(DFF1 .LE. TOL) THEN
      GOTO 299
    ELSE
      GOTO 200
    ENDIF

c  299    IF(SH .GE. 2) THEN
c      EVL=EVL+RQSH
c      WRITE(10,214) RQSH
c  214      FORMAT(/,3X,\'Reyleigh Quotient Shift = \',F10.4)
c    ENDIF

  299    WRITE (10,211) EVL
  211    FORMAT (/,3x,\'STANDARD INVERSE METHOD\',//,
     1      3X,\'***** EIGEN VALUE *****\',/,5X,F10.4)
    WRITE (10,212)
  212    FORMAT (/,3X,\'***** EIGEN VECTOR *****\')
      DO J=1,N
      WRITE (10,213) X(J)
  213    FORMAT (5X,F10.4)
    ENDDO
c    IF(SH .GE. 2) THEN
c     CALL SHIFTI (A,N,NP,X,RQSH)
c    ENDIF
c
c    ENDDO              ! SHIFT ITERATION FOR SHIFT NUMBER
    RETURN
    END



뭘로 돌려야하는거 ???
이게 소스는 맞는거냐;