C++ Amp Matrix Multiplikation vs MatLab

https://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html

;-)

Würde vermuten dass die Matlab implementierung SIMD Instruktionen wie 
AVX verwendet

Ah klingt interessant, schaue ich mir mal an. Scheint ja zum Standard C 
zu gehören.

Ich meine dieses AVX

Möchte keine andere Bibliothek nehmen, will ja was lernen :D

ERic schrieb:
> Ich meine dieses AVX

AVX ist gewiss nicht Standard-C, aber viele x86-Compiler bieten es als 
Erweiterung an. Es gibt aber schon so viele sehr schnelle Libraries für 
lineare Algebra, welche solche Erweiterungen schon nutzen, dass man 
einen sehr guten Grund braucht um es selbst zu implementieren...

Du kannst dir aber anschauen wie Eigen das implementiert

Ja muss ich wohl nochmal schauen. Es soll später auch auf einem ARM 
Processor laufen und. Will mir nix verbauen, indem ich an irgendwelche 
libs gebunden bin. Allerdings erstmal auf Windows. Werde es jetzt mal 
mit diesen avx probieren, sowas gibt's ja auch für arm, neon glaube 
nennt sich das. Ausserdem kann ich eigenen Code besser anpassen. Das 
Beispiel oben war sicher nicht optimal, dient erstmal um herauszufinden 
wie ich schneller multiplizieren kann. Sind die matrizenn kleiner kann 
mann das amp gleich weglassen, da viel langsamer als Standard. Ich will 
auch nicht besser sein als die fertigen libs, allerdings schneller als 
0,5s.

GRÜßE Eric.

Eric schrieb:
> sowas gibt's ja auch für arm, neon glaube nennt sich das

Naja, das sind halt zwei verschiedene Technologien. ARM hat dazu noch 
verschiedene Vector Extensions, NEON ist nur eine davon. Wenn du deinen 
Code für AVX schreibst, musst du ihn für NEON & Co komplett neu 
schreiben. Der Vorteil von Libraries ist, dass dieser Aufwand da schon 
getan ist. Außerdem ist das ganze Thema schon gut erforscht und eben 
durch Libraries abgedeckt; wer da noch was dran verbessern kann (also 
vermutlich Dr. der Mathematik) würde hier nicht im Forum fragen...

Eigen z.B. ist eine Header only library und auch für ARM optimiert. Es 
gibt absolut keinen Grund da irgendwas selber zum implementieren außer 
für den Lerneffekt.

Ja muss ich sehen und den Aufwand mal abschätzen. Allerdings werde ich 
es einmal selber schreiben, die Woche investiere ich mal (dümmer werde 
ich nicht davon) . Vielleicht reicht mir die Performance dann schon, 
wenn nicht suche ich ne Alternativloesung. So viele Befehle sind es ja 
auch nicht bei den avx modulo4 und paar anweisungen ( zumindest bei 
neon) matrixmultiplikation .  Bin ja nur etechnik-Ing., kenne mich nicht 
so gut aus mit den Rechnerunits, baue die Dinger ja nicht.. Jetzt habe 
ich eh erstmal Urlaub. Trotzdem Danke für die Tipps.



Grüße,

Eric

Habe eben den kompletten Algorithmus schon fertig (etwas komplexer als 
nur die Eigenwerte) , und scheitert nur noch an der 
Ausführungsgeschwindigkeit.
Ja ist ein Cortex m, hat aber keine Floating Point Unit. Na gut, jetzt 
aber Urlaub.

Grüße

Andreas M. schrieb:
> Übrigens, die 0.5s/Multiplikation für 1000x1000 Zufallsmatrizen schafft
> bei mir Python/numpy (vermutlich über libblas) auf einem Core einer 10
> Jahre alten CPU.... (Und zwar ohne AVX)

Glückwunsch! Benutzt numpy bei dir die MKL oder openBLAS? Und hat die 
CPU kein AVX? Oder verbietest du es MKL/BLAS explizit?

