Entwickler-Ecke :: Thema anzeigen - OpenGL kompatible Matrix-Multiplikation mit SSE gesucht

Delete - Sa 20.05.06 23:46
Titel: OpenGL kompatible Matrix-Multiplikation mit SSE gesucht

Hi!

Da ich demnächst ja doch noch mal meine ganzen OpenGL Komponenten veröffentlichen wollte, bin ich gerade dabei, sie etwas aufzupolieren. Wichtig wäre mir natürlich eine SSE Optimierung, weil die ja heutzutage praktisch schon Standard ist.

Was mir vorallem mal fehlt wäre eine Matrix-Matrix Multiplikation.
Eigentlich muss es ja sowas haufenweise geben. Trotzdem habe ich aber bei Google nix passendes finden können...

Problem eigentlich ganz klassisch. Zwei 4x4 Matrizen in Spaltenform...

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:

TGLVectorf4 = array[0..3] of TGLFloat;
TGLMatrixf4 = array[0..3] of TGLVectorf4;

Procedure MatrixMultiply(const M1,M2:TGLMatrixf4; out Result:TGLMatrixf4);
begin
 Result[0,0]:=M1[0,0]*M2[0,0]+M1[1,0]*M2[0,1]+M1[2,0]*M2[0,2]+M1[3,0]*M2[0,3];
 Result[0,1]:=M1[0,1]*M2[0,0]+M1[1,1]*M2[0,1]+M1[2,1]*M2[0,2]+M1[3,1]*M2[0,3];
 Result[0,2]:=M1[0,2]*M2[0,0]+M1[1,2]*M2[0,1]+M1[2,2]*M2[0,2]+M1[3,2]*M2[0,3];
 Result[0,3]:=M1[0,3]*M2[0,0]+M1[1,3]*M2[0,1]+M1[2,3]*M2[0,2]+M1[3,3]*M2[0,3];
 Result[1,0]:=M1[0,0]*M2[1,0]+M1[1,0]*M2[1,1]+M1[2,0]*M2[1,2]+M1[3,0]*M2[1,3];
 Result[1,1]:=M1[0,1]*M2[1,0]+M1[1,1]*M2[1,1]+M1[2,1]*M2[1,2]+M1[3,1]*M2[1,3];
 Result[1,2]:=M1[0,2]*M2[1,0]+M1[1,2]*M2[1,1]+M1[2,2]*M2[1,2]+M1[3,2]*M2[1,3];
 Result[1,3]:=M1[0,3]*M2[1,0]+M1[1,3]*M2[1,1]+M1[2,3]*M2[1,2]+M1[3,3]*M2[1,3];
 Result[2,0]:=M1[0,0]*M2[2,0]+M1[1,0]*M2[2,1]+M1[2,0]*M2[2,2]+M1[3,0]*M2[2,3];
 Result[2,1]:=M1[0,1]*M2[2,0]+M1[1,1]*M2[2,1]+M1[2,1]*M2[2,2]+M1[3,1]*M2[2,3];
 Result[2,2]:=M1[0,2]*M2[2,0]+M1[1,2]*M2[2,1]+M1[2,2]*M2[2,2]+M1[3,2]*M2[2,3];
 Result[2,3]:=M1[0,3]*M2[2,0]+M1[1,3]*M2[2,1]+M1[2,3]*M2[2,2]+M1[3,3]*M2[2,3];
 Result[3,0]:=M1[0,0]*M2[3,0]+M1[1,0]*M2[3,1]+M1[2,0]*M2[3,2]+M1[3,0]*M2[3,3];
 Result[3,1]:=M1[0,1]*M2[3,0]+M1[1,1]*M2[3,1]+M1[2,1]*M2[3,2]+M1[3,1]*M2[3,3];
 Result[3,2]:=M1[0,2]*M2[3,0]+M1[1,2]*M2[3,1]+M1[2,2]*M2[3,2]+M1[3,2]*M2[3,3];
 Result[3,3]:=M1[0,3]*M2[3,0]+M1[1,3]*M2[3,1]+M1[2,3]*M2[3,2]+M1[3,3]*M2[3,3];
end;

So sieht das als normaler Delphicode aus. Ich bräuchte jetzt halt den Assembler-Code der SSE nutzt.
Für Vektoren und ähnliches bekomme ich das auch selber hin, aber das hier wird mir etwas zu kompliziert.

Wenn irgendwo jemand passenden Code gesehen hat (kann auch innerhalb C oder Java sein, ist ja alles asm), dann wäre es toll, wenn er es kurz posten könnte.

THX!

