C***********************************************************************
C                                                                      *
C               INVERSE POWER METHOD                                   *
C                                                                      *
C***********************************************************************
C
C
C
C     TO APPROXIMATE AN EIGENVALUE AND AN ASSOCIATED EIGENVECTOR OF THE
C     N BY N MATRIX A GIVEN A NONZERO VECTOR X:
C
C     INPUT:   DIMENSION N; MATRIX A; VECTOR X; TOLERANCE TOL;
C              MAXIMUM NUMBER OF ITERATIONS N.
C
C     OUTPUT:  APPROXIMATE EIGENVALUE MU; APPROXIMATE EIGENVECTOR
C              X OR A MESSAGE THAT THE MAXIMUM NUMBER OF ITERATIONS
C              WAS EXCEEDED.
C
      DIMENSION A(3,3),X(3),Y(3),NROW(3),B(3)
      CHARACTER NAME*14,NAME1*14,AA*1
      INTEGER INP,OUP,FLAG
      LOGICAL OK
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Inverse Power Method.'
      WRITE(6,*) 'The array will be input from a text file in the'
      WRITE(6,*) ' order: A(1,1), A(1,2), ..., A(1,N), A(2,1),'
      WRITE(6,*) ' A(2,2), ..., A(2,N)..., A(N,1), A(N,2),'
      WRITE(6,*) ' ..., A(N,N) '
      WRITE(6,*) 'Place as many entries as desired on each line,'
      WRITE(6,*) ' but separate entries with at least one blank.'
      WRITE(6,*) 'The initial approximation should follow in the'
      WRITE(6,*) 'same format.'
      OK = .FALSE.
      WRITE(6,*) 'Has the input file been created?'
      WRITE(6,*) 'Enter Y or N - letter within quotes '
      WRITE(6,*) ' '
      READ(5,*)  AA
      IF (( AA .EQ. 'Y' ) .OR.( AA .EQ. 'y' )) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext  contained in quotes'
         WRITE(6,*) 'as example:   ''A:DATA.DTA'' '
         WRITE(6,*) ' '
         READ(5,*)  NAME
         INP = 4
         OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL')
         OK = .FALSE.
19       IF (OK) GOTO 111
         WRITE(6,*) 'Input the dimension n.'
         WRITE(6,*)
         READ(5,*) N
         IF (N .GT. 0) THEN
            READ(INP,*) ((A(I,J), J=1,N),I=1,N)
            READ(INP,*) (X(I),I=1,N)
            OK = .TRUE.
            CLOSE(UNIT=INP)
         ELSE
            WRITE(6,*) 'The number must be a positive integer'
         ENDIF
         GOTO 19
111      OK = .FALSE.
112      IF (OK) GOTO 13
         WRITE(6,*) 'Input the tolerance.'
         WRITE(6,*) ' '
         READ(5,*) TOL
         IF (TOL .GT. 0.0) THEN
            OK = .TRUE.
         ELSE
            WRITE(6,*) 'Tolerance must be positive.'
         ENDIF
         GOTO 112
13       OK = .FALSE.
14       IF (OK) GOTO 15
         WRITE(6,*) 'Input maximum number of iterations.'
         WRITE(6,*)
         READ(5,*) NN
         IF (NN .GT. 0) THEN
            OK = .TRUE.
         ELSE
            WRITE(6,*) 'Number must be a positive integer.'
         ENDIF
         GOTO 14
      ELSE
         WRITE(6,*) 'The program will end so the input file can '
         WRITE(6,*) 'be created. '
         OK = .FALSE.
      ENDIF
15    IF(.NOT. OK) GOTO 400
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'INVERSE POWER METHOD'
C     STEP 1
C     Q COULD BE INPUT INSTEAD OF COMPUTED BY DELETING THE NEXT 7 STEPS
      Q=0
      S = 0
      DO 50 I=1,N
           S= S + X(I)*X(I)
           DO 50 J=1,N
50         Q = Q + A(I,J)*X(I)*X(J)
      Q = Q/S
C     STEP 2
      K =1
      WRITE(OUP,4)
      WRITE(OUP,5) ((A(I,J),J=1,N),I=1,N)
      WRITE(OUP,*) 'INPUT VECTOR'
      J=0
      WRITE(OUP,6) J,(X(I),I=1,N)
      WRITE(OUP,7) Q
C     FORM MATRIX A - Q*I
      DO 10 I=1,N
10         A(I,I) = A(I,I)-Q
C     CALL SUBROUTINE TO COMPUTE MULTIPLIERS M(I,J) AND UPPER TRIANGULAR
C     MATRIX FOR MATRIX A USING GAUSS ELIMINATION WITH MAXIMAL COLUMN
C     PIVOTING--NROW HOLDS THE ORDERING OF ROWS FOR INTERCHANGES
      CALL MULTIP(N,A,NROW)
C     STEP 3
      LP = 1
      DO 20 I=2,N