Ich machte ein paar mehr Tests. Dieser Code ist schon 100 mal schneller 
als Standard Matrix Multiplikation. Ich glaube wenn man es in Blöcke 
aufteilt gibt es einen deutlichen Performance Sprung. Werde diesen 
ANsatz wohl weiterverfolgen, das AVX ist erscheint mir doch zu komplex. 
Und ich war bis jetzt auch immer langsamer als STD MatrixMultiplikation.


[c]

float** StandardMultTansposeBlocked(float **M1, float **M2,float** 
result)
{
    int i, j, k;
    float sum=0;
    for (i = 0; i < MU; i++)
        for (j = 0; j < MU; j++)
        {
            M2[i][j] = M2[j][i];
        }

    for (i = 0; i < MU/2; i++) {
        for (j = 0; j < MU/2; j++) {
            float* M1ik = &M1[i][0];
            float* M2jk = &M2[j][0];
            for (k = 0; k < MU/2; k++) {
                sum += M1ik[k] * M2jk[k];
            }
            result[i][j] = sum;
        }
    }
    sum = 0;
    for (i = MU/2; i < MU; i++) {
        for (j = MU/2; j < MU; j++) {
            float* M1ik = &M1[i][0];
            float* M2jk = &M2[j][0];
            sum = result[i][j];
            for (k = MU/2; k < MU; k++) {
                sum += M1ik[k] * M2jk[k];
            }
            result[i][j] = sum;
        }
    }


    return result;
}
[\c]

Gruß,

Eric

Eric schrieb:
> Ich machte ein paar mehr Tests. Dieser Code ist schon 100 mal
> schneller
Hast du auch getestet, ob das Ergebnis korrekt ist? Diese Zeilen sehen 
für mich auf den ersten Blick problematisch aus:

1

for (i = 0; i < MU; i++)

2

    for (j = 0; j < MU; j++)

3

        M2[i][j] = M2[j][i];

[c] float** StandardMultTansposeBlocked(float **M1, float **M2,float** 
result)
{
    int i, j, k;
    float sum=0;
    for (i = 0; i < MU; i++)
        for (j = 0; j < MU; j++)
        {
            M2[i][j] = M2[j][i];
        }

    for (i = 0; i < MU/2; i++) {
        for (j = 0; j < MU/2; j++) {
            float* M1ik = &M1[i][0];
            float* M2jk = &M2[j][0];
            for (k = 0; k < MU/2; k++) {
                sum += M1ik[k] * M2jk[k];
            }
            result[i][j] = sum;
        }
    }
    sum = 0;
    for (i = MU/2; i < MU; i++) {
        for (j = MU/2; j < MU; j++) {
            float* M1ik = &M1[i][0];
            float* M2jk = &M2[j][0];
            sum = result[i][j];
            for (k = MU/2; k < MU; k++) {
                sum += M1ik[k] * M2jk[k];
            }
            result[i][j] = sum;
        }
    }


    return result;
}
[\c]

hab es nur getestet für eine 4*4 Matrix das ging

oh ja ist noch ein Bug drin :D Sorry

So der Code ist Bugfrei :D. Ja Geschwindikeitsverbesserung Faktor 2 
gegenüber Transposed Multi. Naja davor hatte ich mich wohl zu Früh 
gefreut:D

1

 for (i = 0; i < MU / 2; i++) {

2

        for (j = 0; j < MU; j++) {

3

            float* M1ik = &M1[i][0];

4

            float* M2jk = &M2[j][0];

5

            sum = 0;

6

            for (k = 0; k < MU/2; k++) {

7

                sum += M1ik[k] * M2jk[k];

8

            }

9

            result[i][j] = sum;

10

        }

11

12

    }

13

    for (i = MU/2; i < MU; i++) {

14

        for (j = 0; j < MU; j++) {

15

            float* M1ik = &M1[i][0];

16

            float* M2jk = &M2[j][0];

17

            sum = result[i][j];

18

            for (k = 0; k < MU/2; k++) {

19

                sum += M1ik[k] * M2jk[k];

20

            }

21

            result[i][j] = sum;

22

        }

23

    }

ah falls das einer testet, die matrix2 muss davor transponiert werden. 
Bei größeren Matrizen ist es sogar Faktor 4 schneller :D. So ab N=1000

