Forum: PC-Programmierung C++ Amp Matrix Multiplikation vs MatLab


von ERic (Gast)


Lesenswert?

Hallo, hab mir eine Funktion geschrieben um die Eigenwerte einer Matrix 
zu berechnen.

Funktioniert soweit ganz gut nur die Rechenzeit macht Probleme.
Der Flaschenhals sind die Matrixmultiplikationen. Mache die Eigenwert 
berechnung über Householdermatrizen, Matlab nuzt glaube 
Givensrotationen, was eigentlich etwas länger dauern müsste.
1
Library_Return_Values MatricesOperations::MultParallel(double **matrixA, unsigned int sizeRowA, unsigned int sizeColumnA, double **matrixB, unsigned int sizeRowB, unsigned int sizeColumnB, double **matrixOut, unsigned int*sizeRowOut, unsigned int *sizeColumnOut)
2
  {
3
    unsigned int k = 0,i,j;
4
    
5
6
    if (sizeRowA*sizeColumnB < 250000)
7
    {
8
      return MultTrans(matrixA, sizeRowA, sizeColumnA, matrixB, sizeRowB, sizeColumnB, matrixOut, sizeRowOut, sizeColumnOut);
9
    }
10
11
    accelerator usedAccelerator=GetBestAccelerator();
12
13
    for (i = 0; i < sizeRowA; i++)
14
    {
15
      for (j = 0; j < sizeColumnA; j++)
16
      {
17
        GLB_A[k] = matrixA[i][j];
18
        k++;
19
      }
20
    }
21
22
    k = 0;
23
    for (i = 0; i < sizeRowB; i++)
24
    {
25
      for (j = 0; j < sizeColumnB; j++)
26
      {
27
        GLB_B[k] = matrixB[i][j];
28
        k++;
29
      }
30
    }
31
32
    array_view<double, 2> a(sizeRowA, sizeColumnA, GLB_A);
33
    array_view<double, 2> b(sizeRowB, sizeColumnB, GLB_B);
34
    array_view<double, 2> c(sizeRowA, sizeColumnB, GLB_C); c.discard_data();
35
36
    concurrency::parallel_for_each(usedAccelerator.default_view, c.extent,[=](concurrency::index<2> idx) restrict(amp) 
37
    {
38
      int row = idx[0]; int col = idx[1];
39
      double sum = 0.0f;
40
      for(int i = 0; i < b.extent[0]; i++)
41
      {
42
        sum += a(row, i) * b(i, col);
43
        c[idx] = sum;
44
      }
45
    });
46
    c.synchronize();
47
        
48
    k = 0;
49
    for (i = 0; i < sizeRowA; i++)
50
    {
51
      for (j = 0; j < sizeColumnB; j++)
52
      {
53
        matrixOut[i][j] = GLB_C[k];
54
        k++;
55
      }
56
    }
57
58
    *sizeColumnOut = sizeColumnB;
59
    *sizeRowOut = sizeRowA;    
60
    return RET_OK;
61
  }



Ich habe eine Zeitmessung gemacht. Matrix Größen 1000 C=A*B
Matrix Normale Berechnung: 4,2s
Matrix Normal mit Transpose: 1,7s
Matrix parallel: 0,5s.

Diese Ergebnisse habe ich mit MAtLab verglichen. MatLab braucht nur 
0,07s für eine Berechnung von zwei 1000er Berechnungen. Warum ist meine 
Multiplikation so langsam. Wie bekommt man das schneller hin ?
Im Internet schreiben die (MSDN) immer das so eine Multiplikation 
parallel um Faktor 100 schneller sein soll. Als Beschleuniger nehme ich 
den mit dem größten dedizierten Speicher. Die anderen waren noch 
schlechter. Gibt es da irgendwie einen Trick ?

Grüße Eric

von Programmierer (Gast)


Lesenswert?


von abc (Gast)


Lesenswert?

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

von ERic (Gast)


Lesenswert?

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

von ERic (Gast)


Lesenswert?

Ich meine dieses AVX

von ERic (Gast)


Lesenswert?

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

von Programmierer (Gast)


Lesenswert?

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

von abc (Gast)


Lesenswert?

Du kannst dir aber anschauen wie Eigen das implementiert

von M.K. B. (mkbit)


Lesenswert?

Auf welchem System testet du denn? Evtl. hast du auch Cachemisses, die 
bremsen.

Wenn z.B. deine Matrix ein Array aus Pointern ist und relativ groß ist, 
dann kannst du bei jedem Zeilen- bzw. Spaltenwechsel eine Wartezeit 
haben. Das hängt aber von der konkreten Struktur deiner Daten ab.

Evtl. gibt es Profiler für deine CPU, mit denen man sowas analysieren 
kann.

von Irgend W. (Firma: egal) (irgendwer)


Lesenswert?