20         IF(ABS(X(I)).GT.ABS(X(LP))) LP =I
C     STEP 4
      AMAX = X(LP)
      DO 30 I=1,N
30         X(I) = X(I)/AMAX
C     STEP 5
100   IF (K.GT.NN) GOTO 200
C          STEP 6
           DO 35 I=1,N
35         B(I)=X(I)
C          SUBROUTINE SOLVE RETURNS THE SOLUTION OF (A-Q*I)Y=B IN Y
           CALL SOLVE (N,A,NROW,B,Y)
C          STEP 7
           YMU = Y(LP)
C          STEP 8 AND 9
           LP=1
           DO 40 I=2,N
40              IF(ABS(Y(I)).GT.ABS(Y(LP))) LP=I
           AMAX = Y(LP)
           ERR=0.0
           DO 60 I=1,N
                T=Y(I)/AMAX
                IF(ABS(X(I)-T).GT.ERR) ERR=ABS(X(I)-T)
60              X(I)=T
           WRITE(OUP,*) 'ITER. NO. AND VECTOR'
           WRITE(OUP,6) K,(X(I),I=1,N)
           WRITE(OUP,8) YMU,ERR
C          STEP 10
           IF (ERR.LT.TOL) THEN
                YMU=1/YMU+Q
                WRITE(OUP,3) YMU
C               PROCEDURE COMPLETED SUCCESSFULLY
                GOTO 400
           END IF
C          STEP 11
           K = K+1
      GOTO 100
C     STEP 12
200   CONTINUE
      WRITE (OUP,9)
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
3     FORMAT(1X,'APPROXIMATE EIGENVALUE IS ',E15.8)
4     FORMAT(3X,' MATRIX A ')
5     FORMAT(3(1X,E15.8))
6     FORMAT(1X,I2,3(1X,E15.8))
8     FORMAT(1X,'APPROX. EIGENVALUE',1X,E15.8,' ERROR',E15.8)
7     FORMAT(3X,'Q = ',E15.8)
9     FORMAT(1X,'NO COVERGENCE')
      END
C
      SUBROUTINE MULTIP (N,A,NROW)
C
C     SUBROUTINE MULTIP COMPUTES THE MULTIPLIERS AND ROW ORDERING FOR
C     GAUSS ELIMINATION WITH MAXIMAL COLUMN PIVOTING-THE MULTIPLIERS
C     AND UPPER TRIANGULAR FORM ARE RETURNED IN A AND THE ROW ORDERING
C     IS RETURNED IN NROW
C
      DIMENSION A(N,N),NROW(N)
      DO 1 I=1,N
1     NROW(I) = I
      M = N-1
      DO 2 I=1,M
           IMAX  = I
           J = I+1
           DO 3 IP=J,N
           L1 = NROW(IMAX)
           L2 = NROW(IP)
3          IF(ABS(A(L2,I)).GT.ABS(A(L1,I))) IMAX = IP
           IF(ABS(A(NROW(IMAX),I)).LE.1.0E-30) THEN
                WRITE(6,100)
                STOP
           END IF
           JJ = NROW(I)
           NROW(I) = NROW(IMAX)
           NROW(IMAX) = JJ
           I1 = NROW(I)
           DO 4 JJ=J,N
                J1 = NROW(JJ)
                A(J1,I) = A(J1,I)/A(I1,I)
                DO 5 K=J,N
5               A(J1,K) = A(J1,K)-A(J1,I)*A(I1,K)
4          CONTINUE
2     CONTINUE
      RETURN
100   FORMAT(1X,'MULTIP FAILS')
      END
C
      SUBROUTINE SOLVE(N,A,NROW,X,Y)
C
C     SUBROUTINE SOLVE ACCEPTS THE VECTOR X, THE MULTIPLIERS AND UPPER
C     TRIANGULAR FORM FOR A, THE ROW ORDERING NROW, SOLVES A*Y = X, AND
C     RETURNS THE SOLUTION IN Y
C
      DIMENSION A(N,N),NROW(N),X(N),Y(N)
      M = N-1
      DO 2 I=1,M
           J=I+1
           I1 = NROW(I)
           DO 4 JJ=J,N
                J1 = NROW(JJ)
4          X(J1) = X(J1)-A(J1,I)*X(I1)
2     CONTINUE
      IF(ABS(A(NROW(N),N)).LE.1.0E-30) THEN
           WRITE(6,100)
           STOP
      END IF
      N1 = NROW(N)
      Y(N) = X(N1)/A(N1,N)
      L = N-1
      DO 15 K=1,L
           J = L-K+1
           JJ = J+1
           N2 = NROW(J)
           Y(J)=X(N2)
           DO 16 KK=JJ,N
16         Y(J) = Y(J)-A(N2,KK)*Y(KK)
15    Y(J) = Y(J)/A(N2,J)
      RETURN
100   FORMAT(1X,'SOLVE FAILS')
      END
