Arbeit auf gemeinsamer Matrix klappt nicht

_Eisi

Hallo alle zusammen,

bedingt durch mein Studium versuche ich gerade einen PCG-Solver zu programmieren.
Dafür benötigt ich auch eine Matrizenmultiplikation.
Diese wird mit mehreren Threads berechnet. Das ganze funktioneirt soweit auch, jeder einzelne Thread berechnet seine Ergbenisse richtig.

Um Speicher zu sparen wollte ich die vorhandenen Matrizen als Referenz übergeben (inkl. der Ergebnismatrix) und hier entsteht vermutlich das Problem. Die richtig berechneten Werte werden nicht in die Ergebnismatrix geschrieben.

mapped_matrix<double> mat_mul(mapped_matrix<double> &mat_a, mapped_matrix<double> &mat_b){
	int n_threads = 4; // Anzahl der Threads die gestratet werden
	mapped_matrix<double> mat_c (mat_b.size1(),mat_a.size2());

	int thread_n_cols = (mat_a.size2()/n_threads); // Bestimmung der Spalten die jeder Thread berechnen muss
	if(thread_n_cols<1)
		thread_n_cols = 1;

	unsigned long int c = 0;
	unsigned long int d = c + thread_n_cols;
	for(int i=0;i<n_threads;i++){
		if(mat_a.size2()==1){
			calc.create_thread(boost::bind (mat_mul_calc, mat_a, mat_b, mat_c, 0, 1));
			break;
		}
		cout << "Thread:" << (i+1) << " ;c_start:" << c << " ;c_end:" << d << endl;
		calc.create_thread(boost::bind (mat_mul_calc, mat_a, mat_b, mat_c, c, d));
		c = d;
		d = c + thread_n_cols;
		if(d>mat_a.size2()){
			calc.create_thread(boost::bind (mat_mul_calc, mat_a, mat_b, mat_c, c, mat_a.size2()));
			break;
		}
	}

	calc.join_all();

	return mat_c;
}

void mat_mul_calc(mapped_matrix<double> &mat_a, mapped_matrix<double> &mat_b,mapped_matrix<double> &mat_c,unsigned long int c_start,unsigned long int c_end){
	for(unsigned long int i = 0; i < mat_c.size1(); i++)
		for(unsigned long int j = c_start; j < c_end; j++)
			for(unsigned long int k = 0; k <= mat_a.size1(); k++)
				mat_c(i,j) += mat_a(i,k) * mat_b(k,j);

}

Die Funktion "matmul" übergibt die richtigen Werte an "main", nur eben von "mat_mul_calc" kommt nicht's zurück.

Die beiden Funktionen sind einem Beispiel aus einem Buch nachempfunden, wo auch mit Referenzen an den eigentlichen Variablen ehrumgearbeitet wird, weshalb ich etwas verwirrt bin warum es hier nicht klappt.

Weis jemand Rat?

Gruß

Bist Du Dir sicher, dass das <= in der for-Schleife so richtig ist? Genaues kann mehr erst sagen, wenn man die implementierung der mapped_matrix-Klasse sieht.

_Eisi

Hi,

erst mal Danke für die Antwort.

Stimmt du hast recht das <= muss ein < sein, dies ändert aber leider nichts an meinem Problem. In der Funktion mat_mul_calc jeweils die mat_c ausgeben lasse, sind alle Ergebnisse berechnet (aber eben nur die, die in diesem einzelen Thread berechnet wurden).

Die Matrix Klasse kommt aus der Boost Bibliothek
http://www.boost.org/doc/libs/1_42_0/libs/numeric/ublas/doc/matrix_sparse.htm

Hast du mal (als Test) mat_mul_calc direkt aufgerufen und die Paramter so gesetzt, dass mit dem einem Aufruf die ganze Multiplikation durchgeführt wird? Dann kannst du schon mal ein paar Fehlerquellen ausschliessen.

Ich habe mir die Implementierung von matrix_sparse nicht angesehen, aber sofern das mit listen o.ä. implementiert wird muss der Zugriff gelockt werden, es sei denn boost garantiert thread-safeness.

