Matrizen / Vectoren - Multiplikationen beschleunigen

Zur Vorstellung von Komponenten und Units für Lazarus
Antworten
Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Wen man Vektoren / Matrizen mutiplizieren will, kann man die um fas das 10 fache beschleunigen, wen man die SSE-Einheit der CPU benutzt.
Das ganze ist ideal für OpenGL-Freunde. Es läuft mit 32 und 64Bit Intel/AMD-CPUs.

Code: Alles auswählen

Type
  TVector4f = array[0..3] of Single;
  TMatrix = array [0..3, 0..3] of Single


Klassisch mit Single:

Code: Alles auswählen

operator * (const mat0, mat1: TMatrix) Res: TMatrix;
var
  i, j, k: integer;
begin
  for i := 0 to 3 do begin
    for j := 0 to 3 do begin
      Res[i, j] := 0.0;
      for k := 0 to 3 do begin
        Res[i, j] += mat1[i, k] * mat0[k, j];
      end;
    end;
  end;
end


SSE-Beschleunigt:

Code: Alles auswählen

function VectorMultiplySSE(const mat: TMatrix; const vec: TVector4f): TVector4f; assembler;
  {$asmmode intel}
asm
         // pointer to the first 4 elements of the array:
         // Move Unaligned Parallel Scalars
         // load the registers with matrix values:
         Movups  Xmm4, [mat+$30] // xmm4 contains 1,5, 9,13
         Movups  Xmm5, [mat+$20] // +16 (4 bytes x 4 elements)
         // xmm5 contains 2,6,10,14
         Movups  Xmm6, [mat+$10] // +32 (8 bytes x 4 elements)
         // xmm6 contains 3,7,11,15
         Movups  Xmm7, [mat+$00] // +48 (12 bytes x 4
         // esi contains the starting address of the vec array
         // [2, 3, 4, 5], so load this vector into xmm0:
         Movups  Xmm0, [vec] // xmm0 contains 2, 3, 4, 5
         // we'll store the final result in xmm2, initialize it to 0:
         Xorps   Xmm2, Xmm2 // xmm2 contains 4x32
         // bits with zeros
         // now we need to multiply first column (xmm4)
         // by the vector (xmm0) and add it to the total (xmm2)
         // copy content of xmm0 into xmm1:
         Movups  Xmm1, Xmm0 // xmm1 now contains
         // [2, 3, 4, 5]
         // each value in xmm1 has the following mask representation:
         // mask value: 11 10 01 00
         // register value: [ 2, 3, 4, 5 ]
         // Shuffle Parallel Scalars
         Shufps  Xmm1, Xmm1, $FF // FF mask is 11 11 11 11
         // xmm1 contains 2, 2, 2, 2
         // Multiply Parallel Scalars
         Mulps   Xmm1, Xmm4 // multiply xmm1 (2,2,2,2)
         // and xmm4 (1,5,9,13)
         // [ 2*1, 2*5, 2*9, 2*13 ]
         // save it in xmm1
         // Add Parallel Scalars
         Addps   Xmm2, Xmm1 // add xmm2 (0, 0, 0, 0)
         // and xmm1 (2, 10, 18, 26)
         // save it in xmm2
         // xmm2 contains [2,10,18,26]
         // we multiplied first column of the matrix by the vector,
         // now we need to repeat these operations for the remaining
         // columns of the matrix
         Movups  Xmm1, Xmm0 // 3, 3, 3, 3
         Shufps  Xmm1, Xmm1, $AA // AA -> 10 10 10 10
         Mulps   Xmm1, Xmm5 // xmm5: 2, 6, 10, 14
         Addps   Xmm2, Xmm1 // 2+6, 10+18, 18+30, 26+42
         Movups  Xmm1, Xmm0 // 4, 4, 4, 4
         Shufps  Xmm1, Xmm1, $55 // 55 -> 01 01 01 01
         Mulps   Xmm1, Xmm6 // xmm6: 3, 7, 11, 15
         Addps   Xmm2, Xmm1 // 8+12, 28+28, 48+44, 68+60
         Movups  Xmm1, Xmm0 // 5, 5, 5, 5
         Shufps  Xmm1, Xmm1, $00 // 00 -> 00 00 00 00
         Mulps   Xmm1, Xmm7 // xmm7: 4, 8, 12, 16
         Addps   Xmm2, Xmm1 // 20+20, 56+40, 92+60, 128+80
 
         // 40 , 96 , 152 , 208
         // write the results to vectorOut
         Movups  [result], Xmm2
