Wie kann ich das theoretische Maximum von 4 FLOPs pro Zyklus erreichen?

Question

Wie kann die theoretische Spitzenleistung von 4 Gleitkommaoperationen (doppelte Genauigkeit) pro Zyklus auf einer modernen x86-64 Intel CPU erreicht werden?

Soweit ich weiß, dauert es drei Zyklen, bis ein SSE add und fünf Zyklen für eine mul auf den meisten modernen Intel-CPUs zu vervollständigen (siehe zum Beispiel Agner Fogs 'Anweisungstabellen' ). Durch das Pipelining kann man einen Durchsatz von einem add pro Zyklus, wenn der Algorithmus mindestens drei unabhängige Summierungen hat. Da dies sowohl für gepackte addpd sowie der Skalar addsd Versionen und SSE-Register können zwei double kann der Durchsatz bis zu zwei Flops pro Zyklus betragen.

Darüber hinaus scheint es (auch wenn ich keine richtige Dokumentation dazu gesehen habe) add und mul können parallel ausgeführt werden, was einen theoretischen maximalen Durchsatz von vier Flops pro Zyklus ergibt.

Es ist mir jedoch nicht gelungen, diese Leistung mit einem einfachen C/C++-Programm zu reproduzieren. Mein bester Versuch führte zu etwa 2,7 Flops/Zyklus. Wenn jemand ein einfaches C/C++- oder Assembler-Programm beisteuern kann, das die Spitzenleistung demonstriert, wäre ich sehr dankbar.

Mein Versuch:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

double stoptime(void) {
   struct timeval t;
   gettimeofday(&t,NULL);
   return (double) t.tv_sec + t.tv_usec/1000000.0;
}

double addmul(double add, double mul, int ops){
   // Need to initialise differently otherwise compiler might optimise away
   double sum1=0.1, sum2=-0.1, sum3=0.2, sum4=-0.2, sum5=0.0;
   double mul1=1.0, mul2= 1.1, mul3=1.2, mul4= 1.3, mul5=1.4;
   int loops=ops/10;          // We have 10 floating point operations inside the loop
   double expected = 5.0*add*loops + (sum1+sum2+sum3+sum4+sum5)
               + pow(mul,loops)*(mul1+mul2+mul3+mul4+mul5);

   for (int i=0; i<loops; i++) {
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
   }
   return  sum1+sum2+sum3+sum4+sum5+mul1+mul2+mul3+mul4+mul5 - expected;
}

int main(int argc, char** argv) {
   if (argc != 2) {
      printf("usage: %s <num>\n", argv[0]);
      printf("number of operations: <num> millions\n");
      exit(EXIT_FAILURE);
   }
   int n = atoi(argv[1]) * 1000000;
   if (n<=0)
       n=1000;

   double x = M_PI;
   double y = 1.0 + 1e-8;
   double t = stoptime();
   x = addmul(x, y, n);
   t = stoptime() - t;
   printf("addmul:\t %.3f s, %.3f Gflops, res=%f\n", t, (double)n/t/1e9, x);
   return EXIT_SUCCESS;
}

Zusammengestellt mit:

g++ -O2 -march=native addmul.cpp ; ./a.out 1000

erzeugt die folgende Ausgabe auf einem Intel Core i5-750, 2,66 GHz:

addmul:  0.270 s, 3.707 Gflops, res=1.326463

Das heißt, nur etwa 1,4 Flops pro Zyklus. Betrachtet man den Assembler-Code mit g++ -S -O2 -march=native -masm=intel addmul.cpp die Hauptschleife scheint ein wenig optimal zu sein.

.L4:
inc    eax
mulsd    xmm8, xmm3
mulsd    xmm7, xmm3
mulsd    xmm6, xmm3
mulsd    xmm5, xmm3
mulsd    xmm1, xmm3
addsd    xmm13, xmm2
addsd    xmm12, xmm2
addsd    xmm11, xmm2
addsd    xmm10, xmm2
addsd    xmm9, xmm2
cmp    eax, ebx
jne    .L4

Ersetzen der skalaren Versionen durch gepackte Versionen ( addpd y mulpd ) würde die Anzahl der Flops verdoppeln, ohne dass sich die Ausführungszeit ändert, und so käme ich auf knapp 2,8 Flops pro Zyklus. Gibt es ein einfaches Beispiel, bei dem vier Flops pro Zyklus erreicht werden?