otze

mapped_matrix ist nicht threadsafe. Allein das macht das schon kaputt. Insbesondere weil du die Threads über die Spalten aufteilst und nicht über die Zeilen und dort permanent Konflikte raufbeschwörst.

Es wäre eh viel schlauer, in der Funktion mat_mul_calc die Zeilen zu splitten, denn die mapped_matrix<double> ist nichts anderes als:

std::vector<std::map<int,double> >

von der Geschwindigkeit her wäre es also bereits unglaublich viel schlauer, den ersten Index zu splitten und dann über die map zu iterieren. Da du dann sagen kannst, dass immer nur ein Thread auf eine map zugreift, kannst du - vielleicht - deine Threadingprobleme los werden. Ich kenne aber die Implementation nicht genauer, ich garantiere für nichts.

Nebenbei ist es dumm, bei einer sparse matrix über alle elemente zu iterieren. Wenn deine Matrix nicht sparse ist, bist du eh Faktor 100 schneller, wenn du auf matrix<double> umsattelst und brauchst die Threads dort gar nicht. Ausserdem kannst du schauen, ob du nicht mit subrange und dem normalen prod in der Funktion arbeiten kannst.

Wventuell solltest du für jeden thread eine eigene mat_c verwenden und die Teile am Ende zusammenfügen. Das würde deine Threadingprobleme auch lösen.

_Eisi

Danke für die vielen Vorschläge, ich werd das die nächsten Tage mal alles durchspielen. Und meld mich dann wieder

EDIT:
Also mit dem direkten Aufruf (und ohne Threads) klappt das ganze.
Liege ich dann mit der Vermutung richtig (wie bereits erwähnt wurde), dass es an der Matrixklasse liegt? Ich werd mir hierzu mal was einfallen lassen.

Das Verwenden einer Sparse-Matrix zusammen mit mehreren Threads macht nur dann Sinn, wenn das Ergebnis wieder eine Sparse Matrix ist. Ansonsten lieber vector<vector<double>> oder etwas aus ublas verwenden.

krümelkacker

Sehr viele Fehler sehe ich da.

Erstmal frage ich mich, wozu Du meinst, zwei Matritzen multiplizieren zu müssen. Der CG-Algorithmus benötigt nur Matrix-Vektor-Multiplikationen.

Dann teilst Du die Arbeit ungeschickt auf. Statt den Threads verschiedene Spalten der Matrix A zuzuordnen, solltest Du die Zeilen aufteilen.

Dann hast Du einen nebenläufigen Zugruff auf eine mapped_matrix über den operator(). Wer sagt, dass das erlaubt ist? Ist es nicht. Du hast ein Datenrennen. operator() bei mapped_matrix verändert eine interne Baum-Datenstruktur. Das geht nicht nebenläufig ohne Konflikte. Außerdem ist der Zugriff über operator() für eine Multiplikation viel zu langsam. Dazu gibt es ublas::axpy, welches die dünne Struktur vernünftig ausnutzen kann. Um eine Teilmatrix auszuwählen, kannst Du einen proxy nehmen, siehe matrix-proxy Header.

_Eisi

Hi, erstmal danke für die Tipps.

Ich meine zwei Matrizen multiplizeiern zu müssen, weil ein Vektor-Matrix Produkt auch nix anderes ist und ich schlicht weg erst mal sehen wollt ob das so funktioniert. Dann bin ich auf das genannte Problem gestoßen und ich bezweifle dass ein Vektor-Matrix Produkt etwas anderes geliefert hätte.

Um den Zugriff auf die Matrix-Daten werd ich mich kümmern.

Gruß

otze

_Eisi schrieb:

Ich meine zwei Matrizen multiplizeiern zu müssen, weil ein Vektor-Matrix Produkt auch nix anderes ist

Vektor-Matrix ist viel viel einfacher. Bei Matrix-Matrix kannst du 1000 Dinge falsch machen (tust du ja gerade).