program ALG092;
{ SYMMETRIC POWER METHOD ALGORITHM 9.2

   To approximate the dominant eigenvalue and an associated
   eigenvector of the n by n symmetric 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 : array [ 1..10 ] of real;
   YMU,ERR,TOL : real;
   FLAG,N,I,J,NN,K : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 30 ];
   INP,OUP : text;
procedure INPUT;
begin
   writeln('This is the Symmetric 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;
   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] );
{                    Because of the way in which the subroutine
                     norm computes the l2 norm of a vector,
                     the initial input of X is into Y and X
                     is initialized at the zero vector.             }
                     for I := 1 to N do read ( INP, Y[I] );
                     for I := 1 to N do X[I] := 0.0;
                     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 capital N for the maximum number of
                           iterations                                     }
                           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(OUP,'SYMMETRIC POWER METHOD');
      writeln(OUP);
      writeln(OUP,'iter   approx        approx eigenvector');
      writeln(OUP,'       eigenvalue');
   end;
procedure NORM;
   var
      XL,T : real;
      I : integer;
   begin
      XL := 0.0;
      for I := 1 to N do XL := XL + sqr( Y[I] );
{     2-norm of Y    }
      XL := sqrt( XL );
      ERR := 0.0;
      if ( XL > ZERO ) then
         begin
            for I := 1 to N do
               begin
                  T := Y[I] / XL;
                  ERR := ERR + sqr(X[I]-T);
                  X[I] := T
               end;
{              X has a 2-norm of 1.0             }
            ERR := sqrt(ERR)
         end
      else
         begin
            write ('A has a zero eigenvalue - select new vector and begin ');
            writeln ('again ');
            OK := false
         end
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
            OUTPUT;
{           STEP 1                                                     }
            K := 1;
{           Norm computes the norm of the vector X and
            returns X divided by its norm in the variable Y            }
            NORM;
            if ( OK ) then
               begin
{                 STEP 2                                               }
                  while ( ( K <= NN ) and OK ) do
                     begin
{                       STEPS 3 AND 4                                  }
                        YMU := 0.0;
                        for I := 1 to N do
                           begin
                              Y[I] := 0.0;
                              for J := 1 to N do
                                 Y[I] := Y[I] + A[I,J] * X[J];
                              YMU := YMU + X[I] * Y[I]
                           end;
{                       STEP 5 (This step is accomplished in
                                subroutine norm.)                      }
{                       STEP 6                                         }
                        NORM;
                        write ( OUP,K,' ', YMU:12:8,'   ');
                        for I := 1 to N do write ( OUP, X[I]:12:8 );
                        writeln(OUP);
                        if ( OK ) then
                           begin
{                             STEP 7                                   }
                              if ( ERR < TOL ) then
                                 begin
{                                   procedure completed successfully   }
                                    writeln ( OUP ); writeln ( OUP );
                                    write(OUP,'The eigenvalue = ',YMU:12:8);
                                    writeln(OUP,' to tolerance = ',TOL);
                                    writeln(OUP,'on iteration number = ',K);
                                    writeln ( OUP );
                                    writeln(OUP,'Unit eigenvector is:');
                                    writeln ( OUP ) ;
                                    for I := 1 to N do write(OUP,X[I]:12:8);
                                    writeln(OUP);
                                    OK := false
                                 end
                              else
{                                STEP 8                                }
                                 K := K + 1
                           end
                     end;
{                 STEP 9                                               }
                  if ( K > NN ) then
                     writeln ('No convergence within ', NN, ' iterations ')
               end;
               close(OUP)
         end;
   end.