Resolucion de sistemas mediante JACOBI

El siguiente programa muestra uno de los metodos iterativos más conocidos (posiblemente el que más, junto con relajación y Gauss-Seidel), denominado JACOBI. Este programa calcula la solución de sistemas de ecuaciones lineales, siempre y cuando se cumpla el criterio de convergencia del método, que en este caso es que ρ(B_J), que corresponde al máximo de los valores absolutos de las raíces de la ecuación característica de la matriz B_J (det(B_J - λI)) sea menor que 1. (Una condición suficiente pero no necesaria es que la matriz sea diagonalmente domiante, que es más facil de ver). El programa pide un numero de iteraciones a llevar a cabo, de modo que cuanto mayor sea esta cantidad, mayor será la precisión de la solución, siempre y cuando converja. Por lo demás, no me he molestado en hacer que la solución se guarde en un fichero de texto, pero vamos, poco tiempo llevaría. Aquí les dejo el programa:

program JACOBI(input,output);
const
  maxdim=300;
type
  Matriz=record
    comp:array[1..maxdim,1..maxdim] of real;
    dim:integer
  end;
  Vector=record
    comp:array[1..maxdim] of real;
    dim:integer
  end;
  tpContador=1..MAXINT;
var
  c:tpContador;
  Solucion,K,B:Vector;
  M,J,Dinv,L,U:Matriz;{Matriz de iteración del método de Jacobi}
  iteraciones:tpContador;

procedure CargaDatos(VAR M:Matriz; VAR B:Vector; VAR Iteraciones:tpContador);
  var
    i,h:tpContador;
  begin
    writeln;
    writeln('VISITANOS EN:   http://www.matematicapascal.blogspot.com/');writeln;
    writeln('Resuelve sistemas mediante JACOBI');
    writeln;writeln;
    write('Introduzca la dimension de la matriz del sistema: ');readln(M.dim);
    B.dim:=M.dim;
    for i:=1 to M.dim do begin
      writeln;
      for h:=1 to M.dim do begin
        write('    Introduzca M(',i,',',h,'): '); readln(M.comp[i,h])
      end
    end;
    writeln; writeln;
    writeln('Ahora introduzca los terminos independientes:');
    for i:=1 to B.dim do begin
      write('Introduzca B(',i,'): ');readln(B.comp[i])
    end;
    write('Por ultimo introduzca el numero de iteraciones a llevar a cabo: ');readln(iteraciones)
  end;

procedure EscribeMatriz(M:Matriz; B,Solucion:Vector);
{Guarda los datos en un fichero}
  var
    i,j:tpContador;
  begin
    writeln('El sistema introducido es: ');
    writeln;
    for i:=1 to M.dim do begin
      for j:=1 to M.dim do
        write(M.comp[i,j]:7:3);
          write('  |  ',B.comp[i]:7:3); writeln
    end;
    writeln;writeln;
    writeln('La solucion del sistema es: ');
    writeln;
    for i:=1 to Solucion.dim do
      writeln(Solucion.comp[i]:7:3);
    writeln
  end;

procedure MultiplicaMats(VAR x:Matriz; y,z:Matriz);
  var
    i,h,k:tpContador;
  begin
    x.dim:=y.dim;
    for i:=1 to y.dim do
      for h:=1 to y.dim do begin
        x.comp[i,h]:=0;
        for k:=1 to y.dim do x.comp[i,h]:=x.comp[i,h]+y.comp[i,k]*z.comp[k,h]
      end
  end;

procedure MultiplicaMatVect(VAR x:Vector; y:Matriz; z:Vector);
  var
    i,j:tpContador;
  begin
    x.dim:=y.dim;
    for i:=1 to y.dim do begin
      x.comp[i]:=0;
      for j:=1 to y.dim do
        x.comp[i]:=x.comp[i]+y.comp[i,j]*z.comp[j]
    end
  end;

procedure SumaMat(VAR x:Matriz; y,z:Matriz);
  var
    i,h:tpContador;
  begin
    x.dim:=y.dim;
    for i:=1 to y.dim do
      for h:=1 to y.dim do
        x.comp[i,h]:=y.comp[i,h]+z.comp[i,h]
  end;

procedure RestaVect(VAR x:Vector; y,z:Vector);
  {Calcula Y-Z y lo guarda en X}
  var
    i:tpContador;
  begin
    x.dim:=y.dim;
    for i:=1 to x.dim do
      x.comp[i]:=y.comp[i]-z.comp[i]
  end;

procedure SeleccionaDinvLyU(M:Matriz; VAR L,U,Dinv:Matriz; VAR B:Vector; VAR J:Matriz; VAR K:Vector);
  var
    Maux:Matriz;
    i,h:tpContador;
  begin
    Dinv.dim:=M.dim;
    L.dim:=M.dim;
    U.dim:=M.dim;
    for i:=1 to M.dim do
      for h:=1 to M.dim do
        if i=h then begin
          Dinv.comp[i,h]:=(1/M.comp[i,h]);
          U.comp[i,h]:=0;
          L.comp[i,h]:=0
        end
        else if i>h then
        begin
          Dinv.comp[i,h]:=0;
          U.comp[i,h]:=0;
          L.comp[i,h]:=M.comp[i,h]
        end
        else begin
          Dinv.comp[i,h]:=0;
          L.comp[i,h]:=0;
          U.comp[i,h]:=M.comp[i,h]
        end;
    SumaMat(Maux,L,U);
    MultiplicaMats(J,Dinv,Maux);
    MultiplicaMatVect(K,Dinv,B);
  end;
{Ya estan seleccionadas las matrices Dinv, L y U}

procedure IniciaVectorSolucion(M:Matriz; VAR Solucion:Vector);
  var
    i:tpContador;
  begin
    Solucion.dim:=M.dim;
    for i:=1 to Solucion.dim do
      Solucion.comp[i]:=0
  end;

procedure Paso(J:Matriz; VAR K,Solucion:Vector);
  var
    Vaux:Vector;
  begin
    MultiplicaMatVect(Vaux,J,Solucion);
    RestaVect(Solucion,K,Vaux);
  end;

{comienzo del programa}
begin
  CargaDatos(M,B,Iteraciones);{hasta aqui bien}
  IniciaVectorSolucion(M,Solucion);
  SeleccionaDinvLyU(M,L,U,Dinv,B,J,K);
  for c:=1 to Iteraciones do
    Paso(J,K,Solucion);
  EscribeMatriz(M,B,Solucion);
  readln
end.

Salud!!