Und jetzt die Krönung, super schnell, aber nicht durch getestet

1

int i, j, k;

2

    float sum=0;

3

    for (i = 0; i < MU; i++)

4

        for (j = 0; j < MU; j++)

5

        {

6

           M2[i][j] = M2[j][i];

7

           result[i][j] = 0;

8

        }

9

10

11

int loop = MU / 16;

12

13

    for (int count = 0; count <= MU; count += loop)

14

    {

15

        for (i = count; i < loop + count; i++) {

16

            for (j = 0; j < MU; j++) {

17

                float* M1ik = &M1[i][0];

18

                float* M2jk = &M2[j][0];

19

                sum = 0;

20

                for (k = 0; k < loop; k++) {

21

                    sum += M1ik[k] * M2jk[k];

22

                }

23

                result[i][j] = sum;

24

            }

25

26

        }

27

    }

sum = result[i][j]; nicht sum=0. Habe nur die Geanuigkeit für 4*4 
Matrizen überprüft, und dann einfach nur MU beliebig erhöht aus 512 zum 
Beispiel, da ist dann Faktor 100 drin. Aber ob es stimmt keine AHnung, 
vieleicht wird irgendwas durch den Compiler gut optimiert

nee doch zu früh gefreut, da ist noch etwas faul, es wird nicht alles 
richtig berechnet bei großen Matrizen

Eric schrieb:
> for (i = 0; i < MU; i++)
>         for (j = 0; j < MU; j++)
>            M2[i][j] = M2[j][i];

Wie oben schon geschrieben, das stinkt kilometerweit gegen den Wind.

Eric schrieb:
> ah falls das einer testet, die matrix2 muss davor transponiert werden.
Vorsicht! Transponieren ist eine mathematische Operation. Das ist nur 
zufällig in deinem betrachteten Fall das was benötigt wird. Die Funktion 
verlangt ein spezielles Speicherlayout der Matrizen. Das ist nicht das 
gleiche.

Du solltest dich dringend von deinen double** verabschieden, wenn du 
keinen guten Grund für deren Nutzung hast.

Mit Blöcken holt man mit diesen Code nur 10 Prozent heraus gegenüber 
Transposed. Etwas besser als nur Transposed, aber nicht gut genug für 
mich.
War noch ein Bug, die innerste Schleife muss auch bis MU laufen. Naja, 
10 Prozent sind besser als nix.

Walter T. schrieb:
> mh schrieb:
>> Vorsicht! Transponieren ist eine mathematische Operation.
>
> Naja, in Tensorschreibweise ist es ein Tauschen der Zeilen- und
> Spaltenindizes. Das passt hier schon.

Tensoren haben erstmal nichts mit Zeilen und Spalten zu tun und 
transponieren von Tensoren ist auch interessant. Belassen wir es erstmal 
bei Matrizen.
Du hast auch praktischerweise den Rest meines Beitrags nicht zitiert, 
der erklärt, warum ich darauf hinweise, dass transponieren eine 
mathematische Operation ist.

Walter T. schrieb:
> mh schrieb:
>> Tensoren haben erstmal nichts mit Zeilen und Spalten zu tun
>
> Speicherlayout hat auch nichts mit Zeilen und Spalten zu tun. Zumindest
> für den Nutzer nicht. Für den Kernspeicherwickler schon.

Von was redest du? Was hat ein Kernspeicherentwickler mit 
Matrizenmultiplikation zu tun? Falls du es nicht verstanden hast, mit 
Speicherlayout meine ich die Anordnung der Daten im Speicher, also sowas 
wie Row-Major oder Column-Major

Walter T. schrieb:
> Das war allerdings auch noch zu einer Zeit, als 8 MB knapp nicht in den
> L3 Cache passten (und die Matrizen waren auch ein wenig größer, hatten
> aber dafür eine schöne Bandstruktur).

Wenn es sich um 1000x1000 Matrizen handelt, reichen bei diesem "naiven" 
Ansatz auch 8MB L3 Cache nicht aus, da mindestens 2 dieser Matrizen im 
Cache sitzen wollen. Mal davon abgesehen gewinnt man keinen 
Geschwindigkeitsrekord, wenn man die ganze Zeit auf den L3 angewiesen 
ist...