ERic schrieb:
> Möchte keine andere Bibliothek nehmen, will ja was lernen :D
Dann wirst du dir wohl deine low level Bibliothek selber schreiben 
müssen und dort die passenden X86-Befehle mit inline Assembler reinbauen 
müssen.

- https://de.wikipedia.org/wiki/Advanced_Vector_Extensions
- http://0x80.pl/notesen/2019-02-02-autovectorization-gcc-clang.html

.
- Zuvor aber erstmal den Assemblercode von deiner jetzigen Version 
anschauen ob da eventuell schon was verwendet wird.
- Oder auch mal nachschauen ob es irgendwo einen Schalter gibt damit es 
überhaupt verwendet wird.
- Mal nachschauen ob überhaupt die richtige Prozessor Generation 
eingestellt ist, nicht da das was Uraltes angenommen wird das das noch 
garnicht unterstützt.

von Walter T. (nicolas)


Lesenswert?

ERic schrieb:
> Matrix Größen 1000 C=A*B

Also sparse? Das wird wirklich knackig, die seit Jahrzehnten optimierten 
Libraries zu knacken.

von Andreas M. (amesser)


Lesenswert?

Wieso Pointer-Pointer für die Matrizen? Das ist unnötig und versaut die 
Lokalität. Dazu dann das unnötig umkopieren. Parallelisieren wird man 
auch eher mehrere Multiplikationen als innerhalb der gleichen 
Multiplikation. Das so Aufzuteilen ist Maximal ineffektiv, da können die 
L1 Caches nicht vernünftig arbeiten, weil jeder deiner "parallelen" 
Tasks auf einem zufälligen Core landet und damit potentiell die falschen 
Daten im Cache liegen. Wenn man es sequentiell macht, dann liegen die 
Daten der letzten Zeile/Spalte noch im L1 und können direkt wieder 
verwendet werden, wenn es aber ein anderer Core ist....

Ü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)

: Bearbeitet durch User
von Eric (Gast)


Lesenswert?

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.

von Programmierer (Gast)


Lesenswert?

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

von Walter T. (nicolas)


Lesenswert?

Walter T. schrieb:
> Also sparse? Das wird wirklich knackig, die seit Jahrzehnten optimierten
> Libraries zu knacken.

Das war natürlich Unsinn. Das geht ja noch als volle Matrix knapp in den 
L3-Cache.

Programmierer schrieb:
> ARM hat dazu noch
> verschiedene Vector Extensions, NEON ist nur eine davon.

Ich hoffe, wir reden nicht von den kleinen Cortex M. Die ganze 
FPU-Geschichte läuft bei denen ja in einfacher Genauigkeit. Die 
Standard-Verfahren aus der Numerik-Literatur hauchen da aber bei 
mittelmäßigen Konditionszahlen schon ihren Nutzen aus. Also entweder 
doppelte Genauigkeit und auf die Extensions verzichten, oder 
Vorkonditionierung betreiben.

von abc (Gast)


Lesenswert?

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.

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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?

von Eric (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

[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]

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

oh ja ist noch ein Bug drin :D Sorry

von Eric (Gast)


Lesenswert?

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
    }

von Eric (Gast)


Lesenswert?

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

von Walter T. (nicolas)


Lesenswert?

Eric schrieb:
> Bei größeren Matrizen ist es sogar Faktor 4 schneller

Das klingt für eine derart naive Implementierung erst einmal merkwürdig. 
Aber OK, bei einer Standard-Implementierung ist die Transponierung schon 
mit drin, MU ist eine Variable und keine Compilezeitkonstante, 
rechteckige Matrizen funktionieren auch und meistens ist ja auch double.

Wie sieht die Geschwindigkeit aus, wenn MU nicht mehr konstant ist?

: Bearbeitet durch User
von Eric (Gast)


Lesenswert?

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
    }

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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.

von Walter T. (nicolas)


Lesenswert?

mh schrieb:
> Vorsicht! Transponieren ist eine mathematische Operation.

Naja, in Tensorschreibweise ist es ein Tauschen der Zeilen- und 
Spaltenindizes. Das passt hier schon.

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

Und ist der richtige Zeitpunkt, schon seit vorgestern automatische 
Modultests zu haben. Zum Glück unterstützt Matlab das endlich (seit 
2014) auch.

: Bearbeitet durch User
von Eric (Gast)


Lesenswert?

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.

von mh (Gast)


Lesenswert?

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.

von Walter T. (nicolas)


Lesenswert?

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.

Was mich erstaunt ist, dass eine doppelt indizierte Schreibweise schnell 
sein soll. Als ich mich damals mit dem Thema beschäftigt habe, war eine 
der ersten Massnahmen zur Beschleunigung von Skalarprodukten auf eine 
einfache Indizierung zu wechseln, damit sich der Optimizer mit loop 
unrollig austoben kann und nicht so viel Zeit mit der Prüfung der 
Schleifenbedingung verbringt.

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

