{$N+}
(* Program nacte matici A ze souboru matice.txt a rozlozi ji pomoci LU dekompozice *)
(*vysledek vypise na obrazovku *)
program LUDekompozice;
const
  RozmerMatice=5;
type
  TypMatice = double;
  Matice=array[1..RozmerMatice,1..RozmerMatice] of TypMatice;

var
  A,L,U:Matice;
  soucin:Matice;

procedure nactiMatici(jmenoSouboru:string; var m:Matice);
var
  soubor:Text;
  i,j:integer;
begin
  assign(soubor,jmenoSouboru);
  reset(soubor);

  for i := 1 to RozmerMatice do
  begin
    for j:= 1 to RozmerMatice do
    begin
      read(soubor,m[i,j]);
    end;
    readln(soubor);
  end;

  close(soubor);
end;

procedure zobrazMatici(m:Matice);
var
  i,j:integer;
begin
  for i := 1 to RozmerMatice do
  begin
    for j:= 1 to RozmerMatice do
    begin
      write(m[i,j]:8:3);
      write(' ');
    end;
    writeln;
  end;
end;

procedure nasobMatice(m1,m2:Matice; var vysledek:Matice);
var
	i,j,k:integer;
begin
  for i := 1 to RozmerMatice do
  begin
    for j := 1 to RozmerMatice do
    begin
      vysledek[i,j] := 0;
      for k := 1 to RozmerMatice do vysledek[i,j] := vysledek [i,j] + m1[i,k]*m2[k,j];
    end;
  end;
end;

procedure spoctiLUDekompozici(A:Matice; var L,U:Matice);
var
  radek,sloupec,k,diagonala:integer;
  soucet:TypMatice;
begin
  (* napred matice vynulujeme *)
  for radek := 1 to RozmerMatice do
  begin
    for sloupec := 1 to RozmerMatice do begin
        L[radek,sloupec] := 0.0;
        U[radek,sloupec] := 0.0;
    end;
    (* a do diagonaly  nastavime na 1 *)
    L[radek,radek] := 1.0;
  end;

  for diagonala := 1 to RozmerMatice do (* pro prvek na diagonale *)
  begin
    (* projedeme radek do konce *)
    for sloupec := diagonala to RozmerMatice do
    begin
      soucet := 0.0;
      for k := 1 to diagonala-1 do
      begin
        soucet := soucet + L[diagonala,k] * U[k,sloupec];
      end;
      U[diagonala,sloupec] := A[diagonala,sloupec] - soucet;
    end;
    (* a nyni projedeme sloupecek do konce, bez diagonaly *)
    for radek := diagonala + 1 to RozmerMatice do
    begin
      soucet := 0.0;
      for k := 1 to diagonala-1 do
      begin
        soucet := soucet + L[radek,k] * U[k,diagonala];
      end;
      L[radek,diagonala] := (A[radek,diagonala]-soucet)/U[diagonala,diagonala];
    end;
  end;

end;

begin
  nactiMatici('matice.txt',A);

  (* rozlozime *)
  spoctiLUDekompozici(A,L,U);
  (* overime zpetne roznasobeni *)
  nasobMatice(L,U,soucin);
  (* vypiseme puvodni matici *)
  writeln('A:');
  zobrazMatici(A);
  readln;
  (* vypiseme matici L *)
  writeln('L:');
  zobrazMatici(L);
  readln;
  (* vypiseme matici U *)
  writeln('U:');
  zobrazMatici(U);
  readln;
  (* vypiseme kontrolni soucin LU *)
  writeln('Soucin:');
  zobrazMatici(soucin);
  readln;

end.