ich korrigiere nochmal, baut man eine Release Version ist der Code doch 
um einiges schneller. 2,5 s braucht man mit der transposed matrix Mult 
und 0,85s mit der oberen Version. Es war doch keine Zeitverschwendung. 
Allerdings geht das nur in Release Conf.

Testmatrizen waren zwei 1000x1000 Matrizen.

Eric schrieb:
> Testmatrizen waren zwei 1000x1000 Matrizen.

Was ist der Inhalt und wie testest du das Ergebnis?

hab die Matrizen gefüllt mit A[i][j]=(float)i*(float)j Index.

Bin jetzt ungefähr so schnell wie dieses AMP C++, jedenfalls auf meinen 
Rechner. und dann die ersten 10 Indizes vergleichen

also bei dem ersten 10*10 Werte habe ich als Ergebnis verglichen und die 
Laufzeit gemessen.

1

2

3

for (int i = 0; i < MU; i++)

4

    {

5

        for (int j = 0; j < MU; j++)

6

        {

7

            A[i][j] = (float)i * (float)j;

8

            B[i][j] = (float)i * (float)j;

9

        }

10

    }

11

12

13

    tstart = clock();

14

    StandardMultTansposeBlocked(A, B, D);

15

    time1 += clock() - tstart;     // end

16

    time1 = time1 / CLOCKS_PER_SEC;  // rescale to seconds

17

    cout << "  time transpose blocked = " << time1 << " sec." << endl;

18

19

20

float** StandardMultTansposeBlocked(float **M1, float **M2,float** result)

21