Nettes kleines Programm von Mysticial; hier sind meine Ergebnisse (die allerdings nur ein paar Sekunden lang liefen):

• gcc -O2 -march=nocona : 5,6 Gflops von 10,66 Gflops (2,1 Flops/Zyklus)
• cl /O2 , openmp entfernt: 10,1 Gflops von 10,66 Gflops (3,8 Flops/Zyklus)

Das scheint alles etwas kompliziert zu sein, aber meine bisherigen Schlussfolgerungen:

• gcc -O2 ändert die Reihenfolge der unabhängigen Gleitkommaoperationen mit mit dem Ziel, abwechselnd addpd y mulpd wenn möglich. Das Gleiche gilt für gcc-4.6.2 -O2 -march=core2 .

• gcc -O2 -march=nocona scheint die Reihenfolge der Gleitkommaoperationen wie in dem C++-Quelltext definiert ist.

• cl /O2 der 64-Bit-Compiler aus dem SDK für Windows 7 rollt die Schleife automatisch ab und scheint zu versuchen, die Operationen zu arrangieren so anzuordnen, dass Gruppen von drei addpd wechseln sich ab mit drei mulpd (zumindest auf meinem System und für mein einfaches Programm).

• Meine Kern i5 750 ( Nehalem Architektur ) mag es nicht, wenn man abwechselnd addiert und muliert und scheint nicht in der Lage nicht in der Lage, beide Operationen parallel auszuführen. Wenn sie jedoch in 3er-Gruppen zusammengefasst wird, funktioniert sie plötzlich wie von Zauberhand.

• Andere Architekturen (möglicherweise Sandy Bridge und andere) scheinen in der Lage zu sein, add/mul ohne Probleme parallel auszuführen wenn sie sich im Assemblercode abwechseln.

• Auch wenn es schwer zuzugeben ist, aber auf meinem System cl /O2 leistet bei der Low-Level-Optimierung von Operationen für mein System eine viel bessere Arbeit und erreicht für das obige kleine C++-Beispiel fast die Spitzenleistung. Ich habe zwischen 1,85-2,01 Flops/Zyklus (ich habe clock() in Windows verwendet, was nicht so präzise ist. Ich denke, ich muss einen besseren Timer verwenden - danke Mackie Messer).

• Das Beste, was ich mit gcc war das manuelle Abrollen und Anordnen von Schleifen Additionen und Multiplikationen in Dreiergruppen. Mit g++ -O2 -march=nocona addmul_unroll.cpp Ich bekomme bestenfalls 0.207s, 4.825 Gflops was einem Wert von 1,8 Flops/Zyklus entspricht womit ich jetzt recht zufrieden bin.

Im C++-Code habe ich die for Schleife mit:

   for (int i=0; i<loops/3; i++) {
       mul1*=mul; mul2*=mul; mul3*=mul;
       sum1+=add; sum2+=add; sum3+=add;
       mul4*=mul; mul5*=mul; mul1*=mul;
       sum4+=add; sum5+=add; sum1+=add;

       mul2*=mul; mul3*=mul; mul4*=mul;
       sum2+=add; sum3+=add; sum4+=add;
       mul5*=mul; mul1*=mul; mul2*=mul;
       sum5+=add; sum1+=add; sum2+=add;

       mul3*=mul; mul4*=mul; mul5*=mul;
       sum3+=add; sum4+=add; sum5+=add;
   }

Die Baugruppe sieht nun wie folgt aus:

.L4:
mulsd    xmm8, xmm3
mulsd    xmm7, xmm3
mulsd    xmm6, xmm3
addsd    xmm13, xmm2
addsd    xmm12, xmm2
addsd    xmm11, xmm2
mulsd    xmm5, xmm3
mulsd    xmm1, xmm3
mulsd    xmm8, xmm3
addsd    xmm10, xmm2
addsd    xmm9, xmm2
addsd    xmm13, xmm2
...

Answer 1

Ich habe genau diese Aufgabe schon einmal erledigt. Dabei ging es aber hauptsächlich darum, den Stromverbrauch und die CPU-Temperaturen zu messen. Der folgende Code (der ziemlich lang ist) erreicht nahezu das Optimum auf meinem Core i7 2600K.

Das Wichtigste dabei ist die massive manuelle Schleifenentrollung sowie die Verschachtelung von Multiplizierungen und Additionen...

Das vollständige Projekt ist auf meinem GitHub zu finden: https://github.com/Mysticial/Flops

Warnung:

