program ALG093;
{  INVERSE POWER METHOD

   To approximate an eigenvalue and an associated eigenvector of the
   n by n matrix A given a nonzero vector x:

   INPUT:   Dimension n; matrix A; vector x; tolerance TOL;
            maximum number of iterations N.

   OUTPUT:  Approximate eigenvalue MU; approximate eigenvector x
            or a message that the maximum number of iterations was
            exceeded.
}
const
   ZERO = 1.0E-20;
var
   A : array [ 1..10, 1..10 ] of real;
   X,Y,B : array [ 1..10 ] of real;
   NROW : array [1..10] of integer;
   T,AMAX,YMU,ERR,TOL,Q,S : real;
   FLAG,N,I,J,NN,K,LP : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 14 ];
   INP, OUP : text;
procedure INPUT;
   begin
      writeln('This is the Inverse Power Method.');
      OK := false;
      writeln ('The array will be input from a text file in the order: ');
      writeln('A(1,1), A(1,2), ..., A(1,n), A(2,1), A(2,2), ..., A(2,n),');
      write ('..., A(n,1), ');
      writeln ('A(n,2), ..., A(n,n) '); writeln; writeln;
      write ('Place as many entries as desired on each line, but separate ');
      writeln ('entries with ');
      writeln ('at least one blank. ');
      writeln ('The initial approximation should follow in same format.' );
      writeln; writeln;
      writeln ('Has the input file been created? - enter Y or N. ');
      readln ( AA );
      if ( AA = 'Y' ) or ( AA = 'y' ) then
         begin
            writeln ('Input the file name in the form - drive:name.ext, ');
            writeln ('for example: A:DATA.DTA ');
            readln ( NAME );
            assign ( INP, NAME );
            reset ( INP );
            OK := false;
            while ( not OK ) do
               begin
                  writeln ('Input the dimension n. ');
                  readln ( N );
                  if ( N > 0 ) then
                     begin
                        for I := 1 to N do
                           for J := 1 to N do read ( INP, A[I,J] );
                        for I := 1 to N do read ( INP, X[I]);
{                       use X1 for XO                                      }
                        close ( INP );
                        while ( not OK ) do
                           begin
                              writeln ('Input the tolerance. ');
                              readln ( TOL );
                              if ( TOL > 0.0 ) then OK := true
                              else
                                 writeln ('Tolerance must be positive number.')
                           end;
                        OK := false;
                        while ( not OK ) do
                           begin
                              write ('Input maximum number of iterations ');
                              writeln ('- integer.');
                              readln ( NN );
{                             use NN for N                             }
                              if ( NN > 0 ) then OK := true
                              else
                                 writeln ('Number must be positive integer. ')
                           end
                     end
                  else writeln ('The dimension must be a positive integer. ')
               end
         end
      else writeln ('The program will end so the input file can be created. ')
   end;
   procedure OUTPUT;
      begin
         writeln ('Choice of output method: ');
         writeln ('1. Output to screen ');
         writeln ('2. Output to text file ');
         writeln ('Please enter 1 or 2. ');
         readln ( FLAG );
         if ( FLAG = 2 ) then
            begin
               writeln ('Input the file name in the form - drive:name.ext, ');
               writeln('for example:   A:OUTPUT.DTA');
               readln ( NAME );
               assign ( OUP, NAME )
            end
         else  assign ( OUP, 'CON' );
         rewrite ( OUP );
         writeln('INVERSE POWER METHOD');
         writeln ( OUP );
      end;
procedure MULTIP;
   var
      K,I,M,IMAX,J,IP,L1,L2,JJ,I1,J1,N1 : integer;
   begin
      for I := 1 to N do NROW[I] := I;
      OK := true;
      I := 1;
      M := N - 1;
      while ( ( I <= M ) and OK ) do
         begin
            IMAX := I;
            J := I + 1;
            for IP := J to N do
               begin
                  L1 := NROW[IMAX];
                  L2 := NROW[IP];
                  if ( abs( A[L2,I] ) > abs( A[L1,I] ) ) then IMAX := IP
               end;
            if ( abs( A[NROW[IMAX],I] ) <= ZERO ) then
               begin
                  OK := false;
                  write ('A - Q * I is singular, Q = ', Q, ' is an ');
                  writeln ('eigenvalue ')
               end
            else
               begin
                  JJ := NROW[I];
                  NROW[I] := NROW[IMAX];
                  NROW[IMAX] := JJ;
                  I1 := NROW[I];
                  for JJ := J to N do
                     begin
                        J1 := NROW[JJ];
                        A[J1,I] := A[J1,I] / A[I1,I];
                        for K := J to N do
                           A[J1,K] := A[J1,K] - A[J1,I] * A[I1,K]
                     end
               end;
            I := I + 1
         end;
      if ( abs( a[NROW[N],N] ) <= ZERO ) then
         begin
            OK := false;
            writeln ('A - Q * I is singular, Q = ', Q, ' is an eigenvalue ')
         end
   end;
