C++ zu Delphi - (Inverse Matrix) liefert falsches Ergebnis

[MisterNiceGuy]

Hi, ich habe gerade versucht eine Funktion von C++ in Delphi zu übersetzen. Leider liefert mir die Delphi-Übersetzung ein falsches Ergebnis
obwohl ich mir sicher bin alles richtig gemacht zu haben.

//edit: Gelöst, ich habe die repeat-Schleife einmal zu wenig durchlaufen -.- Der Code unten ist korrigiert und steht zur allgemeinen Verfügung bereit! P.S. C++ Quellcode ist von zeiner.at

Findet ihr den Fehler?

C++ Code:

Code:

int Inverse (double a[][NMAX], double ainv[][NMAX], int n) 
{
  int   i, j;                    // Zeile, Spalte
  int   s;                       // Elimininationsschritt
  int   pzeile;                  // Pivotzeile
  int   fehler = 0;              // Fehlerflag
  double f;                      // Multiplikationsfaktor
  const double Epsilon = 0.01;   // Genauigkeit
  double Maximum;                // Zeilenpivotisierung
  extern FILE *fout;
  int pivot = 1;

  // ergänze die Matrix a um eine Einheitsmatrix (rechts anhängen)
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
    {
      a[i][n+j] = 0.0;
      if (i == j)
        a[i][n+j] = 1.0;
    }
  }
  #if DEBUG
     MatOut (stdout, a, n, 2*n);
  #endif
  
  // die einzelnen Eliminationsschritte
  s = 0;
  do {
     // Pivotisierung vermeidet unnötigen Abbruch bei einer Null in der Diagnonalen und
      // erhöht die Rechengenauigkeit
    Maximum = fabs(a[s][s]);
    if (pivot)
    {
       pzeile = s ; 
       for (i = s+1; i < n; i++)
         if (fabs(a[i][s]) > Maximum) {
           Maximum = fabs(a[i][s]) ;
           pzeile = i;
         }
    }
    fehler = (Maximum < Epsilon);

    if (fehler) break;           // nicht lösbar 

    if (pivot)
    {
      if (pzeile != s)  // falls erforderlich, Zeilen tauschen
      { double h;
        for (j = s ; j < 2*n; j++) {
          h = a[s][j];
          a[s][j] = a[pzeile][j];
          a[pzeile][j]= h;
        }
      }
    }

    // Eliminationszeile durch Pivot-Koeffizienten f = a[s][s] dividieren
    f = a[s][s];
    for (j = s; j < 2*n; j++)
      a[s][j] = a[s][j] / f;

    // Elimination --> Nullen in Spalte s oberhalb und unterhalb der Diagonalen
    // durch Addition der mit f multiplizierten Zeile s zur jeweiligen Zeile i
    for (i = 0; i < n; i++ ) {
      if (i != s) 
      {
        f = -a[i][s];                 // Multiplikationsfaktor
        for (j = s; j < 2*n ; j++)    // die einzelnen Spalten
          a[i][j] += f*a[s][j];       // Addition der Zeilen i, s
      }
    }
    #if DEBUG
      fprintf(stdout, "Nach %1i-tem Eliminationschritt:\n", s+1);
      MatOut (stdout, a, n, 2*n);
    #endif
    s++;
  } while ( s < n ) ;

  if (fehler) 
  {
    fprintf (fout, "Inverse: Matrix ist singulär\n");
    return 0; 
  }
  // Die angehängte Einheitsmatrix Matrix hat sich jetzt in die inverse Matrix umgewandelt
  // Umkopieren auf die Zielmatrix
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      ainv[i][j] = a[i][n+j];
    }
  }
  return 1;  
}

Mein Delphi-Code:

Delphi-Quellcode:

function InverseMatrix(Matrix:TMatrix):TMatrix;
  var i,j,s,pzeile,fehler,n:integer;
      f,Maximum,h:double;
      tempMatrix:TMatrix;
  const Epsilon = 0.01;
begin
  //Using this matrix, we can only form the inverse element if height and
  //width are the same
  if length(Matrix) <> length(Matrix[0]) then
  begin
    result := nil;
    exit;
  end;

  fehler := 0;
  n := length(Matrix);
  setlength(result,n,n);
  setlength(tempmatrix,n,2*n);

  tempMatrix := AddIdentityMatrix(Matrix);

  //Elimination - gehe jede Zeile der Matrix durch...
  s := 0;
  repeat
  begin
    //Nimm den Betrag von Element [s][s]
    Maximum := abs(tempMatrix[s][s]);
    //Weise pzeile die aktuelle Zeile zu
    pzeile := s;
    for i := (s+1) to (n-1) do //gehe alle kommenden zeilen durch...
    begin
      if (abs(tempMatrix[i][s]) > Maximum) then //wenn der betrag des wertes in zeile i und Spalte s größer ist
                                                //als das aktuelle maximum dann weise das neue Maximum zu
                                                //und mache bei zeile i weiter
      begin
        Maximum := abs(tempMatrix[i][s]);
        pzeile := i;
      end;
    end;
    //Nun haben wir die Zeile mit dem höchsten Betrag an Stelle s

    //Wenn das Maximum kleiner als 0.01 ist, dann beende weil determinate 0 ist
    if Maximum < Epsilon then
    begin
      fehler := 1;
      break;
    end;

    //Nur wenn die Zeile mit dem höchsten Betrag an Stelle s nicht die Zeile s ist, dann tausche
    //alle Zahlen aus Zeile s mit Zeile pzeile
    if (pzeile <> s) then
    begin
      h := 0;
      for j := s to (2*n)-1 do  //gehe alle zeichen in der aktuellen spalte durch
      begin
        h := tempMatrix[s][j]; //h := tempMatrix[s][j];
        tempMatrix[s][j] := tempMatrix[pzeile][j]; //tempMatrix[s][j] := tempMatrix[pzeile][j];
        tempMatrix[pzeile][j] := h;  //tempMatrix[pzeile][j] := h;
      end;
    end;

    //Eliminationszeile durch Pivot-Koeffizienten f = Matrix[s][s] dividieren
    //...
    //Jetzt aus der Zeile mit unserem Maximum alle Werte durch den Wert an Stelle s teilen
    f := tempMatrix[s][s];
    for j := s to (2*n)-1 do
      tempMatrix[s][j] := tempMatrix[s][j] / f;


    // Elimination --> Nullen in Spalte s oberhalb und unterhalb der Diagonalen
    // durch Addition der mit f multiplizierten Zeile s zur jeweiligen Zeile i
    //...
    //Addiere in jeder Zeile - außer Zeile s - auf jedes Zeichen den wert [s][j]*f
    for i := 0 to n-1 do
    begin
      if( i <> s) then
      begin
        f := -tempMatrix[i][s];
        for j := s to (2*n)-1 do
          tempMatrix[i][j] := tempMatrix[i][j] + (f*tempMatrix[s][j]);
      end;
    end;

    //showmessage('blubb');

    showmessage('Schritt: '+inttostr(s));
    showmatrix(tempmatrix);
    inc(s);
  end until (s = n);

  if fehler = 1 then
  begin
    showmessage('Inverse: Matrix ist singulär!');
    exit;
  end;

  showmatrix(tempmatrix);
  // Die angehängte Einheitsmatrix Matrix hat sich jetzt in die inverse Matrix umgewandelt
  // Umkopieren auf die Zielmatrix
  for i := 0 to n-1 do
    for j := 0 to n-1 do
      result[i][j] := tempMatrix[i][j+n];
end;

P.S. Meine Funktion zum Hinzufügen der Identitätsmatrix funktioniert einwandfrei!