Wenn Sie sich entscheiden, dies zu kompilieren und auszuführen, achten Sie auf Ihre CPU-Temperaturen!!!
Achten Sie darauf, dass Sie es nicht überhitzen. Und stellen Sie sicher, dass die CPU-Drosselung Ihre Ergebnisse nicht beeinträchtigt!

Außerdem übernehme ich keine Verantwortung für Schäden, die durch die Ausführung dieses Codes entstehen könnten.

Anmerkungen:

• Dieser Code ist für x64 optimiert. x86 hat nicht genug Register, um ihn gut zu kompilieren.
• Dieser Code wurde auf Visual Studio 2010/2012 und GCC 4.6 getestet und funktioniert gut.
  ICC 11 (Intel Compiler 11) hat überraschenderweise Schwierigkeiten, es gut zu kompilieren.

• Diese sind für Vor-FMA-Prozessoren bestimmt. Um Spitzen-FLOPS auf Intel Haswell- und AMD Bulldozer-Prozessoren (und später) zu erreichen, werden FMA-Anweisungen (Fused Multiply Add) benötigt. Dies würde den Rahmen dieses Benchmarks sprengen.

  include <emmintrin.h>

  include <omp.h>

  include <iostream>

  using namespace std;

  typedef unsigned long long uint64;

  double test_dp_mac_SSE(double x,double y,uint64 iterations){ register __m128d r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,rA,rB,rC,rD,rE,rF;

  //  Generate starting data.
  r0 = _mm_set1_pd(x);
  r1 = _mm_set1_pd(y);
  
  r8 = _mm_set1_pd(-0.0);
  
  r2 = _mm_xor_pd(r0,r8);
  r3 = _mm_or_pd(r0,r8);
  r4 = _mm_andnot_pd(r8,r0);
  r5 = _mm_mul_pd(r1,_mm_set1_pd(0.37796447300922722721));
  r6 = _mm_mul_pd(r1,_mm_set1_pd(0.24253562503633297352));
  r7 = _mm_mul_pd(r1,_mm_set1_pd(4.1231056256176605498));
  r8 = _mm_add_pd(r0,_mm_set1_pd(0.37796447300922722721));
  r9 = _mm_add_pd(r1,_mm_set1_pd(0.24253562503633297352));
  rA = _mm_sub_pd(r0,_mm_set1_pd(4.1231056256176605498));
  rB = _mm_sub_pd(r1,_mm_set1_pd(4.1231056256176605498));
  
  rC = _mm_set1_pd(1.4142135623730950488);
  rD = _mm_set1_pd(1.7320508075688772935);
  rE = _mm_set1_pd(0.57735026918962576451);
  rF = _mm_set1_pd(0.70710678118654752440);
  
  uint64 iMASK = 0x800fffffffffffffull;
  __m128d MASK = _mm_set1_pd(*(double*)&iMASK);
  __m128d vONE = _mm_set1_pd(1.0);
  
  uint64 c = 0;
  while (c < iterations){
      size_t i = 0;
      while (i < 1000){
          //  Here's the meat - the part that really matters.
  
          r0 = _mm_mul_pd(r0,rC);
          r1 = _mm_add_pd(r1,rD);
          r2 = _mm_mul_pd(r2,rE);
          r3 = _mm_sub_pd(r3,rF);
          r4 = _mm_mul_pd(r4,rC);
          r5 = _mm_add_pd(r5,rD);
          r6 = _mm_mul_pd(r6,rE);
          r7 = _mm_sub_pd(r7,rF);
          r8 = _mm_mul_pd(r8,rC);
          r9 = _mm_add_pd(r9,rD);
          rA = _mm_mul_pd(rA,rE);
          rB = _mm_sub_pd(rB,rF);
  
          r0 = _mm_add_pd(r0,rF);
          r1 = _mm_mul_pd(r1,rE);
          r2 = _mm_sub_pd(r2,rD);
          r3 = _mm_mul_pd(r3,rC);
          r4 = _mm_add_pd(r4,rF);
          r5 = _mm_mul_pd(r5,rE);
          r6 = _mm_sub_pd(r6,rD);
          r7 = _mm_mul_pd(r7,rC);
          r8 = _mm_add_pd(r8,rF);
          r9 = _mm_mul_pd(r9,rE);
          rA = _mm_sub_pd(rA,rD);
          rB = _mm_mul_pd(rB,rC);
  
          r0 = _mm_mul_pd(r0,rC);
          r1 = _mm_add_pd(r1,rD);
          r2 = _mm_mul_pd(r2,rE);
          r3 = _mm_sub_pd(r3,rF);
          r4 = _mm_mul_pd(r4,rC);
          r5 = _mm_add_pd(r5,rD);
          r6 = _mm_mul_pd(r6,rE);
          r7 = _mm_sub_pd(r7,rF);
          r8 = _mm_mul_pd(r8,rC);
          r9 = _mm_add_pd(r9,rD);
          rA = _mm_mul_pd(rA,rE);
          rB = _mm_sub_pd(rB,rF);
  
          r0 = _mm_add_pd(r0,rF);
          r1 = _mm_mul_pd(r1,rE);
          r2 = _mm_sub_pd(r2,rD);
          r3 = _mm_mul_pd(r3,rC);
          r4 = _mm_add_pd(r4,rF);
          r5 = _mm_mul_pd(r5,rE);
          r6 = _mm_sub_pd(r6,rD);
          r7 = _mm_mul_pd(r7,rC);
          r8 = _mm_add_pd(r8,rF);
          r9 = _mm_mul_pd(r9,rE);
          rA = _mm_sub_pd(rA,rD);
          rB = _mm_mul_pd(rB,rC);
  
          i++;
      }
  
      //  Need to renormalize to prevent denormal/overflow.
      r0 = _mm_and_pd(r0,MASK);
      r1 = _mm_and_pd(r1,MASK);
      r2 = _mm_and_pd(r2,MASK);
      r3 = _mm_and_pd(r3,MASK);
      r4 = _mm_and_pd(r4,MASK);
      r5 = _mm_and_pd(r5,MASK);
      r6 = _mm_and_pd(r6,MASK);
      r7 = _mm_and_pd(r7,MASK);
      r8 = _mm_and_pd(r8,MASK);
      r9 = _mm_and_pd(r9,MASK);
      rA = _mm_and_pd(rA,MASK);
      rB = _mm_and_pd(rB,MASK);
      r0 = _mm_or_pd(r0,vONE);
      r1 = _mm_or_pd(r1,vONE);
      r2 = _mm_or_pd(r2,vONE);
      r3 = _mm_or_pd(r3,vONE);
      r4 = _mm_or_pd(r4,vONE);
      r5 = _mm_or_pd(r5,vONE);
      r6 = _mm_or_pd(r6,vONE);
      r7 = _mm_or_pd(r7,vONE);
      r8 = _mm_or_pd(r8,vONE);
      r9 = _mm_or_pd(r9,vONE);
      rA = _mm_or_pd(rA,vONE);
      rB = _mm_or_pd(rB,vONE);
  
      c++;
  }
  
  r0 = _mm_add_pd(r0,r1);
  r2 = _mm_add_pd(r2,r3);
  r4 = _mm_add_pd(r4,r5);
  r6 = _mm_add_pd(r6,r7);
  r8 = _mm_add_pd(r8,r9);
  rA = _mm_add_pd(rA,rB);
  
  r0 = _mm_add_pd(r0,r2);
  r4 = _mm_add_pd(r4,r6);
  r8 = _mm_add_pd(r8,rA);
  
  r0 = _mm_add_pd(r0,r4);
  r0 = _mm_add_pd(r0,r8);
  
  //  Prevent Dead Code Elimination
  double out = 0;
  __m128d temp = r0;
  out += ((double*)&temp)[0];
  out += ((double*)&temp)[1];
  
  return out;

  }

  void test_dp_mac_SSE(int tds,uint64 iterations){

  double *sum = (double*)malloc(tds * sizeof(double));
  double start = omp_get_wtime();

  pragma omp parallel num_threads(tds)

  {
      double ret = test_dp_mac_SSE(1.1,2.1,iterations);
      sum[omp_get_thread_num()] = ret;
  }
  
  double secs = omp_get_wtime() - start;
  uint64 ops = 48 * 1000 * iterations * tds * 2;
  cout << "Seconds = " << secs << endl;
  cout << "FP Ops  = " << ops << endl;
  cout << "FLOPs   = " << ops / secs << endl;
  
  double out = 0;
  int c = 0;
  while (c < tds){
      out += sum[c++];
  }
  
  cout << "sum = " << out << endl;
  cout << endl;
  
  free(sum);

  }

  int main(){ // (threads, iterations) test_dp_mac_SSE(8,10000000);

  system("pause");

  }

