Algorithmus wie beim einfachen Gauß-Algorithmus. Um numerische Fehler und Divisionen durch 0 möglichst zu vermeiden wird eine Pivotisierung angewandt (Wikipedia, Uni-Bielefeld). Im standard Gauß-Algorithmus wurde in jedem Schritt durch ein Element auf der Hauptdiagonalen der Matrix dividiert. Ist dieser Eintrag allerdings eine Null, dann würde die Division einen Fehler verursachen. Die Idee der Pivotisierung ist nun die Zeilen/Spalten der Matrix geschickt zu vertauschen, um ein möglichst großes Element auf der Hauptdiagonalen stehen zu haben, durch das dann dividiert wird. Das Vertauschen von Zeilen/Spalten eine Matrix (A|b) ist möglich, denn dies ändert die Lösung des Gleichungssystems Ax=b nicht. Mathematisch würde man diese Pivotstrategie durch eine Matrix P darstellen, mit der die Matrix A und die rechte Seite b multipliziert werden:
Ax=b ⇔ (P*A)x = P*b

Als Pivotelement wird in meiner Implementierung das betragsmäßig größte Element eine Spalte ausgewählt. Andere Ansätze wählen das größte einer Zeile oder sogar der gesamten Matrix aus. Beim letzten Fall ist der Aufwand zur Wahl des Pivotelements aber sehr groß, so dass dieser Fall eher selten Verwendung findet. Die Zeilen der Matrix werden dann aber nicht direkt vertauscht, sondern es wird ein Indexvektor mitgeführt, in dem die beiden zugehörigen Elemente dann vertauscht werden.

• A = input(' Matrix A eingeben: ');
• b = input(' Vektor b eingeben: ');
•  
• n=length(b);
•  
• %Permutationsvektor
• p=[1:n];
•  
• % In Dreiecksform umwandeln
• for k=1:n-1
• % Zeilenpivotisierung
• m=k;
• for l=k:n
• if(abs(A(p(l),k)) > abs(A(p(m),k)))
• m=l;
• end
• end
• if(p(m)~=p(k))
• t = p(k);
• p(k) = p(m);
• p(m) = t;
• end
•  
• % Elemination
• for i=k+1:n
• if(A(p(k),k) == 0)
• stop;
• end
• A(p(i),k) = A(p(i),k) / A(p(k),k);
• b(p(i)) = b(p(i)) - b(p(k))*A(p(i),k);
• for j=k+1:n
• A(p(i),j) = A(p(i),j) - A(p(k),j)*A(p(i),k);
• end
• end
• end
•  
• % Rueckwertssubstitution
• x=zeros(n,1);
• b = b(p(1:n));
• A = A(p(1:n),:);
• for i=n:-1:1
• x(i) = (b(i)-A(i,i+1:n)*x(i+1:n)) / A(i,i);
• end
• x

In Zeile 22/23 bricht der Algorithmus ab, falls alle Spaltenelemente NULL sind. Damit ist das Gleichungssystem nicht mehr eindeutig lösbar, denn die Matrix enthält eine Null-Spalte.