program ALG066;
{  CHOLESKI'S ALGORITHM 6.6

   To factor the positive definite n by n matrix A into LL**T,
   where L is lower triangular.

   INPUT:   the dimension n; entries A(I,J), 1<=I, J<=n of A.

   OUTPUT:  the entries L(I,J), 1<=J<=I, 1<=I<=n of L.

   the entries of U = L**T are U(I,J) = L(J,I), I<=J<=n, 1<=I<=n
}
var
   A : array [ 1..10, 1..10 ] of real;
   S : real;
   FLAG,N,I,J,K,NN,JJ,KK : integer;
   OK : boolean;
   AA : char;
   NAME : string [ 30 ];
   INP,OUP : text;
procedure INPUT;
   begin
      writeln('This is Choleski Factorization Method.');
      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),');
      writeln ('..., A(N,1), 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; writeln;
      writeln ('Has the input file been created? - enter Y or N. ');
      readln ( AA );
      OK := false;
      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 - an integer. ');
                  readln ( N );
                  if ( N > 0 ) then
                     begin
                        for I := 1 to N do
                           for J := 1 to N do read ( INP, A[I,J] );
                        OK := true;
                        close ( INP )
                     end
                  else writeln ('The number 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,'CHOLESKI FACTORIZATION');
      writeln(OUP);
      writeln ( OUP, 'The matrix l output by rows: ');
      for I := 1 to N do
         begin
            for J := 1 to I do write ( OUP, '':2, A[I,J]:12:8 );
            writeln ( OUP )
         end;
      close ( OUP )
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
{           STEP 1                                                     }
            A[1,1] := sqrt( A[1,1] );
{           STEP 2                                                     }
            for J := 2 to N do A[J,1] := A[J,1] / A[1,1];
{           STEP 3                                                     }
            NN := N - 1;
            for I := 2 to NN do
               begin
{                 STEP 4                                               }
                  KK := I - 1;
                  S := 0.0;
                  for K := 1 to KK do S := S - A[I,K] * A[I,K];
                  A[I,I] := sqrt( A[I,I] + S );
{                 STEP 5                                               }
                  JJ := I + 1;
                  for J := JJ to N do
                     begin
                        S := 0.0;
                        KK := I - 1;
                        for K := 1 to KK do S := S - A[J,K] * A[I,K];
                        A[J,I] := ( A[J,I] + S ) / A[I,I]
                     end
               end;
{           STEP 6                                                     }
            S := 0.0;
            for K := 1 to NN do S := S - A[N,K] * A[N,K];
            A[N,N] := sqrt( A[N,N] + S );
{           STEP 7                                                     }
            OUTPUT
         end
   end.