Ausgabe (1 Thread, 10000000 Iterationen) - Kompiliert mit Visual Studio 2010 SP1 - x64 Release:

Seconds = 55.5104
FP Ops  = 960000000000
FLOPs   = 1.7294e+010
sum = 2.22652

Der Rechner ist ein Core i7 2600K @ 4,4 GHz. Theoretische SSE-Spitze ist 4 Flops * 4,4 GHz = 17,6 GFlops . Mit diesem Code wird Folgendes erreicht 17,3 GFlops - nicht schlecht.

Ausgabe (8 Threads, 10000000 Iterationen) - Kompiliert mit Visual Studio 2010 SP1 - x64 Release:

Seconds = 117.202
FP Ops  = 7680000000000
FLOPs   = 6.55279e+010
sum = 17.8122

Die theoretische SSE-Spitze beträgt 4 Flops * 4 Kerne * 4,4 GHz = 70,4 GFlops. Tatsächlich ist 65,5 GFlops .

---

Gehen wir noch einen Schritt weiter. AVX...

#include <immintrin.h>
#include <omp.h>
#include <iostream>
using namespace std;

typedef unsigned long long uint64;

double test_dp_mac_AVX(double x,double y,uint64 iterations){
    register __m256d r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,rA,rB,rC,rD,rE,rF;

    //  Generate starting data.
    r0 = _mm256_set1_pd(x);
    r1 = _mm256_set1_pd(y);

    r8 = _mm256_set1_pd(-0.0);

    r2 = _mm256_xor_pd(r0,r8);
    r3 = _mm256_or_pd(r0,r8);
    r4 = _mm256_andnot_pd(r8,r0);
    r5 = _mm256_mul_pd(r1,_mm256_set1_pd(0.37796447300922722721));
    r6 = _mm256_mul_pd(r1,_mm256_set1_pd(0.24253562503633297352));
    r7 = _mm256_mul_pd(r1,_mm256_set1_pd(4.1231056256176605498));
    r8 = _mm256_add_pd(r0,_mm256_set1_pd(0.37796447300922722721));
    r9 = _mm256_add_pd(r1,_mm256_set1_pd(0.24253562503633297352));
    rA = _mm256_sub_pd(r0,_mm256_set1_pd(4.1231056256176605498));
    rB = _mm256_sub_pd(r1,_mm256_set1_pd(4.1231056256176605498));

    rC = _mm256_set1_pd(1.4142135623730950488);
    rD = _mm256_set1_pd(1.7320508075688772935);
    rE = _mm256_set1_pd(0.57735026918962576451);
    rF = _mm256_set1_pd(0.70710678118654752440);

    uint64 iMASK = 0x800fffffffffffffull;
    __m256d MASK = _mm256_set1_pd(*(double*)&iMASK);
    __m256d vONE = _mm256_set1_pd(1.0);

    uint64 c = 0;
    while (c < iterations){
        size_t i = 0;
        while (i < 1000){
            //  Here's the meat - the part that really matters.

            r0 = _mm256_mul_pd(r0,rC);
            r1 = _mm256_add_pd(r1,rD);
            r2 = _mm256_mul_pd(r2,rE);
            r3 = _mm256_sub_pd(r3,rF);
            r4 = _mm256_mul_pd(r4,rC);
            r5 = _mm256_add_pd(r5,rD);
            r6 = _mm256_mul_pd(r6,rE);
            r7 = _mm256_sub_pd(r7,rF);
            r8 = _mm256_mul_pd(r8,rC);
            r9 = _mm256_add_pd(r9,rD);
            rA = _mm256_mul_pd(rA,rE);
            rB = _mm256_sub_pd(rB,rF);

            r0 = _mm256_add_pd(r0,rF);
            r1 = _mm256_mul_pd(r1,rE);
            r2 = _mm256_sub_pd(r2,rD);
            r3 = _mm256_mul_pd(r3,rC);
            r4 = _mm256_add_pd(r4,rF);
            r5 = _mm256_mul_pd(r5,rE);
            r6 = _mm256_sub_pd(r6,rD);
            r7 = _mm256_mul_pd(r7,rC);
            r8 = _mm256_add_pd(r8,rF);
            r9 = _mm256_mul_pd(r9,rE);
            rA = _mm256_sub_pd(rA,rD);
            rB = _mm256_mul_pd(rB,rC);

            r0 = _mm256_mul_pd(r0,rC);
            r1 = _mm256_add_pd(r1,rD);
            r2 = _mm256_mul_pd(r2,rE);
            r3 = _mm256_sub_pd(r3,rF);
            r4 = _mm256_mul_pd(r4,rC);
            r5 = _mm256_add_pd(r5,rD);
            r6 = _mm256_mul_pd(r6,rE);
            r7 = _mm256_sub_pd(r7,rF);
            r8 = _mm256_mul_pd(r8,rC);
            r9 = _mm256_add_pd(r9,rD);
            rA = _mm256_mul_pd(rA,rE);
            rB = _mm256_sub_pd(rB,rF);

            r0 = _mm256_add_pd(r0,rF);
            r1 = _mm256_mul_pd(r1,rE);
            r2 = _mm256_sub_pd(r2,rD);
            r3 = _mm256_mul_pd(r3,rC);
            r4 = _mm256_add_pd(r4,rF);
            r5 = _mm256_mul_pd(r5,rE);
            r6 = _mm256_sub_pd(r6,rD);
            r7 = _mm256_mul_pd(r7,rC);
            r8 = _mm256_add_pd(r8,rF);
            r9 = _mm256_mul_pd(r9,rE);
            rA = _mm256_sub_pd(rA,rD);
            rB = _mm256_mul_pd(rB,rC);

            i++;
        }

        //  Need to renormalize to prevent denormal/overflow.
        r0 = _mm256_and_pd(r0,MASK);
        r1 = _mm256_and_pd(r1,MASK);
        r2 = _mm256_and_pd(r2,MASK);
        r3 = _mm256_and_pd(r3,MASK);
        r4 = _mm256_and_pd(r4,MASK);
        r5 = _mm256_and_pd(r5,MASK);
        r6 = _mm256_and_pd(r6,MASK);
        r7 = _mm256_and_pd(r7,MASK);
        r8 = _mm256_and_pd(r8,MASK);
        r9 = _mm256_and_pd(r9,MASK);
        rA = _mm256_and_pd(rA,MASK);
        rB = _mm256_and_pd(rB,MASK);
        r0 = _mm256_or_pd(r0,vONE);
        r1 = _mm256_or_pd(r1,vONE);
        r2 = _mm256_or_pd(r2,vONE);
        r3 = _mm256_or_pd(r3,vONE);
        r4 = _mm256_or_pd(r4,vONE);
        r5 = _mm256_or_pd(r5,vONE);
        r6 = _mm256_or_pd(r6,vONE);
        r7 = _mm256_or_pd(r7,vONE);
        r8 = _mm256_or_pd(r8,vONE);
        r9 = _mm256_or_pd(r9,vONE);
        rA = _mm256_or_pd(rA,vONE);
        rB = _mm256_or_pd(rB,vONE);

        c++;
    }

    r0 = _mm256_add_pd(r0,r1);
    r2 = _mm256_add_pd(r2,r3);
    r4 = _mm256_add_pd(r4,r5);
    r6 = _mm256_add_pd(r6,r7);
    r8 = _mm256_add_pd(r8,r9);
    rA = _mm256_add_pd(rA,rB);

    r0 = _mm256_add_pd(r0,r2);
    r4 = _mm256_add_pd(r4,r6);
    r8 = _mm256_add_pd(r8,rA);

    r0 = _mm256_add_pd(r0,r4);
    r0 = _mm256_add_pd(r0,r8);

    //  Prevent Dead Code Elimination
    double out = 0;
    __m256d temp = r0;
    out += ((double*)&temp)[0];
    out += ((double*)&temp)[1];
    out += ((double*)&temp)[2];
    out += ((double*)&temp)[3];

    return out;
}

