PROGRAM POWERMETHOD2
C=======================================================================
C       THIS PROGRAM FINDS ALL EIGENVALUES AND VECTORS FOR A GIVEN [A]
C       USING POWER METHOD.
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( MXN=10, MXITERA=100 )
      DIMENSION A(MXN,MXN), X(MXN),Y(MXN),EIGENVEC(MXN,MXN),EIGEN(MXN),
     *          B(MXN,MXN), AA(MXN,MXN)
C=======================================================================
      EPS = 1.D-13
C=======================================================================
      OPEN (1,FILE='POWER.OUT', STATUS='UNKNOWN')
C=======================================================================
C------ INITIALIZATION
      N = 3
      A(1,1) = 3.D0
      A(1,2) = 2.D0
      A(1,3) = 1.D0
      A(2,1) = 2.D0
      A(2,2) = 2.D0
      A(2,3) = 1.D0
      A(3,1) = 1.D0
      A(3,2) = 1.D0
      A(3,3) = 1.D0
      WRITE (1,*) 'GIVEN MATRIX [A]'
      DO I = 1 , N
      WRITE (1,*) ( A(I,J),J = 1 , N )
      END DO
C-------- SAVE [A]
      DO I = 1 , N
      DO J = 1 , N
      AA(I,J) = A(I,J)
      END DO
      END DO
C
      WRITE (1,'(16H MODE EIGENVALUE, 99(I2,1H)))') ( I,I=1,N)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      DO MODE = 1 , N
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      IF ( MODE .GT. 1 ) THEN
      DO I = 1 , N
      DO J = 1 , N
      A(I,J) = A(I,J) -
     * EIGEN(MODE-1)*EIGENVEC(I,MODE-1)*EIGENVEC(J,MODE-1)
      END DO
      END DO
      END IF
C
      EIGEN(MODE) = 1.D0
      DO I = 1 , N
      X(I) = 0.D0
      END DO
      X(1) = 1.D0
      WRITE (*,*) MODE,EIGEN(MODE), (X(I),I=1,N)
      WRITE (1,*) MODE,EIGEN(MODE), (X(I),I=1,N)
      DIFF = 1.D0
      NITERA = 0
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      DO WHILE ( DIFF .GT. EPS .AND. NITERA .LE. MXITERA )
      NITERA = NITERA + 1
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C-------- COMPUTATION OF [A]{X}={Y}
      DO I = 1 , N
      SUM = 0.D0
      DO J = 1 , N
      SUM = SUM + A(I,J)*X(J)
      END DO
      Y(I) = SUM
      END DO
C-------- COMPUTATION OF EIGENVALUE
      SUM = 0.D0
      DO I = 1 , N
      SUM = SUM + Y(I)*Y(I)
      END DO
      EIGEN(MODE) = DSQRT(SUM)
C-------- COMPUTATION OF ERROR
      DIFF = 0.D0
      DO I = 1 , N
      XNEW = Y(I) / EIGEN(MODE)
      DIFF = DIFF + DABS( X(I) - XNEW )
      X(I) = XNEW
      END DO
      DIFF = DIFF / N
C
      WRITE (*,*) MODE, EIGEN(MODE), (X(I),I=1,N)
      WRITE (1,*) MODE, EIGEN(MODE), (X(I),I=1,N)
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      END DO
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      DO I = 1 , N
      EIGENVEC(I,MODE) = X(I)
      END DO
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      END DO
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C=======================================================================
      WRITE (1,*) 'EIGEN VECTOR DOT PRODUCT'
      WRITE (1,*) 'MODE-I MODE-J DOTPRODUCT'
      DO I = 1 , N-1
      DO J = I+1,N
      DOTPRDCT = 0.D0
      DO K = 1 , N
      DOTPRDCT = DOTPRDCT + EIGENVEC(K,I)*EIGENVEC(K,J)
      END DO
      WRITE (1,*) I, J, DOTPRDCT
      END DO
      END DO
C=======================================================================
C------ REEVALUATION OF [A] USING EIGENVALUES AND EIGENVECTORS
      WRITE(1,*) 'REFORMING [A] USING EIGENVALUES AND VECTORS'
      DO I = 1 , N
      DO J = 1 , N
      B(I,J) = 0.D0
      END DO
      END DO
      DO MODE = 1 , N
      DO I = 1 , N
      DO J = 1 , N
      B(I,J) = B(I,J) + EIGEN(MODE)*EIGENVEC(I,MODE)*EIGENVEC(J,MODE)
      END DO
      END DO
      END DO
      DO I = 1 , N
      WRITE (1,*) ( B(I,J),J = 1 , N )
      END DO
C=======================================================================
C------- PRODUCT OF [A][U]
      DO I = 1 , N
      DO J = 1 , N
      SUM = 0.D0
      DO K = 1 , N
      SUM = SUM + AA(I,K)*EIGENVEC(K,J)
      END DO
      B(I,J) = SUM
      END DO
      END DO
C------- PRODUCT OF [U]TRANSPOSE[A][U]
      DO I = 1 , N
      DO J = 1 , N
      A(I,J) = 0.D0
      END DO
      END DO
      DO I = 1 , N
      DO J = 1 , N
      DO K = 1 , N
      A(I,J) = A(I,J) + EIGENVEC(K,I)*B(K,J)
      END DO
      END DO
      END DO
      WRITE (1,*) 'PRODUCT OF [U]TRANSPOSE[A][U]'
      DO I = 1 , N
      WRITE (1,*) ( A(I,J),J = 1 , N )
      END DO
C=======================================================================
      STOP 'NORMAL TERMINATION'
      END