program ALG063;
{  GAUSSIAN ELIMINATION WITH SCALED COLUMN PIVOTING

   To solve the n by n linear system

   E1:  A[1,1] X[1] + A[1,2] X[2] +...+ A[1,n] X[n] = A[1,n+1]
   E2:  A[2,1] X[1] + A[2,2] X[2] +...+ A[2,n] X[n] = A[2,n+1]
   :
   .
   EN:  A[n,1] X[1] + A[n,2] X[2] +...+ A[n,n] X[n] = A[n,n+1]

   INPUT:   number of unknowns and equations n; augmented
            matrix A = (A(I,J)) where 1<=I<=n and 1<=J<=n+1.

   OUTPUT:  solution x(1), x(2),...,x(n) or a message that the
            linear system has no unique solution.

}
const
   ZERO = 1.0E-20;
var
   A : array [ 1..10, 1..11 ] of real;
   S : array [ 1..10 ] of real;
   NROW : array [ 1..10 ] of integer;
   AMAX,XM,SUM : real;
   FLAG,N,M,ICHG,I,NN,IMAX,J,JJ,IP,JP,NCOPY,I1,J1,N1,K,N2,LL,KK : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 14 ];
   INP,OUP : text;
procedure INPUT;
   begin
      writeln('This is Gauss Elimination with Scaled Column Pivoting.');
      writeln ('The array will be input from a text file in the order: ');
      writeln('A(1,1), A(1,2), ..., A(1,N+1), A(2,1), A(2,2), ..., A(2,N+1),');
      writeln ('..., A(N,1), A(N,2), ..., A(N,N+1) '); writeln;
      write ('Place as many entries as desired on each line, but separate ');
      writeln ('entries with ');
      writeln ('at least one blank. ');
      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 number of equations - an integer. ');
                  readln ( N );
                  if ( N > 0 ) then
                     begin
                        for I := 1 to N do
                           for J := 1 to N + 1 do read ( INP, A[I,J] );
                        OK := true;
                        close ( INP )
                     end
                  else writeln ('The number must be a positive integer. ')
               end
         end
      else
         begin
            writeln('The program will end so the input file can be created.');
            OK := false
         end
   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.OUT');
            readln ( NAME );
            assign ( OUP, NAME )
         end
      else assign ( OUP, 'CON' );
      rewrite(OUP);
      writeln(OUP,'GAUSSIAN ELIMINATION WITH SCALED COLUMN PIVOTING');
      writeln(OUP);
      writeln ( OUP, 'The reduced system - output by rows: ');
      for I := 1 to N do
         begin
            for J := 1 to N do write ( OUP, '':2, A[I,J]:12:8 );
            writeln ( OUP )
         end;
      writeln ( OUP ); writeln ( OUP );
      writeln ( OUP, 'Has solution vector: ');
      for I := 1 to N do
         begin
            J := NROW[I];
            write ( OUP, A[J,M]:12:8 )
         end;
      writeln ( OUP );
      writeln ( OUP, 'with ',ICHG:3,' row interchange(s) ');
      writeln ( OUP );
      writeln ( OUP, 'The rows have been logically re-ordered to: ');
      for I := 1 to N do write ( OUP, NROW[I]:3 ); writeln(OUP);
      close ( OUP )
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
            M := N + 1;
{           STEP 1                                                     }
            for I := 1 to N do
               begin
                  S[I] := abs( A[I,1] );
{                 initialize row pointer                               }
                  NROW[I] := I;
                  for J := 1 to N do
                     if ( abs( A[I,J] ) > S[I] ) then
                        S[I] := abs( A[I,J] );
                  if ( S[I] <= ZERO ) then OK := false
               end;
            NN := N - 1;
            ICHG := 0;
            I := 1;
{           STEP 2                                                     }
{           elimination process                                        }
            while ( OK ) and ( I <= NN ) do
               begin
{                 STEP 3                                               }
                  IMAX := NROW[I];
                  AMAX := abs( A[IMAX,I] ) / S[IMAX];
                  IMAX := I;
                  JJ := I + 1;
                  for IP := JJ to N do
                     begin
                        JP := NROW[IP];
                        if ( abs( A[JP,I] ) / S[JP] > AMAX ) then
                           begin
                              AMAX := abs( A[JP,I] ) / S[JP];
                              IMAX := IP
                           end
                     end;
{                    STEP 4                                            }
{                    system has no unique solution                     }
                     if ( AMAX <= ZERO ) then OK := false
                     else
                        begin
{                          STEP 5                                      }
{                          simulate row interchange                    }
                           if ( NROW[I] <> NROW[IMAX] ) then
                              begin
                                 ICHG := ICHG + 1;
                                 NCOPY := NROW[I];
                                 NROW[I] := NROW[IMAX];
                                 NROW[IMAX] := NCOPY
                              end;
{                          STEP 6                                      }
                           I1 := NROW[I];
                           for J := JJ to N do
                              begin
                                 J1 := NROW[J];
{                                STEP 7                                }
                                 XM := A[J1,I] / A[I1,I];
{                                STEP 8                                }
                                 for K := JJ to M do
                                    A[J1,K] := A[J1,K] - XM * A[I1,K];
                                 A[J1,I] := XM
                              end
                        end;
                     I := I + 1
                  end;
               if ( OK ) then
                  begin
{                    STEP 9                                            }
                     N1 := NROW[N];
                     if ( abs( A[N1,N] ) <= ZERO ) then OK := false
{                    system has no unique solution                     }
                     else
                        begin
{                          STEP 10                                     }
{                          start backward substitution                 }
                           A[N1,M] := A[N1,M] / A[N1,N];
{                          STEP 11                                     }
                           for K := 1 to NN do
                              begin
                                 I := NN - K + 1;
                                 JJ := I + 1;
                                 N2 := NROW[I];
                                 SUM := A[N2,N+1];
                                 for KK := JJ to N do
                                    begin
                                       LL := NROW[KK];
                                       SUM := SUM - A[N2,KK] * A[LL,M]
                                    end;
                                 A[N2,M] := SUM / A[N2,I]
                              end;
{                          STEP 12                                     }
{                          procedure completed successfully            }
                           OUTPUT
                        end
                  end;
               if ( not OK ) then writeln ('System has no unique solution ')
         end
   end.