end;
 
 
function MatrixMultiplySSE(const M1, M2: TMatrix): TMatrix; assembler;
asm
         // load M1 into xxm entirely as we will need it more than once
         Movups   Xmm4, [M1+$00] // movaps
         Movups   Xmm5, [M1+$10]
         Movups   Xmm6, [M1+$20]
         Movups   Xmm7, [M1+$30]
         // compute Result[0]
         Movss    Xmm0, [M2]
         Shufps   Xmm0, Xmm0, $00 //xmm0 = 4x M2[0,0]
         Mulps    Xmm0, Xmm4
         Movss    Xmm1, [M2+$04]
         Shufps   Xmm1, Xmm1, $00 //xmm1 = 4x M2[0,1]
         Mulps    Xmm1, Xmm5
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$08]
         Shufps   Xmm1, Xmm1, $00 //xmm1 = 4x M2[0,2]
         Mulps    Xmm1, Xmm6
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$0c]
         Shufps   Xmm1, Xmm1, $00 //xmm1 = 4x M2[0,3]
         Mulps    Xmm1, Xmm7
         Addps    Xmm0, Xmm1
         Movups   [Result+$00], Xmm0 // movntps
         // compute Result[1]
         Movss    Xmm0, [M2+$10]
         Shufps   Xmm0, Xmm0, $00
         Mulps    Xmm0, Xmm4
         Movss    Xmm1, [M2+$14]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm5
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$18]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm6
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$1c]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm7
         Addps    Xmm0, Xmm1
         Movups   [Result+$10], Xmm0
         // compute Result[2]
         Movss    Xmm0, [M2+$20]
         Shufps   Xmm0, Xmm0, $00
         Mulps    Xmm0, Xmm4
         Movss    Xmm1, [M2+$24]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm5
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$28]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm6
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$2c]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm7
         Addps    Xmm0, Xmm1
         Movups   [Result+$20], Xmm0
         // compute Result[3]
         Movss    Xmm0, [M2+$30]
         Shufps   Xmm0, Xmm0, $00
         Mulps    Xmm0, Xmm4
         Movss    Xmm1, [M2+$34]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm5
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$38]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm6
         Addps    Xmm0, Xmm1
         Movss    Xmm1, [M2+$3c]
         Shufps   Xmm1, Xmm1, $00
         Mulps    Xmm1, Xmm7
         Addps    Xmm0, Xmm1
         Movups   [Result+$30], Xmm0
end;


Wen man dieses Test-Programm laufen lässt, stell man fas eine 10fache Beschleunigung fest.

Code: Alles auswählen

procedure TForm1.Button1Click(Sender: TObject);
var
  i: integer;
  m, m0, m1: Tmat4x4;
  t: TTime;
const
  site = 20000001;
 
  procedure Ausgabe;
  var
    i: integer;
  begin
    for i := 0 to 3 do begin
      WriteLn(m0[i, 0]: 4: 2, '  ', m0[i, 1]: 4: 2, '  ', m0[i, 2]: 4: 2, '  ', m0[i, 3]: 4: 2);
    end;
    WriteLn();
  end;
 
  function GetZeit(z: TTime): string;
  begin
    str(z * 24 * 60 * 60: 10: 4, Result);
  end;
 
begin
  m0.Identity;
  m1.Identity;
  m1.RotateC(3.5);
 
  t := now;
  for i := 0 to site do begin
    m0 := m0 * m1;
  end;
  WriteLn('FPU:   ', GetZeit(now - t));
  Ausgabe;
 
 
  m0.Identity;
  m1.Identity;
  m1.RotateC(3.5);
 
  t := now;
  for i := 0 to site do begin
    m0 := MatrixMultiplySSE(m0, m1);
  end;
  WriteLn('SSE:    ', GetZeit(now - t));
  Ausgabe;
 
  WriteLn();
  WriteLn();
  WriteLn();
end;


Die Ausgabe:

Code: Alles auswählen