void test_dp_mac_AVX(int tds,uint64 iterations){

    double *sum = (double*)malloc(tds * sizeof(double));
    double start = omp_get_wtime();

#pragma omp parallel num_threads(tds)
    {
        double ret = test_dp_mac_AVX(1.1,2.1,iterations);
        sum[omp_get_thread_num()] = ret;
    }

    double secs = omp_get_wtime() - start;
    uint64 ops = 48 * 1000 * iterations * tds * 4;
    cout << "Seconds = " << secs << endl;
    cout << "FP Ops  = " << ops << endl;
    cout << "FLOPs   = " << ops / secs << endl;

    double out = 0;
    int c = 0;
    while (c < tds){
        out += sum[c++];
    }

    cout << "sum = " << out << endl;
    cout << endl;

    free(sum);
}

int main(){
    //  (threads, iterations)
    test_dp_mac_AVX(8,10000000);

    system("pause");
}

Ausgabe (1 Thread, 10000000 Iterationen) - Kompiliert mit Visual Studio 2010 SP1 - x64 Release:

Seconds = 57.4679
FP Ops  = 1920000000000
FLOPs   = 3.34099e+010
sum = 4.45305

Die theoretische AVX-Spitze beträgt 8 Flops * 4,4 GHz = 35,2 GFlops . Tatsächlich ist 33,4 GFlops .