{

22

    int i, j, k;

23

    float sum=0;

24

    for (i = 0; i < MU; i++)

25

        for (j = 0; j < MU; j++)

26

        {

27

           M2[i][j] = M2[j][i];

28

           result[i][j] = 0;

29

        }

30

31

    int loop = 1;

32

33

    for (int count = 0; count <= MU; count += loop)

34

    {

35

        for (i = count; i < loop + count; i++) {

36

37

38

            for (j = 0; j < MU; j++) {

39

40

               float* M1ik = &M1[i][0];

41

                float* M2jk = &M2[j][0];

42

                float sum;

43

                sum = result[i][j];                 

44

                for (k = 0; k < MU; k++) {

45

                   sum += M1ik[k] * M2jk[k];

46

                }

47

                result[i][j] = sum;

48

            }

49

50

        }

51

    }

ich denke die innerste Schleife kann man noch mehr optimieren. 
Vielleicht muss man das aufspalten. Naja nur noch Faktor 10 hinter 
Matlab. So schlecht ist das nicht. Vorher Faktor 30. In Summe macht das 
schon einen großen Unterschied, wenn man 100mal 2 große Matrizen 
multiplizieren muss oder so. Nochmal Faktor 2 Verbesserung und es würde 
für meine Anwendung schon reichen und könnte mit eine externe Lib 
ersparen :D

Eric schrieb:
> hab die Matrizen gefüllt mit A[i][j]=(float)i*(float)j Index.

Die Matrizen sind nicht ausreichend, um die Funktion zu testen, da sie 
symmetrisch ist. Transponieren der M2 Matrix sollte also egal sein. Du 
solltest mindestens i * N + j nehmen, wobei N sowas wie max(j) + 1  oder 
pow(10, floor(log10(max(j))) + 1) damit man jederzeit dem wert ansehen 
kann, wo er in der Matrix steht. M1 und M2 sollten sich auch 
unterscheiden, damit du ein versehentliches vertauschen der Matrizen 
erkennen kannst.

Und wie testes du, ob das Ergebnis der Multiplikation richtig ist?

Walter T. schrieb:
> Ich kenne das Forum. Wenn man behauptet, dass Speicher nichts mit Zeilen
> und Spalten zu tun hat, kommt irgendjemand aus einem Loch gekrochen und
> zieht freudig einen Kernspeicher hervor.

Du bist also gerade aus einem Loch gekrochen? Du hast den Kernspeicher 
hervorgezogen ... oder war das nicht freudig?

Walter T. schrieb:
> mh schrieb:
>> Was ist der Inhalt und wie testest du das Ergebnis?
>
> Die Frage nach den Testfällen ist gut. Isnbesondere, wenn sehr große und
> sehr kleine Werte multipliziert-addiert werden.
Es ist erstmal wichtig für den grundlegenden Test, ob seine Funktion 
überhaupt funktioniert.

Diese For Schleife ist richtig, bei der anderen probiere ich gerade noch 
rum

[c]   for (i = 0; i < MU; i++) {
        for (j = 0; j < MU; j++) {
            float* M1ik = &M1[i][0];
            float* M2jk = &M2[j][0];
            sum = 0;
            for (k = 0; k < MU; k++) {
                sum += M1ik[k] * M2jk[k];
            }
            result[i][j] = sum;
        }
    }

[/c}

selbst diese ist schon Faktor 3 besser, das Ergebnis stimmt. Testen mach 
ich irgendwann. Die Matrix muss vorher transponiert werden, sonst stimmt 
das Ergebnis nicht.

Grüße

Andreas M. schrieb:
> Dieser Code transponiert die Matrix aber nicht! Das wurde schon
> mindestens zwei mal geschrieben. Denke daran das Du hier auf der selben
> Matrix schreibst, die Du gerade liest.
>
> Eric schrieb:
>> Es soll später auch auf einem ARM
>> Processor laufen und.
>
> Eric schrieb:
>> Ich will
>> auch nicht besser sein als die fertigen libs, allerdings schneller als
>> 0,5s.
>
> Wenn Dein Code auf einem ARM in unter 0.5s laufen soll, dann ist es
> komplett sinnfrei es auf einer Intelmaschine auf unter 0.5s zu
> optimieren. Du musst es auf dem Zielsystem machen, ARMs funktionieren
> komplett anders als Intel CPUs

klar meine Testmatrix wird transponiert und ich kann in der innersten 
Schleife durch die spalten gehen. Weiss nicht was du meinst. Ich sehe du 
hast das noch nie ausprobiert. Das ist nur Testcode mit dem ich 
rumspiele. Eine Quadratische Matrix wird so transponiert oder was ist 
daran falls.

Auf einem Arm in abgespeckter Form Matritzen kleiner. Step by Step

Mir ist egal ob die Matrix jetzt richtig tansponiert wird spielt für 
meinen Test keine Rolle

Eric schrieb:
> Step by Step

Dann nimm dir mal einen Blatt Papier und einen spitzen Bleistift und 
gehe deinen Algorithmus, der diese Zeilen enthält
Eric schrieb:
>     for (i = 0; i < MU; i++)
>         for (j = 0; j < MU; j++)
>            M2[i][j] = M2[j][i];
Schritt für Schrit durch für die folgenden Matrizen:

1

M1 = {

2

  {110, 111, 112},

3

  {120, 121, 122},

4

  {130, 131, 132},

5

}

6

M2 = {

7

  {210, 211, 212},

8

  {220, 221, 222},

9

  {230, 231, 232},

10

}

ich vergleiche Standard Matrix Multiplikation mit einer "verbesserten" 
Standard Multiplikation.

Eingangsadaten sind gleich, es wird bei beiden identisch Transponiert, 
auch wenn es nicht ganz richtig ist. Ergo kommt auch das gleiche (wenn 
auch falsche ) raus. Ich habe für meinen Test nix gewonnen, selbst wenn 
ich es richtig implementiere. Mich interssiert nur Laufzeit.

Eric schrieb:
> Mich interssiert nur Laufzeit.
Das trifft es wohl ziemlich genau. An konstruktiver Kritik scheinst du 
zumindest nicht interessiert zu sein.

Das war nicht Böse gemeint, aber es ist tatsächlich für meinen Mini-Test 
egal.
Aber keine Angst, bei meinem richtigen Code wird es auf eine zweite 
Matrix geschrieben.

Grüße, naja ich höre jetzt eh erstmal auf. Ich habe Urlaub und mach 
schon wieder zu viel für die Arbeit.

Der Grund, warum Matlab bei vielen Operationen schneller ist als eigene 
C-Implementierungen ist kein Verdienst von Mathworks.
Wenn Du einen Intel-Prozessor hast, wird die Intel-Math-Kernel-Library 
(MKL) geladen, die hoch optimierte BLA-Funktionen enthält. Die MKL macht 
Gebrauch von den  AVX-Funktionen, die der vorhandene Prozessor bietet 
und unterstützt Multicore.
Du  kannst bei Intel ein Package mit der MKL runterladen und diese auch 
aus deinen eigenen Programmen aufrufen.

Ja da werd ich wohl das mkl mal ausprobieren. Schein ja auch für AMD zu 
gehen. Lizenz braucht man dafür ja nicht. Und dann inkludiert meine 
Matrizen Library eben das MKL. Ich denke so mach ich es.

Jo, Grüße

1

A_SUB = (double*)mkl_malloc(m * k * sizeof(double), 64);

2

    B_SUB = (double*)mkl_malloc(k * n * sizeof(double), 64);

3

    C_SUB = (double*)mkl_malloc(m * n * sizeof(double), 64);

4

5

    count = 0;

6

    for (int i = 0; i < m; i++)

7

    {

8

        for (int j=0;j<k;j++ )

9

        {

10

            A_SUB[count] = B[i][j];

11

            B_SUB[count] = B[i][j];

12

            count++;

13

        }

14

    }

15

16

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,  m, n, k, alpha, A_SUB, k, B_SUB, n, beta, C_SUB, n);