FPU:       2.6060
0.55  0.69  0.00  0.00
-0.69  0.55  0.00  0.00
0.00  0.00  1.00  0.00
0.00  0.00  0.00  1.00
 
SSE:        0.2920
0.55  0.69  0.00  0.00
-0.69  0.55  0.00  0.00
0.00  0.00  1.00  0.00
0.00  0.00  0.00  1.00
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Frank Ranis
Beiträge: 201
Registriert: Do 24. Jan 2013, 21:22

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Frank Ranis »

Hallo Mathias ,

das klingt ja wie Zauberei .

In meinen Aerodynamikprogrammen steckt jede Menge Matrixrechnung drinn.
Da geht es darum lineare Gleichungssysteme (Einfluß_Matrix * Vektor_unbekannt = RHS) zu lösen , das mach ich in der Regel mit dem Gauss Eliminations-Verfahren.

Die Matritzen und Vektoren habe aber ganz andere Dimensionen , so in der Größe von 500*500 bis > 2000*2000 , oder noch mehr.

Jetzt hast Du ja in deinem Demo-Code nun ne 4*4-Matrix * Vektor -Multipikation , quasi per Hand runterprogrammiert .

Siehst Du es eine Möglichkeit die Dimensionen beliebig auszudehnen ?
Und auf die Art nen kompletten Gauss-Löser zu erstellen ? , das wäre der Brüller .

Als Input hätte man dann die Einflußmatrix und den Zielvektor (RHS=rechte Seite) und als Result kämen der aufgelöste (Vektor_unbekannt) zurück , also die Lösung.

Als Abklatsch von dem Gauss könnte man dann auch z.B. ein Matrix-Invert bauen .
Mit einer invertierten Matrix * RHS kommt dann ja auch die Lösung heraus .

Gruß

Frank
www.flz-vortex.de

Horst_h
Beiträge: 72
Registriert: Mi 20. Mär 2013, 08:57

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Horst_h »

Hallo,

nimmt man für solche Dimensionen nicht BLAS/BLIS oder LAPACK
http://forum.lazarus.freepascal.org/ind ... c=27730.45 und folgend

Gruß Horst

mischi
Beiträge: 206
Registriert: Di 10. Nov 2009, 18:49
OS, Lazarus, FPC: macOS, 10.13, lazarus 1.8.x, fpc 3.0.x
CPU-Target: 32Bit/64bit

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von mischi »

Wieso nicht gleich BLAS/Lapack benutzen? Mit etwas Glück ist das sogar noch multithreading. Zumindest bei der Benutzung von Python für Simulationen ist das der Weg.
MiSchi macht die fink-Pakete

Frank Ranis
Beiträge: 201
Registriert: Do 24. Jan 2013, 21:22

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Frank Ranis »

Hi,

Horst_h hat geschrieben:Hallo,

nimmt man für solche Dimensionen nicht BLAS/BLIS oder LAPACK
http://forum.lazarus.freepascal.org/ind ... c=27730.45 und folgend

Gruß Horst

mischi hat geschrieben:Wieso nicht gleich BLAS/Lapack benutzen? Mit etwas Glück ist das sogar noch multithreading. Zumindest bei der Benutzung von Python für Simulationen ist das der Weg.


BLAS/Lapack , äh , ja ,