Ausgabe (8 Threads, 10000000 Iterationen) - Kompiliert mit Visual Studio 2010 SP1 - x64 Release:

Seconds = 111.119
FP Ops  = 15360000000000
FLOPs   = 1.3823e+011
sum = 35.6244

Die theoretische AVX-Spitze beträgt 8 Flops * 4 Kerne * 4,4 GHz = 140,8 GFlops. Tatsächlich ist 138,2 GFlops .

---

Nun zu einigen Erklärungen:

Der leistungsrelevante Teil sind natürlich die 48 Anweisungen innerhalb der inneren Schleife. Sie werden feststellen, dass sie in 4 Blöcke zu je 12 Anweisungen unterteilt ist. Jeder dieser 12 Anweisungsblöcke ist völlig unabhängig voneinander und benötigt im Durchschnitt 6 Zyklen zur Ausführung.

Es gibt also 12 Anweisungen und 6 Zyklen zwischen Ausgabe und Verwendung. Die Latenzzeit der Multiplikation beträgt 5 Zyklen, was gerade ausreicht, um Latenzstörungen zu vermeiden.

Der Normalisierungsschritt ist erforderlich, um zu verhindern, dass die Daten über- oder unterlaufen werden. Dies ist erforderlich, da der Nichtstun-Code den Umfang der Daten langsam erhöht/verringert.

