
C     TITLE:  POWER METHOD FOR EIGENVALUE PROBLEM
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER T
      PARAMETER (MAXN=10)
      DIMENSION A(MAXN,MAXN), Z(MAXN), W(MAXN)
      DATA IOUT/8/
C
      OPEN(IOUT,FILE='power.out')
5     PRINT *, ' WHAT IS THE ORDER OF THE MATRIX A?'
      PRINT *, ' GIVE 0 TO STOP.'
      READ *, N
      IF(N .EQ. 0) THEN
          CLOSE(IOUT)
          STOP
      END IF
      PRINT *, ' GIVE THE MATRIX A, ONE ROW AT A TIME.'
      DO 7 I=1,N
        PRINT *, ' GIVE ROW',I
7       READ *, (A(I,J), J=1,N)
      PRINT *, ' GIVE THE INITIAL GUESS FOR THE EIGENVECTOR Z.'
      READ *, (W(I), I=1,N)
      PRINT *, ' WHICH COMPONENT FOR THE EIGENVALUE COMPUTATION?'
      READ *, ICOMP
C
      WRITE(IOUT,*) 'THE COEFFICIENT MATRIX'
      DO 8 I=1,N
        WRITE(IOUT,1005) I
        WRITE(IOUT,1004) (A(I,J),J=1,N)
8       CONTINUE
      WRITE(IOUT,*) ' '
      WRITE(IOUT,*) 'THE INITIAL GUESS'
      WRITE(IOUT,1004) (W(I),I=1,N)
      WRITE(IOUT,*) ' '
C
C     INITIALIZE.
      WMAX = 0.0D0
      DO 10 I=1,N
        IF(ABS(W(I)) .GT. WMAX) THEN
            IMAX = I
            WMAX = ABS(W(I))
        END IF
10      CONTINUE
      DIVISOR = W(IMAX)
      DO 20 I=1,N
20      Z(I) = W(I)/DIVISOR
      LOOP = 0
C
C     BEGIN MAIN LOOP.
25    DO 40 I=1,N
        SUM = 0.0D0
        DO 30 J=1,N
30        SUM = SUM + A(I,J)*Z(J)
40      W(I) = SUM
      WMAX = 0.0D0
      DO 50 I=1,N
        IF(ABS(W(I)) .GT. WMAX) THEN
            IMAX = I
            WMAX = ABS(W(I))
        END IF
50      CONTINUE
      DIVISOR = W(IMAX)
      FLAM = W(ICOMP)/Z(ICOMP)
      DO 60 I=1,N
60      Z(I) = W(I)/DIVISOR
C
C     PRINT RESULTS.
      PRINT *, 'APPROXIMATE EIGENVECTOR.'
      PRINT 1000, (Z(I), I=1,N)
      WRITE(IOUT,*) 'APPROXIMATE EIGENVECTOR.'
      WRITE(IOUT,1000) (Z(I), I=1,N)
1000  FORMAT(6D20.10)
      IF(LOOP .GT. 1) THEN
          DIFF  =FLAM - PASTLAM
          RATIO = DIFF/PASTDF
          PRINT 1001, IMAX, FLAM, DIFF, RATIO
          WRITE(IOUT,1001) IMAX, FLAM, DIFF, RATIO
1001      FORMAT(' COMPONENT=',I2,4X,'LAMBDA=',1PD16.9,4X,
     *           'DIFF =',D9.2,3X,'RATIO =',D9.2)
      ELSE
          IF(LOOP .GT. 0) THEN
              DIFF = FLAM - PASTLAM
              PRINT 1002, IMAX,FLAM,DIFF
              WRITE(IOUT,1002) IMAX,FLAM,DIFF
1002          FORMAT(' COMPONENT=',I2,4X,'LAMBDA=',1PD16.9,4X,
     *               'DIFF =',D10.3)
          ELSE
              PRINT 1003, IMAX,FLAM
              WRITE(IOUT,1003) IMAX,FLAM
1003          FORMAT(' COMPONENT=',I2,4X,'LAMBDA=',1PD16.9)
          END IF
      END IF
      LOOP = LOOP + 1
      PASTDF = DIFF
      PASTLAM = FLAM
C
C     ANOTHER LOOP?
      PRINT *, ' ANOTHER LOOP? (Y/N)'
      READ(*,'(A1)') T
      IF((T .EQ. 'Y') .OR. (T .EQ. 'y')) THEN
          GO TO 25
      ELSE
          GO TO 5
      END IF
1004  FORMAT(1PE19.9,3E19.9)
1005  FORMAT(' ROW ',I2)
      END