Gruß
Brainiac

Delete - Sa 27.05.06 11:41

Also ich hab's jetzt doch einfach selbst geschrieben (und gerade mal 15 Minuten gebraucht).
Eigentlich ist das ganze ja doch recht einfach wenn man sich es genau ansieht.

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:

procedure MatrixMultiplySSE(const M1,M2:TGLMatrixf4; out Result:TGLMatrixf4);
asm
 // load M1 into xxm entirely as we will need it more than once
 movups   xmm4, [eax] // movaps
 movups   xmm5, [eax+$10]
 movups   xmm6, [eax+$20]
 movups   xmm7, [eax+$30]
 // compute Result[0]
 movss    xmm0, [edx]
 shufps   xmm0, xmm0, $00 //xmm0 = 4x M2[0,0]
 mulps    xmm0, xmm4
 movss    xmm1, [edx+$04]
 shufps   xmm1, xmm1, $00 //xmm1 = 4x M2[0,1]
 mulps    xmm1, xmm5
 addps    xmm0, xmm1
 movss    xmm1, [edx+$08]
 shufps   xmm1, xmm1, $00 //xmm1 = 4x M2[0,2]
 mulps    xmm1, xmm6
 addps    xmm0, xmm1
 movss    xmm1, [edx+$0c]
 shufps   xmm1, xmm1, $00 //xmm1 = 4x M2[0,3]
 mulps    xmm1, xmm7
 addps    xmm0, xmm1
 movups   [ecx], xmm0 // movntps
 // compute Result[1]
 movss    xmm0, [edx+$10]
 shufps   xmm0, xmm0, $00
 mulps    xmm0, xmm4
 movss    xmm1, [edx+$14]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm5
 addps    xmm0, xmm1
 movss    xmm1, [edx+$18]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm6
 addps    xmm0, xmm1
 movss    xmm1, [edx+$1c]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm7
 addps    xmm0, xmm1
 movups   [ecx+$10], xmm0
 // compute Result[2]
 movss    xmm0, [edx+$20]
 shufps   xmm0, xmm0, $00
 mulps    xmm0, xmm4
 movss    xmm1, [edx+$24]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm5
 addps    xmm0, xmm1
 movss    xmm1, [edx+$28]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm6
 addps    xmm0, xmm1
 movss    xmm1, [edx+$2c]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm7
 addps    xmm0, xmm1
 movups   [ecx+$20], xmm0
 // compute Result[3]
 movss    xmm0, [edx+$30]
 shufps   xmm0, xmm0, $00
 mulps    xmm0, xmm4
 movss    xmm1, [edx+$34]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm5
 addps    xmm0, xmm1
 movss    xmm1, [edx+$38]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm6
 addps    xmm0, xmm1
 movss    xmm1, [edx+$3c]
 shufps   xmm1, xmm1, $00
 mulps    xmm1, xmm7
 addps    xmm0, xmm1
 movups   [ecx+$30], xmm0
end;

Trotzdem danke für den Link.
Auf Prefetching hab ich jetzt mal verzichtet. Ich hab zum testen 1000000x willkürliche Matrizen aus einem Array miteinander multiplizieren lassen und die Zeit in Takten gestoppt.
Mit prefetching hab ich dabei eher eine negative Auswirkung feststellen können.
Scheinbar weiß mein P4 wohl bei solch einfachen aufgaben am besten selbst, was er als nächstes braucht.

Was noch performance bringen könnte, wäre wenn die Matrizen 'aligned' im Speicher liegen würden, also die Adressen durch 64 teilbar wären. Weiß jemand, wie man Delphi dazu bringen kann, sowas korrekt anzulegen???

Horst_H - Sa 27.05.06 13:16

Hallo,

reicht nicht teilbar durch 16?
Mit dem FastMemoryManager4 soll dies moeglich sein, bei grossen(>2.5 k) ist dies immer der Fall.

Delphi-Quelltext

1:
2:
3:

{FastMM is actually three memory managers in one: 
small (<2.5K), medium (< 260K) and large (> 260K) blocks 
are managed separately.}

in FastMM4Options.inc

Zitat:

{---------------------------Miscellaneous Options-----------------------------}

{Enable this define to align all blocks on 16 byte boundaries so aligned SSE
instructions can be used safely. If this option is disabled then some of the
smallest block sizes will be 8-byte aligned instead which may result in a
reduction in memory usage.
Medium and large blocks are always 16-byte aligned irrespective of this setting.}
{.$define Align16Bytes}

http://sourceforge.net/project/showfiles.php?group_id=130631

Gruss Horst

Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2024 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!