Matrix invertieren in C++



Wie man eine Matrix mit Hilfe des Gauß-Jordan-Algorithmus invertiert habe ich bereits gezeigt. Nun habe ich versucht den vorgestellten Algorithmus in C++ umzusetzen. Als Ergebnis kam eine überschaubare Funktion mit einer Hilfsfunktion um zwei Zeilen in einer Matrix zu vertauschen.

/*
	Vertauscht zwei Zeilen in einer NxM Matrix.

	line1, line2 - Index der Zeilen, die vertauscht werden sollten.

	return: false, falls line1 oder line2 nicht in der Matrix liegen.
*/
template <size_t N,size_t M>
bool swapLine(double mat[N][M], size_t line1, size_t line2)
{
	if(line1 >= N || line2 >= N)
		return false;

	for(size_t i = 0; i < M; ++i)
	{
		double t = mat[line1][i];
		mat[line1][i] = mat[line2][i];
		mat[line2][i] = t;
	}
	
	return true;
}

/*
	Invertiert eine NxN Matrix mit Hilfe des Gauß-Jordan-Algorithmus.

	mat - Matrix die Invertiert werden soll.
	inv - Die Inverse der Matrix mat.

	return: false, falls die Matrix nicht invertierbar ist.
*/
template <size_t N>
bool invertMatrix(const double mat[N][N], double inv[N][N])
{
	// Eine Nx2N Matrix für den Gauß-Jordan-Algorithmus aufbauen
	double A[N][2*N];
	for(size_t i = 0; i < N; ++i)
	{
		for(size_t j = 0; j < N; ++j)
			A[i][j] = mat[i][j];
		for(size_t j = N; j < 2*N; ++j)
			A[i][j] = (i==j-N) ? 1.0 : 0.0;
	}

	// Gauß-Algorithmus.
	for(size_t k = 0; k < N-1; ++k)
	{
		// Zeilen vertauschen, falls das Pivotelement eine Null ist
		if(A[k][k] == 0.0)
		{
			for(size_t i = k+1; i < N; ++i)
			{
				if(A[i][k] != 0.0)
				{
					swapLine<N, 2*N>(A,k,i);
					break;
				}
				else if(i==N-1)
					return false; // Es gibt kein Element != 0
			}
		}

		// Einträge unter dem Pivotelement eliminieren
		for(size_t i = k+1; i < N; ++i)
		{
			double p = A[i][k]/A[k][k];		
			for(size_t j = k; j < 2*N; ++j)
				A[i][j] -= A[k][j]*p;
		}
	}

	// Determinante der Matrix berechnen
	double det = 1.0;
	for(size_t k = 0; k < N; ++k)
		det *= A[k][k];

	if(det == 0.0)	// Determinante ist =0 -> Matrix nicht invertierbar
		return false;

	// Jordan-Teil des Algorithmus durchführen
	for(size_t k = N-1; k > 0; --k)
	{
		for(int i = k-1; i >= 0; --i)
		{
			double p = A[i][k]/A[k][k];
			for(size_t j = k; j < 2*N; ++j)
				A[i][j] -= A[k][j]*p;
		}
	}

	// Einträge in der linker Matrix auf 1 normieren und in inv schreiben
	for(size_t i = 0; i < N; ++i)
	{
		const double f = A[i][i];
		for(size_t j = N; j < 2*N; ++j)
			inv[i][j-N] = A[i][j]/f;
	}

	return true;
}

Ich habe die Funktion mit ein paar Matrizen ausprobiert und die Implementierung scheint zu funktionieren. Weitere Tests wären aber trotzdem zu empfehlen.

Wie schnell die Implementierung ist, habe ich nicht getestet. Es gibt sicherlich Optimierungsmöglichkeiten. Wer Spaß daran hat, kann ja eine verbesserte Version in Kommentaren schreiben.

Eine CPP-Datei mit ein paar Tests könnt ihr hier runterladen.

Viel Spaß damit! =)




2 Kommentare zu “Matrix invertieren in C++”

  1. Freddy Kayam 29. April 2013 um 10:10 Uhr

    Ich bedanke mich recht herzlich!

    Ich hab seit einer Weile nach einem simplen Code zur Inverierung von matrizen gesucht und endlich gefunden!

  2. Maximam 29. April 2013 um 10:37 Uhr

    Bitte! :)

    Falls im Code Fehler drinstecken, dann bitte ich um eine Rückmeldung.

Trackback URI | Kommentare als RSS

Einen Kommentar schreiben

XHTML: Du kannst folgende Tags verwenden: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <sub> <sup>

Hinweis: Ich behalte mir das Recht vor solche Kommentare, die Beleidigungen oder rechtswidrige Inhalte beinhalten erst nach einer Editierung freizugeben oder kommentarlos zu löschen. Ähnliches gilt auch für Kommentare die offensichtlich nur der Suchmaschinenoptimierung dienen.