3. Matrix invertieren mit ublas

Matrix invertieren mit ublas

  • C

    Nee, ich weiß was der Fehler ist und woran er liegt - ich weiß nur keinen Ausweg 🙂

    Also, weiter oben habe ich den Code eingfügt, ich will eine Matrix über LU-Zerlegung invertieren. Der Code lässt sich zu Zeile 9 kompilieren:

    // Perform LU-factorization
    int res = lu_factorize(A, pm);

    Der Grund dafür ist, dass die LU-Zerlegung nur für reelle Matrizen definiert ist und meine zu invertierende eine komplexe Matrix ist. Und da weiß ich nicht wie ich weiter vorgehen kann.

    Also .. jemand eine Idee dazu?

    0
  • P

    Dann musst du einen anderen Algo nehmen, um die Matrix zu invertieren. Wofür brauchst du denn die Elemente der invertierten Matrix?

    0
  • K

    Cordula schrieb:

    Der Grund dafür ist, dass die LU-Zerlegung nur für reelle Matrizen definiert ist

    Wenn Du mit "LU-Zerlegung" die Implementierung einer LU-Zerlegung meinst, mag das sein. Aber komplexe Matritzen lassen sich natürlich auch so zerlegen.

    Cordula schrieb:

    Also .. jemand eine Idee dazu?

    Zur Fehlermeldung kann ich nicht viel sagen, da fehlt mir der Kontext. Es scheint so, als ob irgendwo ein Vergleichsoperator benötigt wird, der nicht vorhanden ist. Das kann ich mir aber ohne weiteren Kontext nicht erklären, da std::complex<T> zumindest == und != anbietet. Du hast also entweder etwas falsch gemacht oder die Implementierung unterstützt keine komplexen Zahlen. Ich würde auf ersteres Tippen. Falls ich damit Recht haben sollte, wär es von Vorteil, wenn Du selbst versuchen würdest, ein neues, kleines Programm zu erzeugen, was nur eine LU-Zerlegung einer komplexen Matrix berechnet und es dann komplett zeigst, falls es nicht läuft. Liegt der Fehler wirdklich bei lu_factorize, kannst Du auch in der Boost.uBLAS Mailingliste nachfragen, ob jemand einen Patch dafür parat hat.

    0
  • C

    Hallo zusmmen,

    Vielen Dank für eure Antworten, aber es wäre vielleicht sinnvoll wenn ihr das bisher gepostete überfliegen würdet, bevor ihr das Gleiche hinschreibt was ich schon selbst erklärt habe. Das Problem ist dass die ublas-Implementierung der LU-Zerlegung keine komplexen Matrizen unterstützt. Hier hab ich mir ein Work-Around geschrieben, der die komplexe Matrix in Real- und Imaginärteil zerlegt, die Inverse ausrechnet und wieder in einer komplexen Matrix zusammenflickt - für den Fall dass jemand in Zukunft das gleiche Problem hat.

    bool InvertMatrix(mapped_matrix<complex<long double> >& input_matrix, mapped_matrix<complex<long double> > inv_matrix){

	typedef permutation_matrix<long double> pmatrix;

	int Size = input_matrix.size1();

	// To invert the complex input matrix: decompose the matrix into a matrix with the double size and structure as following:
	//  Re(Matrix)  Im(Matrix)
	// -Im(Matrix)  Re(Matrix)
	mapped_matrix<long double> Dec_Real_Im(2*Size, 2*Size);

	// Save Real and Imag parts of the entries in own matrices
	// Create matrices
	mapped_matrix<long double > M_Real(Size, Size);
	mapped_matrix<long double > M_Imag(Size, Size);

	// Iterate over entries and fill matrices
	for (int i = 0; i < Size; i++){
		for (int j = 0; j < JacobiSize; j++){
			complex<long double> Entry_Complex = input_matrix(i, j);
			M_Real(i, j) = Entry_Complex.real();
			M_Imag(i, j) = Entry_Complex.imag();
		}
	}

	// Put the Real- and Imag-matrices together in Dec_Real_Im
	// [Re(Matrix)]  Im(Matrix)
	// -Im(Matrix)   Re(Matrix)
	project(Dec_Real_Im, boost::numeric::ublas::slice ((0,0), 1, Size), boost::numeric::ublas::slice ((0,Size-1), 1, Size)) = M_Real;

	//  Re(Matrix)  [Im(Matrix)]
	// -Im(Matrix)   Re(Matrix)
	project(Dec_Real_Im, boost::numeric::ublas::slice ((0,Size), 1, Size), boost::numeric::ublas::slice ((0,2*Size-1), 1, Size)) = M_Imag;

	//  Re(Matrix)    Im(Matrix)
	// [-Im(Matrix)]  Re(Matrix)
	project(Dec_Real_Im, boost::numeric::ublas::slice ((Size,0), 1, Size), boost::numeric::ublas::slice ((Size-1,0), 1, Size)) = M_Imag*(-1.L);

	//  Re(Matrix)   Im(Matrix)
	// -Im(Matrix)  [Re(Matrix)]
	project(Dec_Real_Im, boost::numeric::ublas::slice ((Size,Size), 1, Size), boost::numeric::ublas::slice ((Size,2*Size-1), 1, Size)) = M_Real;

	// Create a permutation matrix for the LU-factorization
	pmatrix pm(Dec_Real_Im.size1());

	// Create a temporary matrix for the inverse of the extended matrix
	mapped_matrix<long double> inv_matrix_temp (2*Size, 2*Size);

	// Perform LU-factorization
	int res = lu_factorize(Dec_Real_Im, pm);
	if (res != 0){
		return false;
	}

	// Create identitiy matrix of 'inverse'
	inv_matrix_temp.assign(identity_matrix<long double >(Dec_Real_Im.size1()));

	// Backsubstitute to get the inverse
	lu_substitute(Dec_Real_Im, pm, inv_matrix_temp);

	// Cut the upper half of the matrix and put the complex values together
	// !! Operation Reduction possible !!
	//sensitivity_matrix = project(inv_matrix, boost::numeric::ublas::slice((0,), 1, Size)
	for (int i = 0; i < Size, i++){
		for (int j = 0; j < Size, j++){
			long double S_Real = inv_matrix_temp(i, j);
			long double S_Imag = inv_matrix_temp(i, j + Size);
			inv_matrix(i, j) = complex<long double> (S_Real, S_Imag);
		}
	}

	return true;
}
    0
  • B

    Hallo,

    hast du mal das hier gelesen?
    http://lists.boost.org/MailArchives/ublas/2007/09/2299.php

    0
  • C

    Hab ich, danke 🙂 Bei mir ist das Installieren neuer Bibliotheken sehr umständlich, von daher will ich es nach Möglichkeit vermeiden. Aber klar, wenn man's machen kann 🙂

    0
  • C

    Noch einmal der Post:

    Beim Ausführen bekomm ich eine Meldung über einen 'gefährlichen' Cast:

    warning C4244: 'argument' : conversion from 'const boost::numeric::ublas::type_traits<long double>::value_type' to 'unsigned int', possible loss of data

    Diese kommt von der Zeile mit der LU-factorization:

    // Perform LU-factorization
	int res = lu_factorize(*MatrixPointer, pm);

    Heißt das, dass lu_factorize nur für unsigned int definiert ist? Was kann ich sonst für long doubles nehmen?

    0
  • S

    Der Parameter pm ist vom Typ

    permutation_matrix<long double>

    lu_factorize(*MatrixPointer, pm);

    Durch den Aufruf wird der Typ implizit in unsigned int umgewandelt. Da dabei die Nachkommastellen abgeschnitten werden, warnt dich dein Compiler, dass da möglicherweise Daten verloren gehen 😉

    0
  • C

    Stimmt, das hab ich übersehen.

    Ist das nur die permutation_matrix deren Typ umgewandelt wird?

    Wie man bestimmt schon sehen kann, blick ich da nicht so 100% durch .. Also bitte bitte antwortet etwas ausführlicher .. 😃

    0
  • S

    Wie die Fehlermeldung es bereits sagt, ist es nur die permutation_matrix...

    Eigentlich ist das kein besonderer Grund zur Sorge, aber ich würde an deiner Stelle nochmal nachgucken, ob das auch die richtigen Parameter sind, die diese Funktion erwartet. Da ich keine Doku zu dieser Funktion gefunden habe, kann ich das momentan nicht nachprüfen. Denn für mich ist es ziemlich merkwürdig eine Matrix auf einen unsigned int zu casten...

    Vielleicht muss als 2. Parameter etwas völlig anderes übergeben werden, oder die Funktion ist für etwas anderes, als du denkst 😃

    0
Anmelden zum Antworten