Matrix invertieren mit ublas

Cordula

Hallo zusammen

Ich bin dabei eine Matrix über die LU-Zerlegung von boost zu invertieren (Code hier: http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?LU_Matrix_Inversion)

Nun, in der Zeile

permutation_matrix<complex<long double>> pmatrix(input_matrix.size1());

// Create a working copy of the input matrix
mapped_matrix<complex<long double> > dummy_matrix (input_matrix);

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

bekomme ich folgenden Fehler zurück:

error C2146: syntax error : missing ';' before identifier 'pm'
error C3861: 'pm': identifier not found

Da ich jetzt schon paar Stunden daran sitze und nichts finde, dacht ich es wäre sinnvoll, wenn noch jemand drüber schaut

Viele Dank schonmal.

Bevor andere Fragen kommen: Das Invertieren lässt sich nicht umgehen.

Welchen Typ soll pm haben? pmatrix selber ist ja eine Variable, die kannst du nicht als Typ angeben.

krümelkacker

Übrigens: Don't invert that matrix.

Cordula

Ja, das ist mir klar. Aber wo das genau so in dem verfügbaren Code steht .. hab ich es einfach übernommen. Wie sollte es denn richtig heißen?

typedef permutation_matrix<std::size_t> pmatrix;
 	// create a working copy of the input
 	matrix<T> A(input);
 	// create a permutation matrix for the LU-factorization
 	pmatrix pm(A.size1());

 	// perform LU-factorization
 	int res = lu_factorize(A,pm);
        if( res != 0 ) return false;

 	// create identity matrix of "inverse"
 	inverse.assign(ublas::identity_matrix<T>(A.size1()));

 	// backsubstitute to get the inverse
 	lu_substitute(A, pm, inverse);

Oder soll das irgendwas anderes bedeuten?

Übrigens: Ich brauche die Elemente der invertierten Matrix, die bekomme ich nicht wenn ich die Matrix nciht invertiere.

krümelkacker

Im zitierten "offiziellen" Beispiel ist pmatrix ein typedef. Bei Dir war es ein Objekt.

> Übrigens: Ich brauche die Elemente der invertierten Matrix, die bekomme ich nicht wenn ich die Matrix nciht invertiere.

Ist klar. Man kommt nicht immer drum herum. Ich muss hier auch gerade die Diagonale einer Inversen einer positiv definiten Matrix berechnen lassen. Die ist aber leider etwas größer, so dass ich das iterativ über den CG-Algorithmus mache.

Cordula

Wow, da hast du dir was vorgenommen, wa ..

Danke, wirklich nicht der Typ Fehler aus denen man lernt Da ist noch was:

Meine zu invertierende Matrix ist komplex, und da bekomme ich Probleme mit der LU-Factorization

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

error C2784: 'bool boost::operator ==(const boost::array<T,N> &,const boost::array<T,N> &)' : could not deduce template argument for 'const boost::array<T,N> &' from 'std::complex<long double>'
see declaration of 'boost::operator =='

Ne Idee wie ich das lösen kann?

Cordula

Weiß wirklich keiner was?

SeppJ

Gib mal ein bisschen Kontext, am besten ein (bis auf den Fehler) compilierbares Minimalbeispiel. Ich kann dir sagen, was der Fehler bedeutet (das kann Google auch), aber ohne ein bisschen mehr Hintergrund nichts über die genaue Ursache sagen.

Cordula

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?

pumuckl

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

krümelkacker

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.

Cordula

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;
}

Braunstein

Hallo,

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

Cordula

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

Cordula

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?

Soley

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

Cordula

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 ..

Soley

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