Factorización LU

Método para la obtención de la factorización de una matriz A de modo que A=LU tal que L sea una matriz superior y U una matriz inferior. La utilidad en este tipo de factorizaciones reside en la facilidad que supone resolver un sistema triangular, además de tener un coste menor que métodos como el de Gauss.

program factorizacionLU(input,output);
const
  tanMax=20;
type
  tpMatriz=array[1..tanMax,1..tanMax] of real;
  tpDimension=1..maxint;
var
  A,L,U:tpMatriz;
  tan:tpDimension;

procedure Inicializa(VAR x:tpMatriz);
  var
    i,j:tpDimension;
  begin
    for i:=1 to tan do
      for j:=1 to tan do
        x[i,j]:=0
  end;

procedure leedatos(VAR x:tpMatriz; VAR tan:tpDimension);
  var
    i,j:tpDimension;
  begin
    write('Tamaño de la matriz: ');readln(tan);
    writeln;
    for i:=1 to tan do begin
      for j:=1 to tan do begin
        write('Termino(',i,',',j,'): '); readln(x[i,j])
      end;
      writeln
    end
  end;

procedure multiplica(VAR x:tpMatriz; y,z:tpMatriz);
  {multiplica Y*Z en este orden y lo guarda en X}
  var
    i,j,k:tpDimension;
  begin
    for i:=1 to tan do
      for j:=1 to tan do begin
        x[i,j]:=0;
        for k:=1 to tan do x[i,j]:=x[i,j]+y[i,k]*z[k,j]
      end
  end;

procedure resta(VAR x:tpMatriz; y,z:tpMatriz);
  var
    i,j:tpDimension;
  begin
    for i:=1 to tan do
      for j:=1 to tan do
        x[i,j]:=y[i,j]-z[i,j]
  end;

procedure factoriza(VAR x,y:tpMatriz; A:tpMatriz);
  var
    i,paso:tpDimension;
    aux,Vu,Vl:tpMatriz;
  begin
    Inicializa(Vu);Inicializa(Vl);Inicializa(x);Inicializa(y);
    for paso:=1 to tan do begin
      for i:=1 to tan do begin
        y[paso,i]:=A[paso,i];
        Vu[1,i]:=y[paso,i];
        x[i,paso]:=A[i,paso]/A[paso,paso];
        Vl[i,1]:=x[i,paso]
      end;
      multiplica(aux,Vl,Vu);
      resta(A,A,aux)
    end
  end;

procedure escribematriz(m:tpMatriz; tan:tpDimension);
  var
    i,j:tpDimension;
  begin
    writeln;
    for i:=1 to tan do begin
      for j:=1 to tan do
        write(m[i,j]:7:3);
      writeln
    end;
    writeln;writeln
  end;

{comienzo programa}
begin
  Inicializa(A);
  leedatos(A,tan);
  writeln('La matriz introducida es: ');
  escribematriz(A,tan);
  factoriza(L,U,A);
  writeln('La matriz U es:');
  escribematriz(U,tan);
  writeln('La matriz L es:');
  escribematriz(L,tan);
  readln
end.



Salud.