procedure SOLVE;
   var
      M,I,J,I1,JJ,J1,N1,N2,L,K,KK : integer;
   begin
      M := N - 1;
      for I := 1 to M do
         begin
            J := I + 1;
            I1 := NROW[I];
            for JJ := J to N do
               begin
                  J1 := NROW[JJ];
                  B[J1] := B[J1] - A[J1,I] * B[I1]
               end;
         end;
      N1 := NROW[N];
      Y[N] := B[N1] / A[N1,N];
      L := N - 1;
      for K := 1 to L do
         begin
            J := L - K + 1;
            JJ := J + 1;
            N2 := NROW[J];
            Y[J] := B[N2];
            for KK := JJ to N do Y[J] := Y[J] - A[N2,KK] * Y[KK];
            Y[J] := Y[J] / A[N2,J]
         end
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
            OUTPUT;
{           STEP 1                                                             }
{           Q could be input instead of computed by deleting
            the next 7 steps                                           }
            Q := 0.0;
            S := 0.0;
            for I := 1 to N do
               begin
                  S := S + X[I] * X[I];
                  for J := 1 to N do Q := Q + A[I,J] * X[I] * X[J]
               end;
            Q := Q / S;
            writeln('Q is ',Q);
            writeln('input new Q if desired');
            readln(Q);
            writeln(OUP,'Iteration  Eigenvalue  Eigenvector');
{           STEP 2                                                     }
            K := 1;
            for I := 1 to N do A[I,I] := A[I,I] - Q;
{           Call subroutine to compute multipliers M(I,J) and
            upper triangular matrix for matrix A using Gauss
            elimination with maximal column pivoting.
            NROW holds the ordering of rows for interchanges           }
            Multip;
            if ( OK ) then
               begin
{                 STEP 3                                               }
                  LP := 1;
                  for I := 2 to N do
                     if ( abs( X[I] ) > abs( X[LP] ) ) then LP := I;
{                 STEP 4                                               }
                  AMAX := X[LP];
                  for I := 1 to N do X[I] := X[I] / (AMAX);
{                 STEP 5                                               }
                  while ( ( K <= NN ) and OK ) do
                     begin
{                       STEP 6                                         }
                        for I := 1 to N do B[I] := X[I];
{                       Subroutine solve returns the solution of
                        ( A - Q * I )Y = B in Y                        }
                        SOLVE;
{                       STEP 7                                         }
                        YMU := Y[LP];
{                       STEP 8 AND 9                                   }
                        LP := 1;
                        for I := 2 to N do
                           if ( abs( Y[I] ) > abs( Y[LP] ) ) then LP := I;
                        AMAX := Y[LP];
                        ERR := 0.0;
                        for I := 1 to N do
                           begin
                              X[I] := X[I];
                              T := Y[I] / AMAX;
                              if ( abs( X[I] - T ) > ERR ) then
                                 ERR := abs( X[I] - T );
                              X[I] := T
                           end;
                        YMU := 1 / YMU + Q;
{                       STEP 10                                        }
                        writeln(OUP,K:3,' ',YMU:12:8);
                        for I := 1 to N do write(OUP,X[I]:12:8);
                        writeln(OUP);
                        if ( ERR < TOL ) then
                           begin
                              OK := false;
                              write(OUP,'Eigenvalue = ',YMU:12:8);
                              writeln(OUP,' to tolerance = ',TOL);
                              writeln(OUP,'obtained on iteration number ',K);
                              writeln ( OUP );
                              writeln (OUP,'Unit eigenvector is : ');
                              for I := 1 to N do write ( OUP, X[I]:12:8 );
                           end
                        else
{                          STEP 11                                     }
                           K := K + 1
                     end;
                  if K > NN then writeln('Method fails');
               end;
            close ( OUP )
         end
{  STEP 12                                                             }
   end.