gibt es dazu ein paar lauffähige LAZ-Demos ? , Plattformunabhängig (Ich frage wegen den DLL's die da zu saugen sind ) .

Gruß

Frank
www.flz-vortex.de

Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Jetzt hast Du ja in deinem Demo-Code nun ne 4*4-Matrix * Vektor -Multipikation , quasi per Hand runterprogrammiert .

SSE ist für 4x4-Matrizen ausgelegt, welche man für OpenGL braucht. Die Register sind auch 128Byte gross, als 4 Single im Vektor.

Vielleicht kannst es trotzdem verwenden, wen man es zusammen bastelt.
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Warf
Beiträge: 1907
Registriert: Di 23. Sep 2014, 17:46
OS, Lazarus, FPC: Win10 | Linux
CPU-Target: x86_64

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Warf »

Horst_h hat geschrieben:BLAS/Lapack , äh , ja ,

gibt es dazu ein paar lauffähige LAZ-Demos ? , Plattformunabhängig (Ich frage wegen den DLL's die da zu saugen sind ) .

Gruß

Frank


Lazarus Demos zu finden könnte eventuell etwas schwerer werden, aber es sollte auf jeden fall genug C beispiele geben, und die kann man grundsätzlich recht einfach übernehmen.

Aber wenn man es wirklich schnell haben will macht man den Spaß sowiso nicht auf der CPU sondern auf der GPU. cuBLAS von Nvidia läuft selbst auf schlechten GPU's immernoch doppelt so schnell wie auf einer modernen CPU.

mischi
Beiträge: 206
Registriert: Di 10. Nov 2009, 18:49
OS, Lazarus, FPC: macOS, 10.13, lazarus 1.8.x, fpc 3.0.x
CPU-Target: 32Bit/64bit

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von mischi »

Frank Ranis hat geschrieben:gibt es dazu ein paar lauffähige LAZ-Demos ? , Plattformunabhängig (Ich frage wegen den DLL's die da zu saugen sind ) .

Auf macOS sind BLAS/Lapack immer mit dabei, bei Linux weiss ich es nicht so genau, vermute aber auch.
MiSchi macht die fink-Pakete

Horst_h
Beiträge: 72
Registriert: Mi 20. Mär 2013, 08:57

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Horst_h »

Hallo,

http://forum.lazarus.freepascal.org/ind ... #msg172277
dort wird auch eine sourceforge adresse angegeben, aber die neueste Version mit Windows DLL ist wohl.
https://sourceforge.net/projects/openblas/files/v0.2.19
Ich habe jetzt v0.3 compiliert für Linux, aber der Aufruf hat sich geändert, da taucht ein blas_arg_t auf, da sind wohl 4 pdouble statt 3 , kann ich jetzt nicht finden
[c]
 
//common_level3.h:619: int dgemm_thread_nn(blas_arg_t *, BLASLONG *, BLASLONG *, double *, double *, BLASLONG);
typedef struct {
void *a, *b, *c, *d, *alpha, *beta;
BLASLONG m, n, k, lda, ldb, ldc, ldd;
 
#ifdef SMP
void *common;
BLASLONG nthreads;
#endif
 
#ifdef PARAMTEST
BLASLONG gemm_p, gemm_q, gemm_r;
#endif
 
#ifdef PREFETCHTEST
BLASLONG prea, preb, prec, pred;
#endif
 
} blas_arg_t;
 
[/c]
dann passt das nicht mehr :-(

Code: Alles auswählen

 
procedure dgemms (TRANSA,TRANSB: char;
                  m,n,k: int;
                  alpha: double;
                  A    : pdouble;
                  LDA  : int;
                  B    : pdouble;
                  LDB  : int;
                  beta : double;
                  C    : pdouble;
                  LDC  : int); cdecl; external 'libopenblas';



Aber vielleicht tut es die Win64 version
Bei mir läuft zum Beispiel ein Benchmark:
https://github.com/xianyi/OpenBLAS/blob ... Y/cgemm.py

Code: Alles auswählen

 
From: 128 To: 2048 Step=128 Loops=1
        SIZE                    Flops                                   Time
       128x128 :                 5837.307688 MFlops                 0.002874 sec
       256x256 :                 3653.882049 MFlops                 0.036733 sec
       384x384 :                 4101.778470 MFlops                 0.110436 sec
       512x512 :                 2420.294678 MFlops                 0.443641 sec
       640x640 :                 3680.661634 MFlops                 0.569776 sec
       768x768 :                 3381.641964 MFlops                 1.071633 sec
       896x896 :                 3778.228415 MFlops                 1.523091 sec
     1024x1024 :                 1854.743797 MFlops                 4.631332 sec
     1152x1152 :                 3248.714164 MFlops                 3.764748 sec
     1280x1280 :                 2466.665306 MFlops                 6.801578 sec
     1408x1408 :                 2684.847063 MFlops                 8.317224 sec
     1536x1536 :                 1857.497928 MFlops                15.607570 sec
     1664x1664 :                 2584.962270 MFlops                14.259219 sec
     1792x1792 :                 2267.826058 MFlops                20.299917 sec
     1920x1920 :                 2565.324269 MFlops                22.072494 sec

Gruß Horst

Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Dies ist noch eine interessante Seite, die zeigt, welch Erweiterung für was gut ist.
https://www.pc-erfahrung.de/hardware/pr ... -mehr.html

Es könnte sein, das AVX in FPC integriert wird:
http://forum.lazarus.freepascal.org/ind ... t.html#new
Zuletzt geändert von Mathias am Mo 25. Jun 2018, 18:36, insgesamt 1-mal geändert.
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Ich habe angefangen, ein kleines Tutorial für SSE zu schreiben.

http://wiki.freepascal.org/SSE/de#SSE
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Die Matritzen und Vektoren habe aber ganz andere Dimensionen , so in der Größe von 500*500 bis > 2000*2000 , oder noch mehr.


Es sind sogar Vektoren mit 8 Elementen möglich, könnte vielleicht für dich nützlich sein. Auf den neusten CPUS sind sogar 16 möglich, aber da macht der Assembler von FPC (noch) nicht mit.

Code: Alles auswählen

function Vec_Add256(const v0, v1: TVector8f): TVector8f; assembler;
asm
         Vmovups  ymm0, [v0]
         Vmovups  ymm1, [v1]
         Vmulps  ymm3, ymm1, ymm0
         Vmovups  [Result], ymm3
end
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Frank Ranis
Beiträge: 201
Registriert: Do 24. Jan 2013, 21:22

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Frank Ranis »

Hallo Mathias,

Mathias hat geschrieben:
Die Matritzen und Vektoren habe aber ganz andere Dimensionen , so in der Größe von 500*500 bis > 2000*2000 , oder noch mehr.


Es sind sogar Vektoren mit 8 Elementen möglich, könnte vielleicht für dich nützlich sein. Auf den neusten CPUS sind sogar 16 möglich, aber da macht der Assembler von FPC (noch) nicht mit.

Code: Alles auswählen

function Vec_Add256(const v0, v1: TVector8f): TVector8f; assembler;
asm
         Vmovups  ymm0, [v0]
         Vmovups  ymm1, [v1]
         Vmulps  ymm3, ymm1, ymm0
         Vmovups  [Result], ymm3
end


danke für deine Infos .

Habe noch mal darüber nachgedacht , selbst wenn man das für große Dimensionen hin bekommt hat immer noch das Problem mit den Prozessoren.
Bei Intel / AMD funzt es dann , aber wie sieht es aus bei Rechnern mit anderen Prozessoren und den asm-Befehlen ?
Dann müßte man doch wieder alles neu schreiben , oder denke ich da völlig falsch ?

Gruß

Frank
www.flz-vortex.de

Horst_h
Beiträge: 72
Registriert: Mi 20. Mär 2013, 08:57

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Horst_h »

Hallo,

dann muss man alles neu schreiben oder benutzt OpenBlas oder ähnliches, das für alle möglichen CPU´s und Plattformen zu kompilieren ist und von CPU/GPU-Herstellern hochoptimert wird und gibt die entsprechenden .dll/.so oder sonstwie mit.
https://de.wikipedia.org/wiki/Basic_Lin ... ubprograms
Dein Programm muss ja auch für die entsprechende Platform compiliert werden.Du kannst natürlich auch eine weitverbreitete Scriptsprache nehmen ( siehe das Phyton Beispiel ) und ( open- ) blas aufrufen.Die aufwendige Rechnerei wäre hochoptimiert, die Verwaltung und Ausgabe immer noch schnell genug.
Aber man muss das Programm auch internationalisieren, nicht nur CPU`s auch Menschen sprechen unterschiedliche Sprachen.Was war da noch mit Pfingsten... ;-)

Gruß Horst

Mathias
Beiträge: 6146
Registriert: Do 2. Jan 2014, 17:21
OS, Lazarus, FPC: Linux (die neusten Trunk)
CPU-Target: 64Bit
Wohnort: Schweiz

Re: Matrizen / Vectoren - Multiplikationen beschleunigen

Beitrag von Mathias »

Bei Intel / AMD funzt es dann , aber wie sieht es aus bei Rechnern mit anderen Prozessoren und den asm-Befehlen ?

Ich glaube kaum, das du diese komplexen Berechnungen auf einem Raspi machen willst.
Mit Lazarus sehe ich grün
Mit Java und C/C++ sehe ich rot

Antworten