INVERSA DE UNA MATRIZ (mediante operaciones elementales)

El siguiente código está escrito en Pascal. En este otro enlace se lo dejo en C.

Codigo en C.

El programa se basa en el programa de Resolución de sistemas mediante Gauss-Jordan, lo único que las operaciones elementales que se le van realizando a la matriz original también se le realizan a la matriz identidad del mismo rango.

program CalcInversa (input,output);
uses crt;


const
  MaxDim=1000;
  error=0.00000001;{Valor por debajo del cual el programa considerara 0}
type
  tpContador=1..Maxint;
  tpMatriz=record
    nums:array[1..MaxDim,1..MaxDim] of real;
    dim:tpContador
  end;
  tpNom=string[50];
var
  matriz:tpMatriz;
  opcion:tpContador;
  determinado:boolean;


procedure EscribeTitulo;
  begin
    writeln('--CALCULO DE LA INVERSA--':40);
    writeln('Calcula mediante el producto con matrices elementales la matriz inversa..');
    writeln('Si no existe la inversa de la matriz introducida el programa se lo indicara.');
    writeln('De ahora en adelante, a la hora de seleccionar una opcion: 0=SI   y   Cualquier otra tecla=NO')
  end;


procedure EscribirMatriz(matriz:tpMatriz);
  var
    c1,c2:tpContador;
  begin
    for c1:=1 to matriz.dim do begin
        for c2:=1 to (matriz.dim+1) do
          write(matriz.nums[c1,c2]:10:2);
        writeln
    end;
    writeln
  end;


procedure EscribirInversa(matriz:tpMatriz);
  var
    c1,c2:tpContador;
  begin
    for c1:=1 to matriz.dim do begin
        for c2:=(matriz.dim+1) to 2*matriz.dim do
          write(matriz.nums[c1,c2]:10:2);
        writeln
    end;
    writeln
  end;


procedure PideDatos(VAR matriz:tpMatriz);
  var
    c1,c2:tpContador;
  begin
    write('Introduzca la dimension de la matriz: ');readln(matriz.dim);clrscr;EscribeTitulo;
      for c1:=1 to matriz.dim do
        for c2:=1 to matriz.dim do begin
          if c1=c2 then
            matriz.nums[c1,c2+matriz.dim]:=1
          else
            matriz.nums[c1,c2+matriz.dim]:=0;
          write('Introduzca el termino (',c1,',',c2,'): ');readln(matriz.nums[c1,c2]);clrscr;
          EscribeTitulo
        end;
    EscribeTitulo;writeln;
    writeln;writeln('La matriz introducida es la siguiente:');
    EscribirMatriz(matriz);
  end;


procedure Invierte(VAR matriz:tpMatriz;VAR determinado:boolean);
  var
    paso,c1,c2:tpContador;
    PivCorrect:boolean;
    pivote,aux:real;
  begin{0}
    for paso:=1 to matriz.dim do begin{1}
      PivCorrect:=false;
      c1:=paso;
      while (not PivCorrect) and (c1<=matriz.dim) do
        If abs(matriz.nums[c1,paso])>error then
          PivCorrect:=true
        else
          c1:=c1+1;
      If PivCorrect then begin{3}
        pivote:=matriz.nums[c1,paso];
        for c2:=paso to (matriz.dim*2) do
          if c1<>paso then begin
            aux:=matriz.nums[paso,c2];
            matriz.nums[paso,c2]:=matriz.nums[c1,c2]/pivote;
            matriz.nums[c1,c2]:=aux
          end else
            matriz.nums[c1,c2]:=matriz.nums[c1,c2]/pivote;
        {Hasta aquí ha sido solo preparar el pivote para hacer ceros por debajo
        el pivote en estos momentos es 1}
      end;{3}
     for c1:=(paso+1) to matriz.dim do begin
       aux:=matriz.nums[c1,paso];
       for c2:=paso to (matriz.dim*2) do
         matriz.nums[c1,c2]:=matriz.nums[c1,c2]-aux*matriz.nums[paso,c2]
     end;
    end;{1}
    {Aqui la matriz ya esta escalonada (se imprime en pantalla). Se comprueba que el sistema sea determinado}
    determinado:=true;
    for c1:=1 to matriz.dim do                   
      if abs(matriz.nums[c1,c1])<error then
        determinado:=false;


    if determinado then begin
      for paso:=matriz.dim downto 1 do begin
        pivote:=matriz.nums[paso,paso];
        matriz.nums[paso,paso]:=1;
        for c2:=matriz.dim+1 to 2*matriz.dim do
          matriz.nums[paso,c2]:=matriz.nums[paso,c2]/pivote;
        for c1:=(paso-1) downto 1 do begin
          aux:=matriz.nums[c1,paso];
          for c2:=paso to 2*matriz.dim do
            matriz.nums[c1,c2]:=matriz.nums[c1,c2]-matriz.nums[paso,c2]*aux
        end
      end;
      writeln('La matriz invertida es: ');EscribirInversa(matriz); writeln {Aqui la matriz ya esta diagonalizada}
    end
    else
      writeln('La matriz no tiene inversa')
  end;{0}


procedure GuardaFich(matriz:tpMatriz);
  var
    c1,c2:tpContador;
    fichero:text;
    nomfich:tpNom;
  begin
    write('Introduzca el nombre del fichero (sin la extension): ');readln(nomfich);
    nomfich:=nomfich+'.txt';
    assign(fichero,nomfich);
    rewrite(fichero);
    writeln(fichero,'La matriz inversa es:');
    for c1:=1 to matriz.dim do begin
      for c2:=matriz.dim+1 to matriz.dim*2 do
        write(fichero,matriz.nums[c1,c2]:10:3);
      writeln(fichero)
    end;
    close(fichero)
  end;


begin
  opcion:=1; writeln;
  while (opcion=1) do begin
    EscribeTitulo;
    PideDatos(matriz);writeln;
    Invierte(matriz,determinado); writeln;
    If determinado then begin
      write('Quiere guardar la matriz inversa en un fichero? ');readln(opcion);
      If opcion=1 then
        GuardaFich(matriz)
    end;
    write('Quiere volver a trabajar con el programa? ');readln(opcion);
    clrscr {limpia pantalla}
  end
end.

Salud!