Es ist also möglich, eine bessere Leistung zu erzielen, wenn man einfach alle Nullen verwendet und den Normalisierungsschritt weglässt. Da ich den Benchmark jedoch geschrieben habe, um den Stromverbrauch und die Temperatur zu messen, Ich musste sicherstellen, dass die Flops auf "echten" Daten und nicht auf Nullen beruhen. - da die Ausführungseinheiten sehr wohl über eine spezielle Behandlung von Nullen verfügen können, die weniger Strom verbraucht und weniger Wärme erzeugt.

---

Weitere Ergebnisse:

• Intel Core i7 920 @ 3,5 GHz
• Windows 7 Ultimate x64
• Visual Studio 2010 SP1 - x64-Version

Fäden: 1

Seconds = 72.1116
FP Ops  = 960000000000
FLOPs   = 1.33127e+010
sum = 2.22652

Theoretische SSE-Spitze: 4 Flops * 3,5 GHz = 14,0 GFlops . Tatsächlich ist 13,3 GFlops .

Fäden: 8

Seconds = 149.576
FP Ops  = 7680000000000
FLOPs   = 5.13452e+010
sum = 17.8122

Theoretische SSE-Spitze: 4 Flops * 4 Kerne * 3,5 GHz = 56,0 GFlops . Tatsächlich ist 51,3 GFlops .

Meine Prozessortemperaturen erreichten 76°C beim Multi-Thread-Lauf! Wenn Sie diese Tests durchführen, stellen Sie sicher, dass die Ergebnisse nicht durch CPU-Drosselung beeinflusst werden.

---

• 2 x Intel Xeon X5482 Harpertown @ 3,2 GHz
• Ubuntu Linux 10 x64
• GCC 4.5.2 x64 - (-O2 -msse3 -fopenmp)

Fäden: 1

Seconds = 78.3357
FP Ops  = 960000000000
FLOPs   = 1.22549e+10
sum = 2.22652

Theoretische SSE-Spitze: 4 Flops * 3,2 GHz = 12,8 GFlops . Tatsächlich ist 12,3 GFlops .

Fäden: 8

Seconds = 78.4733
FP Ops  = 7680000000000
FLOPs   = 9.78676e+10
sum = 17.8122

Theoretische SSE-Spitze: 4 Flops * 8 Kerne * 3,2 GHz = 102,4 GFlops . Tatsächlich ist 97,9 GFlops .

Answer 2

Es gibt einen Punkt in der Intel-Architektur, den die Leute oft vergessen: Die Dispatch-Ports werden von Int und FP/SIMD gemeinsam genutzt. Das bedeutet, dass Sie nur eine bestimmte Anzahl von FP/SIMD-Bursts erhalten, bevor die Schleifenlogik Blasen in Ihrem Fließkommastrom erzeugt. Mystical hat mehr Flops aus seinem Code herausgeholt, weil er längere Strides in seiner unrolled loop verwendet hat.