okay bin bei 0,16 und damit 10 mal schneller als mit Standard Mult. 
Nutze jetzt MKL auf einem AMD Rechner (glaube auf einen Intel sollte es 
noch schneller laufen). Allerdings geht das parallele Zeug bei mir 
nicht, dgemm geht nur sequentiell. Stell ich auf parallel um bekomm ich 
eine Exception Kernelbase.ddl bzw. thread.1dll  Einstiegspunkt findet er 
nicht. Allerdings beschädigt sind die Dinger denke ich nicht. Habe 
gerade ein Windows Update durchgeführt. Hatte einer auch schon mal das 
Problem mit der KernelBase ?

ah hab es herausegefuncen, parallel geht nur wenn man den Threading 
Building Block aktiviert  :D. Naja Doku hilft :D

dann ist man schon 20 mal schneller bei einer 2000*2000 Matrix :D. Jetzt 
nur noch mit diesen WOrkaround für AMD Prozessoren beschäftigen und dann 
bin ich voll zufrieden. Obwohl die Performance jetzt schon ausreicht :D

Ja das ist mein Testcode, m,k,n sind bei mir zum Test alle gleich . 
Eingangsdaten lasse ich aber als 2D Array. Ist einfacher anzuwenden - 
selbst wenn es paar us länger dauert als ein 1D Array. Obwohl es bei 
Standard Mult bei mir keinen Unterschied gemacht hat 9b 1D oder 2D. Ob 
memcpy schneller ist probiere ich aus. Gibt aber auch ein mkl cpy, 
glaube. Bin noch nicht ganz durch auf der IBM Seite, nur überflogen.

Für Matrizen bis 25×25 scheint auch Standard mult besser. Bei parallel 
ist dieser super Code erst schneller für matrizen grösser 100. Zumindest 
auf meinem Prozessor.

Eine Anekdote zu Mathe Bibliothek von ARM:

Ich hatte mal eine template basierte Matrix mal Vektor Multiplikation 
geschrieben und die Performance gegen die generische Matrix mal Matrix 
Multiplikation in der ARM Bibliothek verglichen (zu berechnende Daten 
waren 256x256 Matrix mal 256 Vektor).

Die ARM Variante war schneller. Aber nur, weil der GCC zu schlecht ist 
meine Schleifen auszurollen. ARM macht dort manuelles Ausrollen. Mit 
manuellen Ausrollen war dann auch mein Code schneller (was zu erwarten 
war, da nicht generisch).

arm_math schrieb:
> Mit
> manuellen Ausrollen war dann auch mein Code schneller (was zu erwarten
> war, da nicht generisch).

Spannend. Mit etwas template-Magic kann man Schleifen-Ausrollen auch 
generisch erzwingen.