Der Gauss-Jordan-Algorithmus bestimmt eine LR-Zerlegung
 
   P A = L R
 
von A. Dabei ist P eine Permutationsmatrix, die durch Zeilen- bzw. Spaltenvertauschungen aus einer Einheitsmatrix hervorgeht, sowie L eine untere und R eine obere Dreiecksmatrix.

Durch Vorwärts- und Rückwärtseinsetzen erhält man aus P, L und R die zu A inverse Matrix.

Als Pseudo Code sieht das Verfahren wie folgt aus:

procedure gaussInverse(n:vectorIndex; (* Größe der Matrix n x n *)
                       const A:matrix; (* Matrix A (Eingabe) *)
                       var Ainv:matrix; (* invertiertes A (Ausgabe) *)
                       var losbar:boolean); (* lösbar ? *)
const epsilon=2.2e-160;
var lu:matrix;
    d:vector; (* Zeilenmaxima *)
    permutation:array[vectorIndex] of vectorIndex;
    i,k,j,j0,m:vectorIndex;
    zeilenmaximum,abslukj,lukk,piv,temp,summe:extended;
    
  procedure tausch(a,b:vectorIndex);
  var i:vectorIndex;
      r:extended;
      z:vector;
  begin
    i:=permutation[a];
    permutation[a]:=permutation[b];
    permutation[b]:=i;
    r:=d[a];
    d[a]:=d[b];
    d[b]:=r;
    z:=lu[a];
    lu[a]:=lu[b];
    lu[b]:=z;
  end;
  
begin
  losbar:=true;
  lu:=A;
  for k:=1 to n do begin
    (* Zeilenmaximum bestimmen und Permutation initialisieren *)
    permutation[k]:=k;
    zeilenmaximum:=0.0;
    for j:=1 to n do begin
      abslukj:=abs(lu[k,j]);
      if zeilenmaximum < abslukj then zeilenmaximum:=abslukj;
    end;
    if zeilenmaximum = 0.0 then losbar:=false; (* singulär *)
    d[k]:=zeilenmaximum;
  end;
  if losbar then begin
    k:=1;
    while losbar
      and (k <= pred(n)) do begin
      piv:=abs(lu[k,k])/d[k];
      j0:=k;
      (* suche aktuelles Pivotelement *)
      for j:=succ(k) to n do begin
        temp:=abs(lu[j,k])/d[j];
        if piv < temp then begin
          piv:=temp; (* piv enthält das Pivotelement *)
          j0:=j; (* j0 dessen Index *)
        end;
      end;
      if piv < epsilon then losbar:=false
      else begin
        if j0 <> k then begin
          tausch(j0,k);
        end;
        lukk:=lu[k,k];
        for j:=succ(k) to n do begin
          if lu[j,k] <> 0.0 then begin
            lu[j,k]:=lu[j,k]/lukk;
            temp:=lu[j,k];
            for m:=succ(k) to n do begin
              lu[j,m]:=lu[j,m]-temp*lu[k,m];
            end;
          end; (* if *)
        end; (* for j *)
      end; (* else *)
      inc(k);
    end; (* while *)
    if abs(lu[n,n]) < epsilon then losbar:=false
    else begin
      for i:=1 to n do begin
        for k:=1 to n do begin
          if i = permutation[k] then Ainv[k,i]:=1.0 (* Einheitsmatrix *)
                                else Ainv[k,i]:=0.0;
          if k > 1 then begin
            for j:=1 to pred(k) do begin
              Ainv[k,i]:=Ainv[k,i]-lu[k,j]*Ainv[j,i];
            end;
          end;
        end;
        for k:=n downto 1 do begin
          summe:=0.0;
          for j:=succ(k) to n do begin
            summe:=summe+lu[k,j]*Ainv[j,i]
          end;
          Ainv[k,i]:=(Ainv[k,i]-summe)/lu[k,k];
        end;
      end;
    end; (* else *)
  end; (* if *)
end; (* gauss *)

procedure matrixVectorProduct(n:vectorIndex; (* Ordnung *)
                              const A:matrix; (* n x n Matrix *)
                              const x:vector; (* Eingabe-Vektor x *)
                              var y:vector); (* Ausgabe-Vektor y=Ax *)
var i,j:vectorIndex;
    s:extended;
begin
  for i:=1 to n do begin
    s:=0;
    for j:=1 to n do begin
      s:=s+x[j]*A[i,j];
    end;
    y[i]:=s;
  end;
end;