Multiplikation mit SSE

[michael_kr]

Habe ein kleines Programm geschrieben, welches (n * n) Matrizzen miteinander multipliziert.
Hier die entsprechende Methode:

// berechnet das Ergebnis der Matrixmultiplikation gemäß Variante SSE2
// **matrixSol : Ergebnismatrix
// **matrixA : erste Matrix
// **matrixB : zweite Matrix
// n : Größe der Matrix
float** mulVarSSE2(float **matrixSol, float **matrixA, float **matrixB, int n) {
	
	__m128 matrixASSE; // Matrix A
	__m128 matrixBSSE; // Matrix B
	matrixSol = (float**) calloc(32, n * sizeof(float*));
	for (int i = 0; i < n; i++) {
		matrixSol[i] = (float*) calloc(32, n * sizeof(float*));     // 
	}
	__m128 **matrixSolSSE = (__m128**) matrixSol;
	
	
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {

			matrixASSE = _mm_set1_ps(matrixA[j][i]);
			for (int k = 0; k < n; k += 4) {
				matrixBSSE = _mm_set_ps(matrixB[i][k + 3], matrixB[i][k + 2], matrixB[i][k + 1], matrixB[i][k]);
		
				matrixSolSSE[j][k] += _mm_mul_ps(matrixASSE, matrixBSSE);
			}
		}
	}
	return matrixSol;
}

Habe es versucht so einzurichten, dass die Lösungen direkt in die entsprechende Lösungmatrix geschrieben werden. (Der Rückgabewert ist noch ein Überbleibsel, darum kümmere ich mich später).
Die Werte werden korrekt berechnet, das Problem liegt darin, dass für jede Zeile immer nur die ersten vier Werte geschrieben werden (also nur für k = 0), der Rest der Zeile besteht aus Nullen...

Kann mir jemand sagen, was ich falsch mache? Die Lösungsmatrix (die, die auch zurück gegeben wird) ist korrekt allokiert, für andere Methoden funktioniert sie, nur eben mit dem Pointer der die "128-Bit Zahl" direkt ins Array schiebt haut es nicht hin.

Danke schonmal im Voraus

 

[aroxx]

ich bin mir nicht sicher ob du sowas wie += mit den intel intrinsics machen darfst.
wie liest du denn deine ergebnis matrix aus? kommt da noch irgendwo ein _mm_store_?

vllt solltest das "_mm_set_ps" in Zeile 22 durch ein "_mm_loadu_ps (float const* mem_addr)" ersetzen. ist besserer stil

 

[michael_kr]

Für die ersten vier Werte jeder Zeile funktioniert das.
Die Ergebnis Matrix kann ja über matrixSol[j] ganz normal ausgelesen werden.
Das speichern in die Matrix übernehmen diese Teile:

__m128 **matrixSolSSE = (__m128**) matrixSol;      // Zeile 14

matrixSolSSE[j][k] += _mm_mul_ps(matrixASSE, matrixBSSE);

Dieses Konstrukt in Zeile 14 hab ich hier gefunden.

 

[aroxx]

#ifndef TIME_SSE
		float xFloat = 1.0f;
		for (int i=0 ; i < length; i++)
		{

			pResult[i] = sqrt(xFloat) / xFloat;	// Even though division is slow, there are no intrinsic functions like there are in SSE
			xFloat += 1.0f;
		}

#endif	// !TIME_SSE

das ist die einzige Stelle wo ich auf der von dir angegeben Seite ein += gefunden habe. das ist floating point addition und nicht sse-register addition.

versuch mal dein

matrixSolSSE[j][k] += _mm_mul_ps(matrixASSE, matrixBSSE);

durch

matrixSolSSE[j][k] = _mm_add_ps(matrixSolSSE[j][k], _mm_mul_ps(matrixASSE, matrixBSSE));

zu ersetzen. ich hoffe dass dir das nicht um die ohren fliegt... wäre heute nicht weihnachten würde ich etwas ausführlicher testen

 

[michael_kr]

das ist mal ne gute Idee. Werde es versuchen

Ergänzung (24. Dezember 2014)