Wenn Sie sich die Architektur der Nehalem/Sandy Bridge ansehen http://www.realworldtech.com/page.cfm?ArticleID=RWT091810191937&p=6 ist es ganz klar, was passiert.

Im Gegensatz dazu sollte es bei AMD (Bulldozer) einfacher sein, Spitzenleistungen zu erreichen, da die INT- und FP/SIMD-Pipes über separate Ausgabeports mit eigenem Scheduler verfügen.

Dies ist nur theoretisch, da ich keinen dieser Prozessoren zum Testen habe.

Answer 3

Zweige können Sie definitiv davon abhalten, theoretische Spitzenleistungen zu erbringen. Erkennen Sie einen Unterschied, wenn Sie manuell einige Schleifen auflösen? Zum Beispiel, wenn Sie 5 oder 10 Mal so viele Operationen pro Schleifeniteration durchführen:

for(int i=0; i<loops/5; i++) {
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
      mul1*=mul; mul2*=mul; mul3*=mul; mul4*=mul; mul5*=mul;
      sum1+=add; sum2+=add; sum3+=add; sum4+=add; sum5+=add;
   }

Answer 4

Mit Intels icc Version 11.1 auf einem 2.4GHz Intel Core 2 Duo erhalte ich

Macintosh:~ mackie$ icc -O3 -mssse3 -oaddmul addmul.cc && ./addmul 1000
addmul:  0.105 s, 9.525 Gflops, res=0.000000
Macintosh:~ mackie$ icc -v
Version 11.1 

Das kommt dem Idealwert von 9,6 Gflops sehr nahe.

EDIT:

Ups, wenn man sich den Assemblercode ansieht, scheint es, dass icc nicht nur die Multiplikation vektorisiert, sondern auch die Additionen aus der Schleife herausgezogen hat. Durch das Erzwingen einer strengeren fp-Semantik ist der Code nicht mehr vektorisiert:

Macintosh:~ mackie$ icc -O3 -mssse3 -oaddmul addmul.cc -fp-model precise && ./addmul 1000
addmul:  0.516 s, 1.938 Gflops, res=1.326463

EDIT2:

Wie gewünscht:

Macintosh:~ mackie$ clang -O3 -mssse3 -oaddmul addmul.cc && ./addmul 1000
addmul:  0.209 s, 4.786 Gflops, res=1.326463
Macintosh:~ mackie$ clang -v
Apple clang version 3.0 (tags/Apple/clang-211.10.1) (based on LLVM 3.0svn)
Target: x86_64-apple-darwin11.2.0
Thread model: posix

Die innere Schleife des Clang-Codes sieht wie folgt aus:

        .align  4, 0x90
LBB2_4:                                 ## =>This Inner Loop Header: Depth=1
        addsd   %xmm2, %xmm3
        addsd   %xmm2, %xmm14
        addsd   %xmm2, %xmm5
        addsd   %xmm2, %xmm1
        addsd   %xmm2, %xmm4
        mulsd   %xmm2, %xmm0
        mulsd   %xmm2, %xmm6
        mulsd   %xmm2, %xmm7
        mulsd   %xmm2, %xmm11
        mulsd   %xmm2, %xmm13
        incl    %eax
        cmpl    %r14d, %eax
        jl      LBB2_4

EDIT3:

Zum Schluss noch zwei Vorschläge: Erstens, wenn Sie diese Art von Benchmarking mögen, sollten Sie die Verwendung des rdtsc Unterricht ist anstelle von gettimeofday(2) . Es ist viel genauer und liefert die Zeit in Zyklen, was Sie normalerweise sowieso interessiert. Für gcc und Freunde können Sie es so definieren:

#include <stdint.h>

static __inline__ uint64_t rdtsc(void)
{
        uint64_t rval;
        __asm__ volatile ("rdtsc" : "=A" (rval));
        return rval;
}

Zweitens sollten Sie Ihr Benchmark-Programm mehrere Male ausführen und die nur beste Leistung . In modernen Betriebssystemen laufen viele Dinge parallel, der Prozessor kann sich in einem Energiesparmodus mit niedriger Frequenz befinden usw. Wenn Sie das Programm wiederholt ausführen, erhalten Sie ein Ergebnis, das dem Idealfall näher kommt.