: Bearbeitet durch User
von mh (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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.

von Eric (Gast)


Lesenswert?

Testmatrizen waren zwei 1000x1000 Matrizen.

von mh (Gast)


Lesenswert?

Eric schrieb:
> Testmatrizen waren zwei 1000x1000 Matrizen.

Was ist der Inhalt und wie testest du das Ergebnis?

von Eric (Gast)


Lesenswert?

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

von Walter T. (nicolas)


Lesenswert?

mh schrieb:
> mit
> Speicherlayout meine ich die Anordnung der Daten im Speicher, also sowas
> wie Row-Major oder Column-Major

Ich auch. Und ich habe keine Skrupel, die Umsortierung zwischen dem 
einen und dem anderen als "transponieren" zu bezeichnen.

mh schrieb:
> Was hat ein Kernspeicherentwickler mit
> Matrizenmultiplikation zu tun?

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.

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.

von Eric (Gast)


Lesenswert?

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
    }

von Eric (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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.

von Eric (Gast)


Lesenswert?

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}

von Andreas M. (amesser)


Lesenswert?

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

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

Eric schrieb:
> Ja ist ein Cortex m, hat aber keine Floating Point Unit. Na gut, jetzt
> aber Urlaub.

LOL! Cortex-M. Drei 1000x1000 Matrizen benötigen je nach verwendetem 
Datentyp (8Bit,16Bit,32Bit/float,64Bit/double) etwa 
2,8MB/5,8MB/11,5MB/23MB RAM. Was für ein Chip mit Cortex-M ist das denn, 
der soviel RAM hat?

Externer RAM? Na dann: Matrix-Mul mit 1000er Matrixen heißt 1Mrd 
Multiplikationen = 2Mrd Loads. Wenn die Multiplikation nur 0.5s dauern 
soll heißt dass, die Daten müssen mit ca. 30GB/s in die ALU gehen. 
Deswegen nochmal die Frage: Was ist das Zielsystem?

von Walter T. (nicolas)


Lesenswert?

Eric schrieb:
> for (int i = 0; i < MU; i++)
>     {
>         for (int j = 0; j < MU; j++)
>         {
>             A[i][j] = (float)i * (float)j;
>             B[i][j] = (float)i * (float)j;
>         }
>     }

Dir ist aber klar, dass Zahlen zwischen 0 und 1000 so ziemlich die 
harmloseste Variante sind, die float zu bieten hat? Interessant ist es, 
wenn auch Zahlen in der Größenordnung 100*eps oder 1e30 dabei sind.

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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

von mh (Gast)


Lesenswert?

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
}

von Eric (Gast)


Lesenswert?

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.

von mh (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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.

von Unwichtig (Gast)


Lesenswert?

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.

von Andreas M. (amesser)


Lesenswert?

Und genau wie Intel die MKL bereitstellt, gibt es das auch von ARM 
fertig optimiert in Form der CMSIS DSP Bibliothek:

https://arm-software.github.io/CMSIS_5/DSP/html/group__groupMatrix.html

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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 ?

von Eric (Gast)


Lesenswert?

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

von Eric (Gast)


Lesenswert?

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

von Klaus W. (mfgkw)


Lesenswert?

Zu deinem Umkopieren von B nach A_SUB bzw. B_SUB:

Du scheinst deine 2D-Felder B etc. als Feld von double* organisiert zu 
haben, wobei die double* dann auf je ein Feld mit den doubles einer 
Zeile zeigen.
cblas_dgemm() dagegen erwartet wohl ein Feld von double, in dem 
hintereinander alle Werte stehen.

1. Wenn du dein Programm so umbaust, daß du dieselbe Organisation hast, 
könntest du entweder jeweils ein ganzes Feld mit memcpy() kopieren oder 
(falls du nicht wirklich eine Kopie brauchst) gleich dein Feld übergeben 
und sogar auf die Kopie verzichten.

2. Selbst wenn du die Organisation deiner Feld lässt als Feld von 
double*..., könntest du zumindest jeweils eine ganze Zeile einer Matrix 
kopieren statt der k Einzelwerte in der inneren Schleife.
Also die innere Schleife ersetzen durch einen memcpy()-Aufruf.

Möglicherweise ist memcpy() schneller als das einzelne Kopieren?

---

Was mir aber noch komisch vorkommt an deinem letzten Quelltext:
A_SUB ist m*k groß, B_SUB aber k*n.
Beide werden aber mit m*k werten gefüllt.
Ist das so richtig?

: Bearbeitet durch User
von Eric (Gast)


Lesenswert?

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.

von arm_math (Gast)


Lesenswert?

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

von Programmierer (Gast)


Lesenswert?

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.

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.