Selbes Ergebnis...

 

[aroxx]

ok dann nur ne einfache debugging idee: versuch das ganze mal mit 16x16 einheitsmatrizen oder mit matrizen, die du von 1-256.0f vorinitialisierst. dann lass dir mal den inhalt deiner sseregister ausgeben und vergleiche, ob das mit deiner papierrechnung übereinstimmt. wenns nicht ganz eilig ist schau ich mir das morgen gern mal genauer an

 

[michael_kr]

Lasse mir die Lösungen mehrerer Varianten ausgeben, die stimmen. Die hatte ich auf Papier ausprobiert. Die ersten vier Werte der SSE2 Variante stimmen auch mit den anderen überein. Für größere Matrizzen (hab es für n = 16 auch mal ausgeben lassen) füllt er auch nur die ersten vier Werte.

Wünsche ein frohes Fest und danke für deine Hilfe an Weihnachten!

Ergänzung (24. Dezember 2014)

Und das eilt nicht, reicht also morgen vollkommen aus

 

[aroxx]

Ok dann morgen Obwohl ich mega neugierig bin was da kaputt geht
Nur noch eine anmerkung: dein code geht vermutlich für matrizen, deren größe kein Vielfaches von 4 ist, kaputt. ist das egal?

 

[michael_kr]

Ja funktioniert nur für vielfache von 4
Bin gerade nicht am PC, vllt poste ich später noch die Konsolen Ausgabe

Ergänzung (25. Dezember 2014)

Hier die versprochene Konsolen Ausgabe:

 

[aroxx]

hab dir meine mailadresse per pm geschrieben. schick mir bitte mal den ganzen code zu (inkl main und i/o), dann versuch ich das problem zu finden

 

[Hancock]

Hmm, also ich seh da n paar Fehler, schon mathematisch ist der Code nicht richtig.

Der Originalalgorithmus geht ja wie

for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
        for (int k = 0; k < n; ++k)
            c[i][j] += a[i][k] + b[k][j];

Die äußeren Schleifen hast du dann ja richtig umgesetzt, aber dann wir es wild.

Die innerste Schleife ist ein Skalarprodukt, das versuchst du zu programmieren, vergisst aber, die Parameter aus A jeweils mitzunehmen und denkst, dass _mm_mul_pd ein 4-er Skalarprodukt macht.

Schau mer mal das innere der i,j Schleifen an:

matrixASSE = _mm_set1_ps(matrixA[j][i]);
for (int k = 0; k < n; k += 4) {
    matrixBSSE = _mm_set_ps(matrixB[i][k + 3], matrixB[i][k + 2], matrixB[i][k + 1], matrixB[i][k]);

    matrixSolSSE[j][k] += _mm_mul_ps(matrixASSE, matrixBSSE);
}

Der Trivialcode ohne SSE wäre (du scheinst andersrum zu indizieren):

for(int k=0;k<n;++k)
    matrixSolSSE[j][k] += matrixA[k][i]*matrixB[j][k];

Jetzt mit SSE, 4 Multiplikationen auf einmal (in 4 Speicherstellen und am Ende zusammensummieren):

//Zeile * Spalte (Indizes vertauscht => A*B = (B^T*A^T)^T => Matrizen tauschen)
register __m128 buf = _mm_setzero_ps();
for (int k = 0; k < n; k += 4) {
    ma = _mm_loadu_ps(&(matrixA[i][k]));
    mb = _mm_set_ps(matrixB[k][j], matrixB[k+1][j], matrixB[k+2][j], matrixB[k+3][j]);
    buf=_mm_add_ps(buf, _mm_mul_ps(ma, mb));
}
matrixSol[i][j] = buf.m128_f32[0] + buf.m128_f32[1] + buf.m128_f32[2] + buf.m128_f32[3];

Das sollte funktionieren.

Noch ne Anmerkung: Du kannst den ersten Parameter so nicht als Rückgabewert verwenden, weil du matrixSol überschreibst, was der Caller aber nicht sieht. (das ganze _alloc mit einem if(!matrixSol) umklammern).